4. Реализация спектрального анализа с использованием бпф
Произведем спектральный анализ данных сигналов. Для того, чтобы избежать явления Гибсона, исходный сигнал умножается на функцию окна. В данном случае будем использовать треугольное окно и окно Хэмминга.
На графиках построены спектры входных сигналов, умноженные на функцию окна.
Окно Чебышева:
% окно Чебышева
window1=chebwin(n,r);
% спектр входного сигнала с использованием окна Чебышева
FT=fft(window1.*y');
FT=FT(1:n/2);
figure(3);
clf, plot(step,abs(FT));
axis([0 25 0 max(abs(FT))]);
ylabel('Magnitude');
xlabel('Frequency (Hertz)');
title('Power spectrum (triangle win)');
На графике видно, что фундаментальная частота нашего сигнала находится в районе 3 Гц. Но сигнал немного зашумлен.
% спектр реального сигнала с использованием окна Чебышева
FT=fft(window1.*yf);
FT=FT(1:n/2);
figure(4);
clf, plot(step,abs(FT));
axis([0 25 0 max(abs(FT))]);
ylabel('Magnitude');
xlabel('Frequency (Hertz)');
title('Power spectrum of real signal (triangle win)');
Здесь трудно определить, что это за сигнал, так как он слишком зашумлен. Но ориентировочно можно сказать, что фундаментальная частота находится в районе 1Гц.
Окно Ханнинга:
% окно Ханнинга
window2=hanning(n);
% спектр входного сигнала с использованием окна Ханнинга
FT=fft(window2.*y');
FT=FT(1:n/2);
figure(5);
clf, plot(step,abs(FT));
axis([0 25 0 max(abs(FT))]);
ylabel('Magnitude');
xlabel('Frequency (Hertz)');
title('Power spectrum (hamming win)');
% спектр реального сигнала с использованием окна Ханнинга
FT=fft(window2.*yf);
FT=FT(1:n/2);
figure(6);
clf, plot(step,abs(FT));
axis([0 25 0 max(abs(FT))]);
ylabel('Magnitude');
xlabel('Frequency (Hertz)');
title('Power spectrum of real signal (hamming win)');
5. Проектирование и реализация цф
Синтез низкочастотного БИХ-фильтра с использованием окна Чебышева 1 рода:
% проектирование и реализация ЦФ
Ws=25; % частота пропускания
Wn = (Ws/Fs)*2; % частота среза
len = 21; % длина фильтра
[b,a] = CHEBY1(len,0.5,Wn);% проектирование НЧ-фильтра Чебышева 1 рода
[H,F] = freqz(b,a,n,Fs);
figure(7);
clf, plot (F,abs(H));
ylabel('Magnitude')
xlabel('Frequency (Hertz) ')
title('Filter Response')
Теперь проведем фильтрацию заданных сигналов.
% фильтрация входного сигнала
format long e
fily=filter(b,1,y);
figure(8);
clf, plot(t,fily);
ylabel('Amplitude (mkV)');
xlabel('Time');
title('Input signal after filtration');
Определим частотный спектр отфильтрованного сигнала:
% спектр отфильтрованного входного сигнала
FT=fft(window1.*fily');
FT=FT(1:n/2);
figure(9);
clf, plot(step,abs(FT));
axis([0 25 0 max(abs(FT))]);
ylabel('Magnitude');
xlabel('Frequency (Hertz)');
title('Power spectrum after filtration');
Аналогично для реального сигнала:
% фильтрация реального сигнала
filyf=filter(b,1,yf');
figure(10);
clf, plot(t,filyf);
ylabel('Amplitude (mkV)');
xlabel('Time');
title('Real signal after filtration');
% спектр отфильтрованного реального сигнала
FT=fft(window1.*filyf');
FT=FT(1:n/2);
figure(11);
clf, plot(step,abs(FT));
axis([0 25 0 max(abs(FT))]);
ylabel('Magnitude');
xlabel('Frequency (Hertz)');
title('Power spectrum of real signal after filtration');
В связи с тем, что сигналы плохо отфильтрованы, то попробуем изменить полосу пропускания фильтра (возьмем её равной 15Гц).
% проектирование и реализация ЦФ
Ws=15; % частота пропускания
Wn = (Ws/Fs)*2; % частота среза
len = 21; % длина фильтра
[b,a] =CHEBY1(len,0.5,Wn);% проектирование НЧ-фильтра Чебышева 1 рода
[H,F] = freqz(b,a,n,Fs);figure(7);
clf, plot (F,abs(H));
ylabel('Magnitude')
xlabel('Frequency (Hertz) ')
title('Filter Response')
Тогда сигналы и их спектры будут иметь следующий вид.
Попробуем теперь взять частоту пропускания равной 5Гц.
% проектирование и реализация ЦФ
Ws=5; % частота пропускания
Wn = (Ws/Fs)*2; % частота среза
len = 21; % длина фильтра
[b,a] =CHEBY1(len,0.5,Wn);% проектирование НЧ-фильтра Чебышева 1 рода
[H,F] = freqz(b,a,n,Fs);
figure(7);
clf, plot (F,abs(H));
ylabel('Magnitude')
xlabel('Frequency (Hertz) ')
title('Filter Response')
Тогда сигналы и их спектры будут иметь следующий вид.
