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版