Фильтр FIR нижних частот - метод окна

В этом примере показано, как спроектировать и реализовать фильтр FIR с помощью двух функций командной строки: fir1 и designfiltи интерактивное приложение Filter Designer.

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

rng default

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

Конструкция фильтра представляет собой фильтр нижних частот КИХ с порядком 20 и частотой отсечки 150 Гц. Используйте окно Кайзера длиной на один образец больше порядка фильтра и $$\beta \mathrm{=}$$3. Посмотритеkaiser для получения подробной информации об окне Кайзера.

Использовать fir1 для проектирования фильтра. fir1 требует нормализованных частот в интервале [0,1], где 1 соответствует $$\delta$$ рад/выборка. Использовать 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')

Фильтр FIR нижних частот с конструктором фильтров

В этом примере показано, как разработать и реализовать фильтр FIR нижних частот с помощью метода окна в интерактивном приложении 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 для проектирования фильтра.

Полосовые фильтры - системы FIR и IIR минимального порядка

В этом примере показано, как проектировать полосовой фильтр и фильтровать данные с помощью эквипплов FIR минимального порядка и фильтров BIR Butterworth. Можно моделировать множество реальных сигналов как суперпозицию колеблющихся компонентов, низкочастотный тренд и аддитивный шум. Например, экономические данные часто содержат колебания, которые представляют собой циклы, наложенные на медленно изменяющийся восходящий или нисходящий тренд. Кроме того, существует аддитивная шумовая составляющая, которая представляет собой комбинацию погрешности измерения и собственных случайных флуктуаций в процессе.

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

Создайте сигнал как суперпозицию двух синусоидальных волн с частотами 1/7 и 1/30 циклов/день. Добавьте член низкочастотного тренда увеличения и $$N\mathrm{(}\mathrm{}\mathrm{}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 циклов/день. Используйте эту информацию в спецификациях для стоп-полос фильтра.

Проектирование эквриппельных фильтров минимального порядка и BIR Butterworth со следующими характеристиками: полоса пропускания от [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);

Сравните порядок фильтров FIR и IIR и отклики развернутой фазы.

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

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

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

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

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

Фильтрация данных с помощью полосового и низкочастотного фильтров.

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

Постройте график оценки PSD выхода полосового БИХ-фильтра. Вы можете заменить yiir с yfir в следующем коде для просмотра оценки PSD выходного сигнала полосового фильтра FIR.

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

plot(f,10*log10(pxx))

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

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

Постройте график первых 120 дней выхода фильтра FIR и IIR.

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')

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

Постройте график выходного сигнала фильтра КИХ нижних частот, наложенного на наложение 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')

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

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

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

Повторите схему формирования сигнала и фильтра нижних частот с помощью 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 длину FIR-фильтра в выборках $$(\mathrm{}10\mathrm{\times }\mathrm{}\mathrm{}\mathrm{}$$0,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 не проявляет задержки из-за фазовой характеристики фильтра FIR.