数字信号处理（二）窗函数法设计FIR数字滤波器

最新推荐文章于 2024-07-30 16:22:24 发布

格桑蓝莲

最新推荐文章于 2024-07-30 16:22:24 发布

阅读量5.8w

点赞数 84

分类专栏：工程类项目文章标签：数字信号处理窗函数法设计滤波器各种滤波器 matlab仿真信号分析处理

本文链接：https://blog.csdn.net/weixin_44566643/article/details/103109690

版权

工程类项目专栏收录该内容

5 篇文章 11 订阅

订阅专栏

数字信号处理（二）窗函数法设计FIR数字滤波器

（一)实验目的

(1) 熟悉矩形窗、汉宁窗、海明窗等常用窗函数。
(2) 掌握用上述窗函数法设计FIR数字滤波器的原理和方法。
(3) 熟悉线性相位FIR数字滤波器特性。
(4) 了解窗函数类型及窗口长度对滤波器特性的影响。

(二)实验原理

无限长冲激响应数字滤波器的优点是可以利用模拟滤波器的设计结果，而模拟滤波器的设计可以查阅大量图表，所以设计方法较为简单，但它有个缺点就是相位非线性。如果需要实现相位线性，则要采用全通网络进行相位校正。也就是说，IIR滤波器很难设计成具有线性相位的。FIR滤波器在保证幅度特性的同时，很容易实现严格的线性相位特性。
设计FIR数字滤波器最常用的方法是窗函数法，又称窗口法。这种方法一般是先给定所要求的理想滤波器的频率响应H_d (jω)，要求设计一个FIR滤波器频率响应H_d (jω)，去逼近理想的频率响应H_d (jω)。然而，窗函数法设计FIR数字滤波器是在时域进行的，因此，必须首先由理想频率响应H_d (jω)的傅里叶反变换推导出对应的单位冲激响应h(n)。
在这里插入图片描述

（三）实验内容与结果分析

1.矩形窗
已知ω_c=π/4,利用矩形窗函数设计低通滤波器。在N=15，N=33两种窗口长度时，所绘制的幅频特性和相频特性如图1-1和1-2所示。
在这里插入图片描述

图1-2 N=33时矩形窗函数滤波器
2.三角窗
已知ω_c=π/4,利用三角窗函数设计低通滤波器。在N=15，N=33两种窗口长度时，所绘制的幅频特性和相频特性如图2-1和2-2所示。

3.汉宁窗
已知ω_c=π/4,利用汉宁窗函数设计低通滤波器。在N=15，N=33两种窗口长度时，所绘制的幅频特性和相频特性如图3-1和3-2所示。
在这里插入图片描述

（四）思考题

（1）如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器? 写出设计步骤。
答：①通过通带截止频率和阻带截止频率确定所设计的滤波器的过渡带宽。
在这里插入图片描述
②通过所需的阻带最小衰减查表确定所选的窗函数类型。
③通过所选的窗函数的过渡带宽公式，确定N值。
④将h_d (n)与窗函数相乘得FIR数字滤波器的单位冲激响应h(n)。

（2）如果要求用窗函数法设计带通滤波器，且给定下、上边带截止频率为ω_1和ω_2，试求理想带通滤波器的单位冲激响应H_d (n)。
答：理想带通滤波器的频率响应为
在这里插入图片描述
利用傅立叶逆变换推导出对应的单位冲激响应

思考题：设计一个线性相位带通滤波器，下阻带截止频率为2Khz,上阻带截止频率为6khz,通带下截止频率为3khz，通带上截止频率为5khz,通带最大衰减为0.2db，阻带最小衰减为55db,抽样频率为20khz
在这里插入图片描述

（五）个人感想与总结

1.通过“窗函数法设计FIR数字滤波器”的实验，我对于滤波器的设计又了更加深入的理解。特别是对于滤波器类型的选择，还有根据所给出的条件进行滤波器的确定还有滤波器窗口长度的确定。
2.对于不同种类的窗口，如：矩形窗、三角窗、汉宁窗、海明窗等都更加熟悉，也知道了它们的表达式和性能特点，以及它们的函数图形。
3.我也通过实验，知道了已知单位冲激响应，要在Matlab中求其幅频特性，并进行作图的方法。对于FIR滤波器的知识也有了更扎实的掌握。对我理解窗函数设计滤波器方法很有帮助，我相信，这次数字信号处理实验的参与仅仅是打开了我利用matlab信号分析与处理的大门，在今后我一定会不断努力加油让matlab真正成为我的一项技能。

（六）附录：程序代码

Ideallp函数：
%ideallp函数 理想线性相位低通滤波器计算
%wc 截止频率 N--理想线性滤波器的长度
function hd=ideallp(wc,N)
tao=(N-1)/2;
n=[0:1:(N-1)];
m=n-tao+eps;
hd=sin(wc*m)./(pi*m);
	矩形窗函数：
%本实验模块用于产生矩形窗
clc;clear;
N=[15,33];
for i=1:2
    figure(i);
    wc=pi/4;%理想低通滤波器的截止频率
    n=[0:N(i)-1];
    wr=rectwin(N(i))';%得到一个矩形窗向量wr，分别对应15个1和33个1
    hd=ideallp(wc,N(i));%理想线性相位低通滤波器的单位冲激响应
    h=hd.*wr;%实际低通滤波器的单位冲激响应
    [Hw,w]=freqz(h,[1],1024,'whole');
    db=20*log10(abs(Hw));
    subplot(2,2,1);
    plot(n,wr);grid;
    title('矩形窗');
    ylabel('w(n)');%画出矩形窗
    subplot(2,2,2);
    stem(n,h);
    title('h(n)波形');%画出实际滤波器单位冲激响应
    subplot(2,2,3)
    plot(w/pi,db);
    title('幅频特性（dB）');
    xlabel('\omega/\pi');
    ylabel('20log|H(e^j^\omega)|(dB)');
    set(gca,'XLim',[0 1]);%将横坐标限制在0到1之间画出幅频特性
    subplot(2,2,4);
    theta=unwrap(angle(Hw));
    plot(w/pi,theta);
    title('相频特性');
    xlabel('\omega/\pi');
    set(gca,'XLim',[0 1]);%将横坐标限制在0到1之间画出相频特性
end
	三角窗：
%本模块用于产生三角窗
N=[15,33];
for i=1:2;
    figure(i);
    wc=pi/4;%理想低通滤波器的截止频率
    n=[0:N(i)-1];
    wr=triang(N(i))';%得到一个三角窗向量wr，分别对应15个1和33个1
    hd=ideallp(wc,N(i));%理想线性相位低通滤波器的单位冲激响应
    h=hd.*wr; %实际低通滤波器的单位冲激响应
    [Hw,w]=freqz(h,[1],1000,'whole');
    db=20*log10(abs(Hw));
    subplot(2,2,1);
    plot(n,wr);hold on;
    scatter(n,wr);
    set(gca,'XLim',[0 N(i)-1]);%限制横坐标画出三角窗
    title('三角窗');
    xlabel('n');
    ylabel('w(n)');
    subplot(2,2,2);
    stem(n,h);
    set(gca,'XLim',[0 N(i)-1]);%画出实际滤波器的单位冲激响应
    title('h(n)波形');
    subplot(2,2,3)
    plot(w/pi,db);
    title('幅频特性');
    xlabel('\omega/\pi');
    ylabel('20log|H(e^j^\omega)|(dB)');%画出实际滤波器的幅频特性
    set(gca,'XLim',[0 1]);
    subplot(2,2,4);
    theta=unwrap(angle(Hw));
    plot(w/pi,theta);
    title('相频特性');
    xlabel('\omega/\pi');
    set(gca,'XLim',[0 1]);%画出实际滤波器的相频特性
end
	汉宁窗：
%本实验模块用于产生汉宁窗
N=[15,33];
for i=1:2
    figure(i);
    wc=pi/4;%理想低通滤波器的截止频率
    n=[0:N(i)-1];
    wr=hann(N(i))';%得到一个汉宁窗向量wr，分别对应15个1和33个1
    hd=ideallp(wc,N(i));%理想线性相位低通滤波器的单位冲激响应
    h=hd.*wr;%实际低通滤波器的单位冲激响应
    [Hw,w]=freqz(h,[1],1000,'whole');
    db=20*log10(abs(Hw));
    subplot(2,2,1);
    plot(n,wr);grid;
    set(gca,'XLim',[0,N(i)-1]);
    title('汉宁窗');
    xlabel('n');
    ylabel('w(n)');%画出汉宁窗向量
    subplot(2,2,2);
    stem(n,h);
    title('h(n)波形');%画出实际滤波器单位冲激响应
    subplot(2,2,3)
    plot(w/pi,db);
    title('幅频特性（dB）');
    xlabel('\omega/\pi');
    ylabel('20log|H(e^j^\omega)|(dB)');
    set(gca,'XLim',[0 1]);%将横坐标限制在0到1之间画出幅频特性
    subplot(2,2,4);
    theta=unwrap(angle(Hw));
    plot(w/pi,theta);
    title('相频特性');
    xlabel('\omega/\pi');
    set(gca,'XLim',[0 1]);%将横坐标限制在0到1之间画出相频特性
end
（5） 思考题程序：
      %思考题：设计一个线性相位带通滤波器，下阻带截止频率为2Khz,上阻带截止频率为6khz,
%通带下截止频率为3khz，通带上截止频率为5khz,通带最大衰减为0.2db，阻带最小衰减为55db,抽样频率为20k
%将模拟频率转为数字频率，分别为
ws1=0.2*pi;
wp1=0.3*pi;
wp2=0.5*pi;
ws2=0.6*pi;
AS=55;
RP=0.2;
AS0=50;
RP0=1;
deltaw=min((wp1-ws1),(ws2-wp2));%计算过渡带宽
N=ceil((AS-7.95)/(2.286*deltaw));%凯泽窗的长度
beta=0.1102*(AS-8.7);%计算凯泽的β
while(AS0<AS)|(RP0>RP)
    N=N+1;
    n=[0:N-1];
    wd=(kaiser(N,beta))';%求凯泽窗函数；
    wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;%理想带通的截止频率
    hd=ideallp(wc2,N)-ideallp(wc1,N);%理想带通滤波器的单位冲激响应
    h=hd.*wd;%实际带通滤波器的单位冲激响应；
    [H,w]=freqz(h,[1],1000,'whole');
    mag=abs(H);
    db=20*log10((mag+eps)/max(mag));
    dw=2*pi/1000;
    RP0=-min(db(wp1/dw+1:1:wp2/dw+1));%检查通带最大衰减
    AS0=-max(db(ws2/dw+1:1:501));
    if N>300
        error('没有找到合适的窗函数');
    end
end
subplot(121);stem(n,h,'.');grid;
title('凯泽窗');xlabel('n');ylabel('w(n)');
axis([0,N,-0.3,0.3]);
subplot(122);
plot(w/pi,db);
title('幅度响应');xlabel('\omega/\pi');
ylabel('20log|H(e^j^\omega)|(db)');
axis([0,1,-80,5]);grid;
[N,RP0,AS,beta];