Добавил:
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

Лабораторная работа 4 / 2008-03-27-19-33-Ольга-4

.doc
Скачиваний:
32
Добавлен:
01.05.2014
Размер:
876.54 Кб
Скачать

Filter 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) эти задержки существенно различаются.