采用窗函数法设计理想低通,高通滤波器,参考北京交通大学陈后金主编的【数字信号处理】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版