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 оборотов в день. Используйте эту информацию в спецификациях для диапазонов остановок фильтра.

Проектируйте конечную импульсную характеристику равновесия и БИХ Фильтров Баттерворта со следующими спецификациями: полоса пропускания от [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 конечную импульсную характеристику фильтр для сравнения. Спецификации фильтра 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 не показывает задержки из-за фазового отклика конечной импульсной характеристики.