Документация

КИХ-фильтр Lowpass – метод окна

Этот пример показывает как разработке и реализации КИХ-фильтр с помощью двух функций командной строки, fir1 и designfilt и интерактивного приложения Filter Designer.

Создайте сигнал использовать в примерах. Сигнал является синусоидой на 100 Гц в дополнении $$N(0,1/4)$$ белый Гауссов шум. Установите генератор случайных чисел на состояние по умолчанию для восстанавливаемых результатов.

rng default

Fs = 1000;
t = linspace(0,1,Fs);
x = cos(2*pi*100*t)+0.5*randn(size(t));

Проект фильтра является КИХ lowpass фильтр с порядком, равным 20 и частота среза 150 Гц. Используйте окно Kaiser с длиной одна выборка, больше, чем порядок фильтра и $$\beta =3$$. Смотрите kaiser для получения дополнительной информации об окне Kaiser.

Используйте fir1, чтобы разработать фильтр. fir1 требует нормированных частот в интервале [0,1], где 1 соответствует $$\pi$$ рад/выборка. Чтобы использовать fir1, необходимо преобразовать все спецификации частоты в нормированные частоты.

Разработайте фильтр и просмотрите ответ значения фильтра.

fc = 150;
Wn = (2/Fs)*fc;
b = fir1(20,Wn,'low',kaiser(21,3));

fvtool(b,1,'Fs',Fs)

Примените фильтр к сигналу и постройте результат в течение первых десяти периодов синусоиды на 100 Гц.

y = filter(b,1,x);

plot(t,x,t,y)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Разработайте тот же фильтр с помощью designfilt. Установите ответ фильтра на 'lowpassfir' и введите спецификации как пары Name,Value. С designfilt можно задать проект фильтра в Гц.

Fs = 1000;
Hd = designfilt('lowpassfir','FilterOrder',20,'CutoffFrequency',150, ...
       'DesignMethod','window','Window',{@kaiser,3},'SampleRate',Fs);

Отфильтруйте данные и постройте результат.

y1 = filter(Hd,x);

plot(t,x,t,y1)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

КИХ-фильтр Lowpass с Filter Designer

Этот пример показывает как разработке и реализации lowpass КИХ-фильтр с помощью метода окна с интерактивным приложением Filter Designer.

y2 = filter(Hd,x);

plot(t,x,t,y2)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Можно также использовать интерактивный инструмент filterBuilder, чтобы разработать фильтр.

Полосовые фильтры – КИХ минимального заказа и БИХ-системы

Этот пример показывает, как разработать данные о полосовом фильтре и фильтре с КИХ минимального заказа equiripple и БИХ-Фильтрами Баттерворта. Можно смоделировать много реальных сигналов как суперпозицию колеблющихся компонентов, низкочастотного тренда и аддитивного шума. Например, экономические данные часто содержат колебания, которые представляют циклы, наложенные на медленное варьирование вверх или тенденцию к понижению. Кроме того, существует аддитивный шумовой компонент, который является комбинацией погрешности измерения и свойственных случайных колебаний процесса.

В этих примерах примите вас демонстрационный некоторый процесс каждый день в течение одного года. Примите, что процесс имеет колебания приблизительно в однонедельных и одномесячных шкалах. Кроме того, существует низкочастотный восходящий тренд в данных и дополнении $$N(0,1/4)$$ белый Гауссов шум.

Создайте сигнал как суперпозицию двух синусоид с частотами 1/7 и 1/30 циклов/день. Добавьте низкочастотный увеличивающийся термин тренда и $$N(0,1/4)$$ белый Гауссов шум. Сбросьте генератор случайных чисел для восстанавливаемых результатов. Данные выбираются на 1 выборке/день. Постройте получившийся сигнал и оценку степени спектральной плотности (PSD).

rng default

Fs = 1;
n = 1:365;

x = cos(2*pi*(1/7)*n)+cos(2*pi*(1/30)*n-pi/4);
trend = 3*sin(2*pi*(1/1480)*n);

y = x+trend+0.5*randn(size(n));

[pxx,f] = periodogram(y,[],[],Fs);

subplot(2,1,1)
plot(n,y)
xlim([1 365])
xlabel('Days')
grid

subplot(2,1,2)
plot(f,10*log10(pxx))
xlabel('Cycles/day')
ylabel('dB')
grid

Низкочастотный тренд появляется в степени спектральная оценка плотности как увеличенная низкочастотная степень. Низкочастотная степень появляется на приблизительно 10 дБ выше колебания в 1/30 циклах/день. Используйте эту информацию в спецификациях для полос задерживания фильтра.

Разработайте КИХ минимального заказа equiripple и БИХ-Фильтры Баттерворта со следующими спецификациями: полоса пропускания от [1/40,1/4] циклов/день и полосы задерживания от [0,1/60] и [1/4,1/2] циклов/день. Установите и затухания полосы задерживания на 10 дБ и допуск пульсации полосы пропускания к 1 дБ.

Hd1 = designfilt('bandpassfir', ...
    'StopbandFrequency1',1/60,'PassbandFrequency1',1/40, ...
    'PassbandFrequency2',1/4 ,'StopbandFrequency2',1/2 , ...
    'StopbandAttenuation1',10,'PassbandRipple',1, ...
    'StopbandAttenuation2',10,'DesignMethod','equiripple','SampleRate',Fs);
Hd2 = designfilt('bandpassiir', ...
    'StopbandFrequency1',1/60,'PassbandFrequency1',1/40, ...
    'PassbandFrequency2',1/4 ,'StopbandFrequency2',1/2 , ...
    'StopbandAttenuation1',10,'PassbandRipple',1, ...
    'StopbandAttenuation2',10,'DesignMethod','butter','SampleRate',Fs);

Сравните порядок КИХ и БИХ-фильтров и развернутых фазовых откликов.

fprintf('The order of the FIR filter is %d\n',filtord(Hd1))

The order of the FIR filter is 78

fprintf('The order of the IIR filter is %d\n',filtord(Hd2))

The order of the IIR filter is 8

[phifir,w] = phasez(Hd1,[],1);
[phiiir,w] = phasez(Hd2,[],1);

figure
plot(w,unwrap(phifir))
hold on
plot(w,unwrap(phiiir))
hold off

xlabel('Cycles/Day')
ylabel('Radians')
legend('FIR Equiripple Filter','IIR Butterworth Filter')
grid

БИХ-фильтр имеет очень более низкоуровневый что КИХ-фильтр. Однако КИХ-фильтр имеет линейный фазовый отклик по полосе пропускания, в то время как БИХ-фильтр не делает. КИХ-фильтр задерживает все частоты в полосе пропускания фильтра одинаково, в то время как БИХ-фильтр не делает.

Кроме того, скорость изменения фазы на единицу частоты больше в КИХ-фильтре, чем в БИХ-фильтре.

Разработайте lowpass КИХ equiripple фильтр для сравнения. Спецификации фильтра lowpass: полоса пропускания [0,1/4] циклы/день, затухание полосы задерживания равняется 10 дБ, и набор допуска пульсации полосы пропускания к 1 дБ.

Hdlow = designfilt('lowpassfir', ...
    'PassbandFrequency',1/4,'StopbandFrequency',1/2, ...
    'PassbandRipple',1,'StopbandAttenuation',10, ...
    'DesignMethod','equiripple','SampleRate',1);

Отфильтруйте данные с фильтрами lowpass и полосой пропускания.

yfir = filter(Hd1,y);
yiir = filter(Hd2,y);
ylow = filter(Hdlow,y);

Постройте оценку PSD полосового БИХ-фильтра вывод. Можно заменить yiir на yfir в следующем коде, чтобы просмотреть оценку PSD КИХ-полосового фильтра вывод.

[pxx,f] = periodogram(yiir,[],[],Fs);

plot(f,10*log10(pxx))

xlabel('Cycles/day')
ylabel('dB')
grid

Оценка PSD показывает, что полосовой фильтр ослабляет низкочастотный тренд и высокочастотный шум.

Постройте первые 120 дней КИХ и БИХ-фильтра вывод.

plot(n,yfir,n,yiir)

axis([1 120 -2.8 2.8])
xlabel('Days')
legend('FIR bandpass filter output','IIR bandpass filter output', ...
    'Location','SouthEast')

Увеличенная задержка фазы КИХ-фильтра очевидна в фильтре вывод.

Постройте lowpass КИХ-фильтр вывод, наложенный на суперпозицию 7-дневных и 30-дневных циклов для сравнения.

plot(n,x,n,ylow)

xlim([1 365])
xlabel('Days')
legend('7-day and 30-day cycles','FIR lowpass filter output', ...
    'Location','NorthWest')

Вы видите в предыдущем графике, что низкочастотный тренд очевиден в фильтре lowpass вывод. В то время как фильтр lowpass сохраняет 7-дневные и 30-дневные циклы, полосовые фильтры выполняют лучше в этом примере, потому что полосовые фильтры также удаляют низкочастотный тренд.

Фильтрация нулевой фазы

Этот пример показывает, как выполнить фильтрацию нулевой фазы.

Повторите генерацию сигнала и проект фильтра lowpass с fir1 и designfilt. Вы не должны выполнять следующий код, если у вас уже есть эти переменные в вашей рабочей области.

rng default

Fs = 1000;
t = linspace(0,1,Fs);
x = cos(2*pi*100*t)+0.5*randn(size(t));

% Using fir1
fc = 150;
Wn = (2/Fs)*fc;
b = fir1(20,Wn,'low',kaiser(21,3));

% Using designfilt
Hd = designfilt('lowpassfir','FilterOrder',20,'CutoffFrequency',150, ...
       'DesignMethod','window','Window',{@kaiser,3},'SampleRate',Fs);

Отфильтруйте данные с помощью filter. Постройте первые 100 точек фильтра вывод наряду с наложенной синусоидой с той же амплитудной и начальной фазой как входной сигнал.

yout = filter(Hd,x);
xin = cos(2*pi*100*t);

plot(t,xin,t,yout)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Input Sine Wave','Filtered Data')
grid

Смотря на начальные 0,01 секунды отфильтрованных данных, вы видите, что вывод задерживается относительно входа. Задержка, кажется, составляет приблизительно 0,01 секунды, который является почти 1/2, длина КИХ просачиваются выборки $$(\mathrm{}10\times 0.\mathrm{}\mathrm{}001)$$.

Эта задержка происходит из-за фазового отклика фильтра. КИХ просачивается, этими примерами является фильтр линейной фазы типа I. Групповая задержка фильтра является 10 выборками.

Постройте групповую задержку с помощью fvtool.

fvtool(Hd,'Analysis','grpdelay')

Во многих приложениях искажение фазы приемлемо. Это особенно верно, когда фазовый отклик линеен. В других приложениях желательно иметь фильтр с нулевым фазовым откликом. Нулевой фазовый отклик не находится технически возможно в непричинном фильтре. Однако можно реализовать фильтрацию нулевой фазы с помощью причинного фильтра с filtfilt.

Отфильтруйте входной сигнал с помощью filtfilt. Постройте ответы, чтобы сравнить фильтр выходные параметры, полученные с filter и filtfilt.

yzp = filtfilt(Hd,x);

plot(t,xin,t,yout,t,yzp)

xlim([0 0.1])
xlabel('Time (s)')
ylabel('Amplitude')
legend('100-Hz Sine Wave','Filtered Signal','Zero-phase Filtering',...
    'Location','NorthEast')

В предыдущей фигуре вы видите, что вывод filtfilt не показывает задержку из-за фазового отклика КИХ-фильтра.