数字信号处理研讨

《数字信号处理》课程研究性学习报告

数字滤波器设计专题研讨

【目的】

(1) 掌握IIR和FIR数字滤波器的设计方法及各自的特点。 (2) 掌握各种窗函数的时频特性及对滤波器设计的影响。

(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。

【研讨题目】 基本题

1．分析矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的频域特性，并进行比较。

【题目分析】分析不同的窗函数的频率特性，可以看出其主瓣宽度与旁瓣宽度的差异，过渡带的宽度的差异。

【仿真结果】

【结果分析】

各种窗特点：各种窗函数都采用了相同的长度，计算fft的长度均为512个点。从结果可以看出矩形窗的主瓣幅度最大，但其的宽度最小，并且其旁瓣幅度也较大。其他的几种窗函数恰恰与矩形窗相反，blackman窗与kaiser窗的旁瓣幅度很小，几乎为0。

【自主学习内容】

1、各种窗函数的调用；

2、利用fftshift函数处理FFT的结果，便于观察。

【阅读文献】

数字信号处理教材及MATLAB教程

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

探讨Kaiser窗时bate值的改变可以改变窗函数的形状，从而达到不同的阻带衰减。bate值为0时，Kaiser窗就是矩形窗，并且随着bate值增加，Kaiser窗在两端的衰减逐渐加大。

【问题探究】

在谱分析中如何选择窗函数，在滤波器设计中如何选择窗函数？

答：根据设计要求的Ωp、Ωs和Ap、As确定滤波器的Ωc和窗函数的类型及其长度N,再确定窗函数的幅度函数和相位函数，计算窗函数的IDTFT得到hd[k],加窗截断而得到h[k]。

【仿真程序】

L=512;N=20; figure(1) w1=zeros(1,100); w2=ones(1,N); w3=zeros(1,100); wh1=[w1 w2 w3];

WH1=fftshift(fft(wh1,L)); w=(0:L-1)-L/2; plot(w,abs(WH1)); title('矩形窗的频域图'); figure(2) wh2=hann(N)';

WH2=fftshift(fft(wh2,L)); w=(0:L-1)-L/2; plot(w,abs(WH2)); title('汉纳窗的频域图'); figure(3) wh3=hamming(N)';

WH3=fftshift(fft(wh3,L)); w=(0:L-1)-L/2; plot(w,abs(WH3));

title('汉明窗的频域图'); figure(4)

wh4=blackman(N)'; WH4=fftshift(fft(wh4,L)); w=(0:L-1)-L/2; plot(w,abs(WH4));

title('blackman的频域图'); figure(5) bate=20;

wh5=kaiser(N,bate); WH5=fftshift(fft(wh5,L)); w=(0:L-1)-L/2; plot(w,abs(WH5)); title('kaiser窗的频域图');

【研讨题目】 基本题

2．（M5-5）在用窗口法设计FIR滤波器时，由于理想滤波器的频幅响应在截频处发生突变，使得设计出的滤波器的频幅响应发生振荡，这个现象被称为Gibbs现象。解决这个问题的一个方案是本书中介绍的用逐步衰减的窗函数。另一个方案是使理想滤波器过渡带为渐变的，如下图所示具有线性过渡带的理想低通滤波器的频率响应，试用窗口法设计逼近该频率响应的FIR滤波器。

HL(ej)spps

题2图

【设计步骤】

（1）根据所需设计的滤波器，确定线性相位滤波器的类型。

（2）确定理想滤波器的幅度函数Ad（Ω）。

（3）确定理想滤波器的相位Φd（Ω）=-0.5MΩ+β。 （4）计算hd[k]的IDTFT。

（5）用窗函数截断hd[k]而得到h[k]。

【单位脉冲响应证明】

试证该滤波器的单位脉冲响应为

c k0,π, hL[k]

2sin(Δk/2)sin(ck) k0Δkπk其中：Δsp，c(ps)/2

【仿真结果】 M1=8

M1=16

M1=64

M1=128

【结果分析】

FIR滤波器增加采样点，即提高分辨率，对于阻带波动没有效果。而滤波器渐变设计提高分辨率能够减少波动。

【自主学习内容】 【阅读文献】 数字信号处理教材 【发现问题】 【问题探究】 【仿真程序】

wp=0.7*pi; ws=0.3*pi; Ap=1; As=30;

N=ceil(7*pi/(wp-ws)); N=mod(N+1,2)+N M=N-1; w=hamming(N);

wc=(wp+ws)/2; k=0:M;

hd=(wc/pi)*sinc(wc*(k-0.5*M)/pi); h=hd'.*w;

omega=linspace(0,pi,512); mag=freqz(h,[1],omega); magdb=abs(mag);

plot(omega/pi,magdb,'b'); grid; w=ws-wp; M1=1; k2=-M1:M1; wc=(wp+ws)/2;

hd=sinc(w*k2/2).*(sin(wc*k2)./(k2.*pi)); hd(M1+1)=wc/pi;

omega2=linspace(0,pi,512); mag2=freqz(hd,[1],omega2); magdb2=abs(mag2); hold on;

plot(omega2/pi,magdb2,'r'); legend('逐步','渐变'); grid on;

【研讨题目】中等题

3．Dhexian.wav是对频率为293.66, 369.99, 440Hz的D大调和弦以8000Hz抽样所得的数字音乐信号，试设计一数字滤波器从和弦中分离出440Hz的音符。要求：

（1）设计IIR数字高通滤波器，通过实验研究P,s，AP,As的选择对滤波效果及滤波器阶数的影响，给出滤波器指标选择的基本原则，确定你认为最合适的滤波器指标。

（2）能否用IIR数字带通滤波器从和弦中分离出440Hz的音符？利用（1）确定的基本原则，给出数字带通滤波器的指标。设计IIR数字带通滤波器，并将结果与高通滤波器比较，给出你的结论。

（3）用窗函数法设计FIR数字高通滤波器，分别利用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗截断。讨论用窗函数法设计FIR数字高通滤波器时如何确定滤波器的指标，比较相同过渡带时用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗设计滤波器的阶数。

（4）采用Parks-McClellan算法，设计FIR数字高通滤波器。试参照（1）确定的最合适的高通滤波器指标，给出FIR数字高通滤波器的指标。将设计结果与（1）中的IIR数字滤波器，从幅度响应、相位响应、滤波器阶数等方面进行比较。

【温磬提示】

在IIR数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。

【设计步骤】

1、首先应先画出原信号的时域和频域图

2、IIR滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，则通过频率变换转换为设计相应的高通，带通等。在设计的全过程的各个步骤，matlab都提供相应的工具箱函数，使得IIR数字滤波器设计变得非常简单。总的来说，我的设计思路主要有以下两种：

思路一：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。

归一化模拟模拟高,带通 模拟域 冲激响应不变法 数字高,带通

低通原型 或带阻 频率变换 双线性变换法 图2-1 先频率变换再离散

或带阻 思路二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。

数字域 归一化模拟数字原型低低通原型 通双线性变换法 频率变换 图2-2 先离散再频率变换

数字高,带通或带阻

【仿真结果】

矩形窗

hann窗

Hamming窗

Blackman窗

Kasier窗

（4） ds=0.001

ds=0.01

ds=0.1

【结果分析】

1、通过仿真可得到，若想取得较好的实验结果，必须选取合适的fp.fs,ap,as值，才能是滤波效果较好。知可选择Wp=0.85pi,Ws=0.78pi,Ap=1,As=50指标进行滤波。

2、知道通过仿真能够进行滤波，对于带通滤波器，选取合适的参数可以完成与高通滤波器同样的滤波效果

3、在相同截取长度下，衰减逐渐增加，同时相应的近似过渡带也逐渐加宽。可以看出在增加衰减的同时会增加过渡带的宽度，所以在设计滤波器时要合理的选择窗函数。而凯泽窗则是完全与实际的衰减幅度类似。可以自主调节所选窗函数的长度。

4、可知，Parks-McClellan算法设计出的滤波器阻带的衰减是等波纹的。且通过改变波动的幅度可以控制衰减的大小。

相位响应比较：通过相位响应我们可以看出BW行IIR滤波器的线性性比较好，而其余的IIR滤波器的相位响应都有非线性失真，即相位非线性。而FIR滤波器在通带区间的相位响应都是直线，即线性的。

幅度响应比较：通过IIR滤波器和FIR滤波器的幅度响应可以看出，不同方法设计的IIR滤波器的通阻带衰减波动不一样，而FIR滤波器都有较大的波动。两者在衰减幅度方面都能达到滤波器设计的要求。都比较符合设计。

阶数方面：通过各个图形的N值可以看出，设计相同效果的滤波器，FIR滤波器的阶数

更少。

【自主学习内容】

【阅读文献】

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

【问题探究】

【仿真程序】

（1）Wp=0.85*pi; Ws=0.78*pi; Ap=1; As=50;

T=2;Fs=1/T; wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T;wp1=1/wp;ws1=1/ws; [N,wc]=buttord(wp1,ws1,Ap,As,'s');[num,den]=butter(N,wc,'s');[numa,dena]=lp2hp(num,den,1);[numd,dend]=bilinear(numa,dena,Fs); w=linspace(0,pi,1024);h=freqz(numd,dend,w);

plot(w/pi,20*log10(abs(h)));axis([0 1 -50 0]);grid;xlabel('Normalized frequency');ylabel('Gain,dB');

（2）M=64;Wp1=0.8*pi;Wp2=0.9*pi; Wp=(Wp1+Wp2)/2;m=0:M/2;

Wm=2*pi*m./(M+1);

mtr1=ceil(Wp1*(M+1)/(2*pi)); mtr2=floor(Wp2*(M+1)/(2*pi))+2; Ad=double([(Wm>=Wp1)&(Wm<=Wp2)]); Ad(mtr1)=0.29; Ad(mtr2)=0.30;

Hd=Ad.*exp(-j*0.5*M*Wm);

Hd=[Hd conj(fliplr(Hd(2:M/2+1)))]; h=real(ifft(Hd)); w=linspace(0.1,pi,1000); H=freqz(h,[1],w);

plot(w/pi,20*log10(abs(H)));grid;

（3） 矩形窗

clear all Fsamp = 10e3; fp = 4.2e3; fs = 3.8e3; Ap = 1; As = 50;

Wp = 2*pi*fp/Fsamp; %¹éÒ»»¯

Ws = 2*pi*fs/Fsamp; N=ceil(1.8*pi/(Wp-Ws)); N=mod(N+1,2)+N; M=N-1;

w=ones(1,N);fprintf('N=%.0f\\n',N); Wc=(Wp+Ws)/2; k=0:M;

hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1; h=hd.*w;

omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=20*log10(abs(mag));

plot(omega/pi,magdb); grid;xlabel('Normalized frequency');ylabel('Gain,dB');

N=23

hann窗

wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;

N=ceil(6.4*pi/(wp-ws));N=mod(N+1,2)+N;M=N-1;w=hann(N);wc=(wp+ws)/2;k=0:M;

hd=-(wc/pi)*sinc(wc*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1; h=hd'.*w;

omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=20*log10(abs(mag));

plot(omega/pi,magdb,'b');grid;w=ws-wp;

N=79

Hamming窗

wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;

N=ceil(6.4*pi/(wp-ws));N=mod(N+1,2)+N;M=N-1;w=hamming(N);wc=(wp+ws)/2;k=0:M;

hd=-(wc/pi)*sinc(wc*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1; h=hd'.*w;

omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=20*log10(abs(mag));

plot(omega/pi,magdb,'b');grid;w=ws-wp;

N=89

Blackman窗 N=143

wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;

N=ceil(6.4*pi/(wp-ws));N=mod(N+1,2)+N;M=N-1;w=balckman(N);wc=(wp+ws)/2;k=0:M;

hd=-(wc/pi)*sinc(wc*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1;

h=hd'.*w;

omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=20*log10(abs(mag));

plot(omega/pi,magdb,'b');grid;w=ws-wp;

Kasier窗 N=74

wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;

N=ceil(6.4*pi/(wp-ws));N=mod(N+1,2)+N;M=N-1;beta=0.1102*(As-8.7);w=kaiser(N,beta);wc=(wp+ws)/2;k=0:M;

hd=-(wc/pi)*sinc(wc*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1; h=hd'.*w;

omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=20*log10(abs(mag)); plot(omega/pi,magdb,'b');grid;w=ws-wp;

（4）

Fp=0.84;Fs=0.76;ds=0.001;dp=ds; f=[Fs Fp];a=[0 1];dev=[ds dp]; [M,fo,ao,w] = remezord(f,a,dev); h = remez(M,fo,ao,w); w=linspace(0,pi,1000); mag=freqz(h,[1],w);

plot(w/pi,20*log10(abs(mag))); xlabel('Normalized frequency'); ylabel('Gain, dB');grid;

相位谱程序

Wp=0.84*pi; Ws=0.76*pi; Ap=1; As=50;

T=2;Fs=1/T; wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T;wp1=1/wp;ws1=1/ws; [N,wc]=buttord(wp1,ws1,Ap,As,'s');[num,den]=butter(N,wc,'s');[numa,dena]=lp2hp(num,den,1);[numd,dend]=bilinear(numa,dena,Fs); w=linspace(0,pi,1024);h=freqz(numd,dend,w); plot(w/pi, angle(h)); grid;xlabel('Normalized

frequency');ylabel('Gain,dB');title(' BWÐÍÏàλÏìÓ¦'); Wp=0.84*pi;Ws=0.76*pi;Ap=1;As=50;

N=ceil(1.8*pi/(Wp-Ws));N=mod(N+1,2)+N;M=N-1;w=ones(1,N);fprintf('N=%.0f\\n',N);

Wc=(Wp+Ws)/2;k=0:M;

hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1; h=hd.*w;

omega=linspace(0,pi,512);mag=freqz(h,[1],omega);plot(omega/pi,angle(mag));

grid;xlabel('Normalized frequency');ylabel('Gain,dB');

【研讨题目】 (附加题)

3．试用频率取样法设计频率响应逼近HDIF(ej)j,用III型IV型线性相位滤波器，并将所得结果进行比较。

【题目分析】

可利用III型滤波器的频率响应函数和频率取样法设计FIR滤波器

【设计步骤】 【仿真结果】

π的FIR数字微分器。分别采

【结果分析】

【自主学习内容】

【阅读文献】

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

【问题探究】

【仿真程序】 III型 M=64;Wp1=-pi;Wp2=pi;m=0:M/2; Wm=2*pi*m./(M+1);

Ad=double([(Wm>=Wp1)&(Wm<=Wp2)]).*Wm; Hd=j*Ad.*exp(-j*0.5*M*Wm); Hd=[Hd conj(fliplr(Hd(2:M/2+1)))]; h=real(ifft(Hd)); w=linspace(-pi,pi,1000); H=freqz(h,[1],w);

plot(w/pi,20*log10(abs(H)));axis([-1 1 -50 20]); grid;

IV型 M=63;Wp1=-pi;Wp2=pi;m=0:M/2; Wm=2*pi*m./(M+1);

Ad=double([(Wm>=Wp1)&(Wm<=Wp2)]).*Wm; Hd=j*Ad.*exp(-j*0.5*M*Wm);

Hd=[Hd conj(fliplr(Hd(2:(M+1)/2)))]; h=real(ifft(Hd)); w=linspace(-pi,pi,1000); H=freqz(h,[1],w);

plot(w/pi,20*log10(abs(H)));grid;