Лабораторная работа 4 / 2008-03-27-19-33-Ольга-4
.docFilter Design
Пакет Filter Design содержит функции, предназначенные для синтеза дискретных фильтров.
Выберем тип АЧХ – фильтр нижних частот Lowpass.
Выберем нерекурсивный фильтр Equiripple – синтез фильтров с равномерными пульсациями АЧХ методом Ремеза.
Порядок фильтра зададим 20 – число коэффициентов фильтра.
Fs – частота дискретизации.
Fpass – конец полосы пропускания.
Fstop - начало полосы задерживания.
Выбрав опцию Magnitude and Phase Response, строим АЧХ и ФЧХ.
Также можно построить реакцию на единичный импульс.
Реакция на ступенчатую функцию.
Экспортируем коэффициенты фильтра в рабочее пространство. Получаем вектор Num, размером 1*21 (21 отсчет в характеристике фильтра).
Num =
Columns 1 through 8
0.0378 0.0141 -0.0034 -0.0294 -0.0488 -0.0437 -0.0030 0.0694
Columns 9 through 16
0.1538 0.2215 0.2474 0.2215 0.1538 0.0694 -0.0030 -0.0437
Columns 17 through 21
-0.0488 -0.0294 -0.0034 0.0141 0.0378
Fs = 1000;%частота дискретизации
POINTS = 128; % число точек в частотной характеристике фильтра
f =freqz(Num,1,POINTS);
A = abs(f);
B = unwrap(angle(f)); % замена фазы f большей pi на ее дополнение до 2*pi
figure(1);
subplot(2,1,1), semilogx (1:POINTS,A),grid ;
subplot(2,1,2), semilogx (1:POINTS,B),grid;
Функция CHIRP предназначена для генерации колебаний с единичной амплитудой, мгновенная частота которых меняется по закону (линейно, квадратично, логарифмический закон).
Математический закон изменения мгновенной частоты:
F(t) = Fo + bt , b = (F1-Fo)/t1 ; 'linear'
F(t) = Fo + bt^2 , b = (F1-Fo)/t1^2 ; 'quadratic'
F(t) = Fo + exp(bt), b = ln(F1-Fo)/t1; 'logarithmic' - противоречит названию, т.к. зависимость частоты о времени экспоненциальная.
f1 = 100;%начальная частота
f2 = 1000;%конечная частота
FS=5000;
dt = 1/FS;
t = 0:dt:1-dt;
N=length(t);
x = chirp(t,f1,t(N),f2,'linear'),title('linear'); % линейно
q = chirp(t,f1,t(N),f2,'quadratic'),title('quadratic'); % квадратично
w = chirp(t,f1,t(N),f2,'logarithmic'),title('logarithmic'); % логарифмический закон
figure(3);
specgram(x,[],Fs),title('linear');
figure(9);
specgram(q,[],Fs),title('quadratic');
figure(10);
specgram(w,[],Fs),title('logarithmic');
Применим построенный фильтр к колебаниям с меняющейся частотой. Для этого воспользуемся функцией FILTER, фильтрация значений в соответствии с вектором коэффициентов Num.
dt = 1/Fs;
t = 0:dt:0.5-dt;
N = length(t);
x = chirp(t,80,t(N),250,'linear');
q = chirp(t,80,t(N),250,'quadratic');
w = chirp(t,80,t(N),250,'logarithmic');
figure(2);
subplot(2,1,1),plot(x) ;
y = filter(Num,1,x);
subplot(2,1,2),plot(y) ;
Для сигнала с линейно меняющейся частотой.
figure(11);
subplot(2,1,1),plot(q) ;
y = filter(Num,1,q);
subplot(2,1,2),plot(y) ;
Для сигнала с квадратично меняющейся частотой.
figure(12);
subplot(2,1,1),plot(w) ;
y = filter(Num,1,w);
subplot(2,1,2),plot(y) ;
Для сигнала с частотой, меняющейся по логарифмическому закону.
Поскольку в первом случае с линейным методом частота меняется быстрее в зависимости от времени, то профильтрованный сигнал быстрее затухает, т.к. фильтр не пропускает высокие частоты. Частота, меняющаяся по логарифмическому закону, растет медленнее от времени, поэтому фильтрованный сигнал медленнее затухает, чем в двух предыдущих случаях.
Создадим два сигнала с разной частотой, сложим их и применим фильтр к полученному сигналу.
F1 = 90;
F2 = 200;
x1 = sin(2*pi*F1*t);
x2 = sin(2*pi*F2*t);
x = x1+x2;
y = filter(Num,1,x);
subplot(2,1,1),plot(x1),title('x1');
subplot(2,1,2),plot(x2),title('x2');
Сигналы x1 и х2.
subplot(2,1,1),plot(x),title('x');
subplot(2,1,2),plot(y),title('filter');
Сумма двух сигналов и сигнал, профильтрованный ФНЧ. После фильтрации осталась только низкочастотная составляющая сигнала.
Задержим сигнал х на 50 отсчетов и применим фильтр.
start = 50;
z = zeros(1,N);
z(start:N) = x(start:N);
y = filter(Num,1,z);
subplot(2,1,1),plot(z),grid,title('delay signal') ;
subplot(2,1,2),plot(y),grid,title('filtered signal');
plot(y); grid; title('filtered signal'), xlim([30 100]);
Сигнал, пропущенный через фильтр, задержан относительно исходного сигнала.
Это связано с фазовой и групповой задержкой в фильтре.
Фильтр FIR1.
n = 40; % порядок фильтра
b = fir1 (n,[50/Fs 250/Fs]); % возвращает вектор коэффициентов длиной n+1=41
% создание полосового фильтра методом окон, в данном случае тип окна не указан, поэтому по умолчанию используется окно Хэмминга.
Частотные характеристики полосового фильтра.
W = freqz(b);
subplot(3,1,1),stem(b),grid,title('impulse response');
subplot(3,1,2),plot(abs(W)),grid,title('magnitude response');
subplot(3,1,3),plot(unwrap(angle(W))),grid,title('phase response');
Реакция на импульсное воздействие, АЧХ и ФЧХ.
В полосе пропускания фазовая характеристика линейна.
Медианные фильтры.
Создадим сигнал, один из участков которого меняется по линейному и синусоидальному закону, и отдельные импульсы, которые нужно профильтровать.
z = zeros(1,N);
u=ones(1,30); % единичный вектор
z(10) = 5;
z(2) = 3;
z(21:30) = 1:10;
z(40) = 7;
z(60) = 6;
y = medfilt1(z,7); % порядок фильтра N=7
% Y(k) это медиана на Z ( k-(N-1)/2 : k+(N-1)/2 ), т.е. Z( k-3 : k+3 )
subplot(2,1,1),stem(z), xlim([0 100]), grid;
subplot(2,1,2),stem(y),grid, xlim([0 100]), title('filtered signal');
На выходе фильтра сигнал сглажен, т.к. число 7 4 раза является медианой на промежутках Z( k-3 : k+3 ): 4,5,6,7,8,9,10; 0,5,6,7,8,9,10; 0,0,6,7,8,9,10; 0,0,0,7,8,9,10, (промежутки упорядочены по возрастанию).
figure(12);
z(21:30) = 5*x1(21:30); % синусоидальный участок сигнала
y = medfilt1(z,7);
subplot(2,1,1), stem(z), xlim([0 100]),grid;
subplot(2,1,2), stem(y), grid, xlim([0 100]), title('filtered signal');
В данном случае число 2.13 является медианой для каждого участка отсчетов синусоиды и синусоидальный сигнал полностью сглаживается.
figure(13);
z(21:30) = 5*х2(21:30); % синусоидальный участок сигнала
y = medfilt1(z,7);
subplot(2,1,1), stem(z), xlim([0 100]), grid;
subplot(2,1,2), stem(y), grid, xlim([0 100]), title('filtered signal');
В данном случае медианами являются 2,94 ; 0; -2,94.
Одиночные импульсы после фильтрования пропадают, поскольку на тех промежутках где они находятся, Z(k-3 : k+3 ), при построении вариационного ряда такие импульсы всегда стоят в конце, т.к. вокруг нулевые значения и никогда не являются медианой ряда.
Медианный фильтр обеспечивает удаление случайных одиночный импульсов.
Рекурсивный фильтр Butterworth.
Wn=0.15; % частота среза
Nr=20; % порядок фильтра
fp=1:500; % отсчеты частот
[b,a] = BUTTER(Nr,Wn); % вектора a, b коэффициентов фильтра
imp=impz(b,a); % импульсная характеристика
step=stepz(b,a); % переходная характеристика
fr=freqz(b,a,fp,Fs); % расчет частотных характеристик фильтра
A = abs(fr);
B = unwrap(angle(fr));
figure(14);
subplot(2,1,1),semilogx (1:500,A),grid ;
subplot(2,1,2),semilogx (1:500,B),grid;
Амплитудно-частотная и фазочастотная характеристики
figure(16);
subplot(2,1,1),plot(imp),title('impulse response'),grid;
subplot(2,1,2),plot(step),title('step response'),grid;
Реакция на импульсную функцию и на единичное ступенчатое воздействие.
Групповая и фазовая задержки фильтра.
tau=grpdelay(b,a,fp,Fs)/Fs; % вычисление групповой задержки
tau2=phasedelay(b,a,fp,Fs)/(2*pi); % вычисление фазовой задержки
figure(17);
plot(fp,tau,'b--',fp,tau2,'b'),grid;
xlabel('freq, Hz'); ylabel('delay, s');
title('group delay and phase delay');
legend('group delay','phase delay');
Фазовая задержка – это задержка несущей частоты, групповая – огибающей.
Пропускаем сигнал х через фильтр.
F1 = 90; % огибающая
F2 = 150; % несущая частота
x1 = sin(2*pi*F1*t);
x2 = sin(2*pi*F2*t);
x = x1+x2;
figure(15);
y = filter(b,a,x);
subplot(2,1,1),plot(x),title('x') ;
subplot(2,1,2),plot(y),title('filter') ;
Сдвиг огибающей относительно несущей частоты обусловлен различием фазовой и групповой задержек на этой частоте. На несущей частоте сигнала (200) эти задержки существенно различаются.