FIR滤波器窗函数设计法详细步骤以及Matlab代码
采用窗函数法设计理想低通，高通滤波器，参考北京交通大学陈后金主编的【数字信号处理】5.2节 窗函数法设计线性相位FIR数字滤波器P164，和P188。
设计步骤如下：
1） 确定滤波器类型，不同的FIR类型可设计不同类型的滤波器，I型可设计LP(低通滤波器)，HP(高通滤波器)，BP(带通滤波器)，BS(带阻滤波器)。
|
Fir I型 |
Fir II型 |
Fir III型 |
Fir IV型 |
|
LP，HP，BP，BS |
LP，BP |
BP |
HP，BP，BS |
2） 确定设计的滤波器的参数
Eg：若要设计一个低通滤波器，fp=20，fs=30；Ap=1，As=40，则3db截频Wc = 2*pi*（fs-fp）/Fs；Fs为采样频率。
当选定某一窗函数时，衰耗Ap和As就已经确定，凯撒窗除外。Ap和As的计算方法可参看另外一篇博客： https://www.cnblogs.com/xhslovecx/p/10118570.html
3） 确定窗函数
|
窗的类型 |
主瓣宽度 |
近似过渡带宽度 |
δp，δs |
Ap（dB） |
As（dB） |
|
矩形窗 |
4pi/N |
1.8pi/N |
0.09 |
0.82 |
21 |
|
Hann |
8pi/N |
6.2pi/N |
0.0064 |
0.056 |
44 |
|
Hamming |
8pi/N |
7pi/N |
0.0022 |
0.019 |
53 |
|
Blackman |
12pi/N |
11.4pi/N |
0.0002 |
0.0017 |
74 |
|
Kaiser |
可调窗，需要确定 β值 |
||||
50<A , β = 0.1102（A-8.7）；
21<=A<=50, β=0.5842(A-21)^0.4 + 0.07886(A-21)；
A<21, β = 0；
4) 确定滤波器的阶数M，首先确定滤波器的长度N。对于除凯撒窗以外的窗函数，N值由以下公式确定：
N>=(窗函数近似过渡带宽度)/（Wp-Ws）
|
Fir I型 |
Fir II型 |
Fir III型 |
Fir IV型 |
|
脉冲响应h[k]为偶对称 |
h[k]为偶对称 |
h[k]为奇对称 |
h[k]为奇对称 |
|
窗函数长度：N=mod(N+1,2)+N |
N=mod(N,2)+N |
N=mod(N+1,2)+N |
N=mod(N,2)+N |
|
阶数M=N-1为偶数 |
M=N-1为奇数 |
M=N-1为偶数 |
M=N-1为奇数 |
若采用Kaiser窗，则
M≈ （A-7.95）÷ 2.285*|Wp-Ws|，A>21。其中，A = -20lg（min（δp，δs））
5) 理想滤波器脉冲信号如下：
hd = (Wc/pi)*sinc(Wc*(k-0.5*M)/pi);%低通
hd = -(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);%高通
6) 加窗：
W = hanning(N); W = hamming(N); W = blackman(N); N = M+1;
W = kaiser(N,beta);
7) 截断
h = hd.*W;
8）滤波
sigFiltered = filter(h,[1],signal);
clear
close all
global Fs
Fs = 360;
load '118m.mat'%mit数据库第118条数据
signal = val(1,100000:111600)/200;
%% 采用FIR I型设计20Hz以下的低通滤波器
fp=14; fs=18;
detap = 0.01;detas = 0.01;
[M,beta] = selectFirFilterN(fp,fs,detap,detas);
N = M+1;
w = kaiser(N,beta);
hd = FIRItypeIdealpulse(fp,fs,N,'low');
h = hd.*w';
% 设计的滤波器
omega = linspace(0,pi,512);
mag = freqz(h,[1],omega);
figure
plot(omega/(2*pi)*Fs,20*log10(abs(mag)));
title('FIR低通(14hz以下)滤波器频率相应');
xlabel('频率');
ylabel('增益(dB)');
%% 采用FIR I型设计8Hz以上的高通滤波器
fp2 = 8; fs2=4;
detap2 = 0.01; detas2 = 0.01;
[M2,beta2] = selectFirFilterN(fp2,fs2,detap2,detas2);
N2 = M2+1;
w2 = kaiser(N2,beta2);
hd2 = FIRItypeIdealpulse(fp2,fs2,N2,'high');
h2 = hd2.*w2';
% 设计的滤波器
omega = linspace(0,pi,512);
mag = freqz(h2,[1],omega);
figure
plot(omega/(2*pi)*Fs,20*log10(abs(mag)));
title('FIR高通(8Hz以上)滤波器频率相应');
xlabel('频率');
ylabel('增益(dB)');
%% 信号滤波
sigFiltered = filter(h,[1],signal);
sigFiltered2 = filter(h2,[1],sigFiltered);
figure
subplot(3,1,1);
plot(signal,'r');
subplot(3,1,2);
hold on
plot(sigFiltered(M/2:end),'b');
subplot(3,1,3);
hold on
plot(sigFiltered2(M/2+M2/2:end),'g');
ylim([-2,2]);
title('信号滤波');
%% 傅里叶变换画出滤波后的频谱
data = FilteredSignal;
M = length(data);
N = M*2-1;
X = fft(data,N);
f = [0:M-1]*Fs/N;
figure
Xabs = abs(fftshift(X));
plot(f(1:end/2),Xabs(M:end-M/2));
title('滤波后的信号频谱');
function [M,beta] = selectFirFilterN(fp,fs,detap,detas)
% 自动选择kaiser窗对应的M和beta值
global Fs
wp = 2*pi*(fp/Fs);
ws = 2*pi*(fs/Fs);
A = -20*log10(min(detap,detas));
M = ceil((A-7.95)/(2.285*abs(wp-ws)));
M = mod(M,2)+M;
% 确定beta的值
if A<21
beta = 0;
elseif A>=21 && A<=50
beta = 0.5842*(A-21)^0.4+0.07886*(A-21);
else
beta = 0.1102*(A-8.7);
end
function hd = FIRItypeIdealpulse(fp,fs,N,type)
%==================================================
% 理想FIR I型低通滤波器,wc是截止角频率，阶数M
%==================================================
global Fs
wp = 2*pi*(fp/Fs);
ws = 2*pi*(fs/Fs);
wc = (wp+ws)/2;
N = mod(N+1,2)+N;
M = N-1;
k = 0:M;
if strcmp(type,'high')
hd = -(wc/pi)*sinc(wc*(k-0.5*M)/pi);
hd(0.5*M+1) = hd(0.5*M+1)+1;
% hd = sinc(k-0.5*M)-(wc/pi)*sinc(wc*(k-0.5*M)/pi);
elseif strcmp(type,'low')
hd = (wc/pi)*sinc(wc*(k-0.5*M)/pi);
else
disp('error');
hd = [];
end
end
参考：陈后金主编 数字信号处理 第2版
