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版