Загрузите MAT-файл chirp. Файл содержит сигнал, y, дискретизированный с частотой Fs = 8192 Гц. Сигнал имеет большую часть степени выше Fs/ 4 = 2048 Гц, или половина частоты Найквиста. Добавьте случайный шум к сигналу.

load chirp
y = y + randn(size(y))/25;
t = (0:length(y)-1)/Fs;

Спроектируйте фильтр верхней конечной импульсной характеристики 34-го порядка, чтобы ослабить компоненты сигнала ниже Fs/4. Задайте нормированную частоту среза 0,48, которая соответствует приблизительно 1966 Гц. Визуализируйте частотную характеристику фильтра.

f = [0 0.48 0.48 1];
mhi = [0 0 1 1];
bhi = fir2(34,f,mhi);

freqz(bhi,1,[],Fs)

Пропустите сигнал щебета. Постройте график сигнала до и после фильтрации.

outhi = filter(bhi,1,y);

figure
subplot(2,1,1)
plot(t,y)
title('Original Signal')
ylim([-1.2 1.2])

subplot(2,1,2)
plot(t,outhi)
title('Higpass Filtered Signal')
xlabel('Time (s)')
ylim([-1.2 1.2])

Смените фильтр с верхнего на lowpass. Используйте тот же порядок и отключение. Снова отфильтровайте сигнал. Результатом в основном является шум.

mlo = [1 1 0 0];
blo = fir2(34,f,mlo);
outlo = filter(blo,1,y);

subplot(2,1,1)
plot(t,y)
title('Original Signal')
ylim([-1.2 1.2])

subplot(2,1,2)
plot(t,outlo)
title('Lowpass Filtered Signal')
xlabel('Time (s)')
ylim([-1.2 1.2])

Спроектируйте lowpass 30-го порядка с нормализованной частотой среза $$0.6\pi$$ рад/образец. Постройте график идеальной частотной характеристики, наложенной на фактическую частотную характеристику.

f = [0 0.6 0.6 1];
m = [1 1 0 0];

b1 = fir2(30,f,m);
[h1,w] = freqz(b1,1);

plot(f,m,w/pi,abs(h1))
xlabel('\omega / \pi')
lgs = {'Ideal','fir2 default'};
legend(lgs)

Перепроектируйте фильтр с помощью сетки интерполяции с 64 точками.

b2 = fir2(30,f,m,64);
h2 = freqz(b2,1);

hold on
plot(w/pi,abs(h2))
lgs{3} = 'npt = 64';
legend(lgs)

Перепроектируйте фильтр с помощью сетки интерполяции с 64 точками и интервала с 13 точками вокруг частоты среза.

b3 = fir2(30,f,m,64,13);
h3 = freqz(b3,1);

plot(w/pi,abs(h3))
lgs{4} = 'lap = 13';
legend(lgs)

Спроектируйте конечная импульсная характеристика со следующей частотной характеристикой:

F1 = 0:0.01:0.18;
A1 = 0.5+sin(2*pi*7.5*F1)/4;

F2 = [0.2 0.38 0.4 0.55 0.562 0.585 0.6 0.78];
A2 = [0.5 2.3 1 1 -0.2 -0.2 1 1];

F3 = 0.79:0.01:1;
A3 = 0.2+18*(1-F3).^2;

Проектируйте фильтр с помощью окна Хэмминга. Задайте порядок фильтра 50.

N = 50;

FreqVect = [F1 F2 F3];
AmplVect = [A1 A2 A3];

ham = fir2(N,FreqVect,AmplVect);

Повторите вычисление с помощью окна Кайзера, которое имеет параметр формы 3.

kai = fir2(N,FreqVect,AmplVect,kaiser(N+1,3));

Перепроектируйте фильтр с помощью designfilt функция. designfilt по умолчанию используется прямоугольное окно. Вычислите нулевую фазовую характеристику фильтра по 1024 точкам.

d = designfilt('arbmagfir','FilterOrder',N, ...
    'Frequencies',FreqVect,'Amplitudes',AmplVect);

[zd,wd] = zerophase(d,1024);

Отображение нулевых характеристик трех фильтров. Наложите идеальный ответ.

zerophase(ham,1)
hold on
zerophase(kai,1)
plot(wd/pi,zd)
plot(FreqVect,AmplVect,'k:')
legend('Hamming','Kaiser','designfilt','ideal')