滤波器分析还是Matlab做更方便,写的特别简单的fir滤波器设计代码。
其实也也就是套用了函数的框架,注释了一点关键内容方便学习而已。
其实网上很多教程,善用搜索吧。
% fir filter
% Better used in GNSS data analysis.
% wp: 通带边界归一化频率
% ws: 阻带边界归一化频率
% x: data list
% Sp:sample point
% flg: 0:FIR; 1:zero phase shift fir
function y = fir_filter(wp, ws, x, Sp, Fs, flg)
if nargin <6
flg = 0;
end
if nargin<5
Fs = 1;
end
if nargin<4
Sp = 512;
end
if nargin<3 || nargin>6
error('error in fir_filter.');
end
wdelta = ws - wp;
% 过渡带宽等于汉宁窗函数的主瓣宽度。
% 求得滤波器所用常函数的最小带宽
% 汉宁窗: 8db; 三角窗: 4db; 海明窗: 8db; 矩形窗: 4db
N = ceil(8*pi/wdelta);
% 截止频率取为通带和阻带的中点
Wn = 0.5*(wp+ws);
b = fir1(N, Wn/pi, hanning(N+1));
% 512:频率采样点, 50采样点
[H,f] = freqz(b,1,Sp,Fs); % 含义分别为分子,分母(FIR为1),插值频率,采样频率
%
figure(1)
subplot(2,1,1);
plot(f,20*log10(abs(H)));
xlabel('Frequency');
ylabel('Amplify');
grid on;
subplot(212)
plot(f,180/pi*unwrap(angle(H)));
xlabel('Frequency');
ylabel('Amplify');
grid on;
% Normal FIR, Phase Delay.
if flg == 1
y = fftfilt(b,x);
end
% Zero-phase-shift Filter
if flg == 0
y = filtfilt(b,1,x);
end
t = 0:1/Fs:(length(x)-1)/Fs;
figure(2)
subplot(211)
plot(t,x)
title('Before FIR');
xlabel('Time/s');
ylabel('Deviation/m');
xlim([0,max(t)]);
subplot(212)
plot(t,y)
title('After FIR');
xlabel('Time/s');
ylabel('Deviation/m');
xlim([0,max(t)]);
end