Проектирование КИХ – фильтра нижних частот с помощью оконного метода
Фильтр с конечной импульсной характеристикой (нерекурсивный фильтр, трансверсальный фильтр, КИХ-фильтр) или FIR-фильтр (FIR сокр. от finite impulse response — конечная импульсная характеристика) — один из видов линейных фильтров, характерной особенностью которого является ограниченность по времени его импульсной характеристики (с какого-то момента времени его выходной сигнал становится точно равным нулю). Такой фильтр называют часто нерекурсивным из-за отсутствия обратной связи. Знаменатель передаточной функции такого фильтра — константа. Может быть реализован как в цифровом виде (аппаратном или программном), так и в аналоговой структуре.
В этом примере показано, как спроектировать и реализовать КИХ-фильтр с помощью двух функций командной строки, 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));
Конструкция фильтра представляет собой КИХ-фильтр нижних частот с порядком, равным 20, и частотой среза 150 Гц. Используйте окно Кайзера с длиной на одну выборку (sample) больше порядка фильтра и β=3.
Окно Кайзера описывается формулой:
где I0 - модифицированная функция Бесселя первого рода нулевого порядка; β -коэффициент определяющий долю энергии, сосредоточенной в главном лепестке спектра оконной функции. Чем больше β тем больше доля энергии, и шире главный лепесток, и меньше уровень боковых лепестков. На практике используются значения от 4 до 9.
Используйте fir1 для создания фильтра.
fir1 требует нормированных частот в интервале [0,1], где 1 соответствует π рад/выборка. Чтобы использовать fir1, необходимо преобразовать все характеристики частот в нормированные частоты.
Создайте фильтр и просмотрите на амплитудную характеристику.
fc = 150;
Wn = (2/Fs)*fc;
b = fir1(20,Wn,"low",kaiser(21,3));
[h,f] = freqz(b,1,[],Fs);
plot(f,mag2db(abs(h)))
xlabel("Frequency (Hz)")
ylabel("Magnitude (dB)")
grid
Примените фильтр к сигналу и постройте график для первых десяти периодов синусоиды частотой 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")
Проектирование низкочастотного КИХ-фильтра с помощью конструктора фильтров
В этом примере показано, как спроектировать и реализовать КИХ-фильтр нижних частот с помощью интерактивного приложения Filter Designer.
- Запустите приложение, введя filterDesigner в командной строке.
- Установите для Response Type (тип характеристики) значение Lowpass (низкочастотный).
- Установите Design Method (метод проектирования) на FIR (КИХ)и выберите метод Window (метод окна).
- В разделе Filter Order (порядок фильтра) выберите Specify order (указать порядок). Задайте порядок 20.
- В разделе Frequency Specifications (частотные характеристики) установите Units (Единицы измерения ) на Hz (Гц), Fs (частота дискредитации) на 1000 и Fc (частота среза) на 150.
- Щелкните Design Filter.
- Выберите File > Export..., чтобы экспортировать КИХ-фильтр в рабочую область MATLAB в виде коэффициентов или объекта фильтра. В этом примере экспортируйте фильтр как объект. Укажите имя переменной как Hd.
- Нажмите " Export ".
- Отфильтруйте входной сигнал в командном окне с помощью экспортированного объекта фильтра. Постройте график для первых десяти периодов синусоиды частотой 100 Гц.
y2 = filter(Hd,x);
plot(t,x,t,y2)
xlim([0 0.1])
xlabel("Time (s)")
ylabel("Amplitude")
legend("Original Signal","Filtered Data")
- Выберите File > Generate MATLAB Code > Filter Design Function , чтобы сгенерировать функцию MATLAB для создания объекта filter с использованием ваших спецификаций.
Вы также можете использовать интерактивный инструмент filterBuilder для создания своего фильтра.
Полосовые фильтры – Системы КИХ и БИХ минимального порядка
В этом примере показано, как спроектировать полосно-пропускающий фильтр и отфильтровать данные с помощью КИХ-фильтров Баттерворта минимального порядка. Многие реальные сигналы можно смоделировать как суперпозицию колебательных компонентов, низкочастотного тренда и аддитивного шума. Например, экономические данные часто содержат колебания, которые представляют собой циклы, накладывающиеся на медленно меняющийся восходящий или нисходящий тренд. Кроме того, существует аддитивный шумовой компонент, который представляет собой сочетание погрешности измерения и присущих процессу случайных колебаний.
Фильтр с бесконечной импульсной характеристикой (Рекурсивный фильтр, БИХ-фильтр) или IIR-фильтр (IIR сокр. от infinite impulse response — бесконечная импульсная характеристика) — линейный электронный фильтр, использующий один или более своих выходов в качестве входа, то есть образующий обратную связь. Основным свойством таких фильтров является то, что их импульсная переходная характеристика имеет бесконечную длину во временной области, а передаточная функция имеет дробно-рациональный вид. Такие фильтры могут быть как аналоговыми, так и цифровыми.
Примерами БИХ-фильтров являются фильтр Чебышёва, фильтр Баттерворта, Фильтр Калмана и фильтр Бесселя.
В этих примерах предположим, что вы ежедневно в течение года снимаете показания какого-то процесса. Предположим, что процесс имеет колебания примерно в недельном и месячном масштабах. Кроме того, в данных наблюдается низкочастотный восходящий тренд и аддитивный 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))
Порядок КИХ-фильтра равен 78
fprintf("The order of the IIR filter is %d\n",filtord(Hd2))
Порядок действия фильтра БИХ равен 8
[phifir,wf] = phasez(Hd1,[],1);
[phiiir,wi] = phasez(Hd2,[],1);
figure
plot(wf,unwrap(phifir))
hold on
plot(wi,unwrap(phiiir))
hold off
xlabel("Cycles/Day")
ylabel("Radians")
legend("FIR Equiripple Filter","IIR Butterworth Filter")
grid
БИХ-фильтр имеет гораздо меньший порядок, чем КИХ-фильтр. Однако КИХ-фильтр имеет линейную фазовую характеристику в полосе пропускания, а БИХ-фильтр — нет. КИХ-фильтр одинаково задерживает все частоты в полосе пропускания, а БИХ-фильтр — нет.
Кроме того, скорость изменения фазы на единицу частоты в КИХ-фильтре выше, чем в БИХ-фильтре.
Для сравнения спроектируйте КИХ-фильтр нижних частот с одинаковой пульсацией. Характеристики фильтра нижних частот: полоса пропускания [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 на выходе полосно-пропускающего КИХ-фильтра.
[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")
Увеличение фазовой задержки в КИХ-фильтре заметно на выходе фильтра.
Для сравнения постройте график выходного сигнала КИХ-фильтра нижних частот, наложенного на 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 длины КИХ-фильтра в выборках (сэмплах) (10×0.001).
Эта задержка обусловлена фазовой характеристикой фильтра. В этих примерах КИХ-фильтр является линейно-фазовым фильтром первого типа. Групповая задержка фильтра составляет 10 отсчётов.
Постройте график групповой задержки.
[gd,f] = grpdelay(Hd,[],Fs);
plot(f,gd)
xlabel("Frequency (Hz)")
ylabel("Group Delay (samples)")
grid
Во многих случаях фазовое искажение допустимо. Это особенно верно, когда фазовая характеристика является линейной. В других случаях желательно иметь фильтр с нулевой фазовой характеристикой. Нулевая фазовая характеристика технически невозможна в некаузальном фильтре. Однако нулевую фазовую фильтрацию можно реализовать с помощью каузального фильтра с 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 нет задержки из-за фазовой характеристики КИХ-фильтра.