采用窗函数法设计理想低通,高通滤波器,及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版

版权声明:本文为xhslovecx原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://www.cnblogs.com/xhslovecx/p/10196851.html