КИХ-фильтр 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 Гц. Используйте окно Кайзера с длиной одна выборка, больше, чем порядок фильтра и $$\beta =3$$. Смотрите 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, длина КИХ просачиваются выборки $$(10\times 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 не показывает задержку из-за фазового отклика КИХ-фильтра.