华北电力大学 数字信号处理课程设计实验报告

华北电力大学实验报告

实验环境实验名称

MATLAB6.5实验一：FFT的应用

1、熟悉MATLAB在数字信号处理中的应用。

2、掌握利用FFT计算序列线性卷积的基本原理及编程实现。

3、掌握对连续信号进行采样的基本原理和方法，并利用FFT对信号进行频谱分析。

1.利用FFT计算线性卷积

1.周期卷积是线性卷积以L为周期的周期延拓序列的主值序列；

2.两个长度为M，N的序列的线性卷积可用长度均为L的循环卷积来代替，其中，L应满足L>M+N-1。而循环卷积可以用FFT来计算，这样我们就可以利用FFT来计算两个序列的线性卷积。

因此，为了使循环卷积与线性卷积的计算结果相等，必须通过补零，将x1(n)和x2（n）有延长到L=N+M-1点。

因此，线性卷积可由下列步骤完成：

（1）将x1(n)和x2(n)都延长到L=M+N-1点；

（2）计算x1(n)的L点FFT，即：X1(k)=FFT[x1(n)]；（3）计算x2(n)的L点FFT，即：X2(k)=FFT[x2(n)];（4）计算Y(k)=X1(k)*X2(k);

（5）计算Y（k）的反变换，即y(n)=IFFT[X1(k)*X2(k)]。

直接计算DFT共需N*N次复数乘法和N(N-1)次复数加法。而FFT仅需计算0.5M次复数乘法和M*N次复数加法。由于在计算机上计算乘法所需的时间比计算加法多得多，所以FFT的运算量比DFT要少的少。2．利用FFT对信号进行谱分析

利用FFT对模拟信号进行谱分析室，应将模拟信号离散化以得到离散时间信号，同时得考虑谱分析中参数的选择。

为避免混叠失真，要求抽样频率fs>2f0（f0连续信号的最高频率），频率分辨率F>离散频率的间隔，记录长度的取样数N（Tp=NT），这三者之间需N>2f0/F。

谱分析的步骤：

首先，利用上面所选参数，在记录长度Tp中对连续时间信号xa(n)进行N点取样，得到离散时间信号x(n)。

然后，利用FFT计算信号的频谱：X(K)=FFT[x(n)]。MATLAB中可以用abs（x）来计算模值。

由于有限长序列补零以后，只是频谱的取样点有所增加，所以不会影响原频谱的分布。

实验目的

实验原理

第1页

1.对于两个序列：x(n)=nR16(n),h(n)=R8(n)

（1）在同一图形窗口中绘出两序列的时域图形。

（2）利用FFT编程计算两序列的线性卷积，绘出其时域图形。

2．利用FFT对信号进行谱分析

对于连续信号xa(t)=cos(2πf1t)+5cos(2πf2t)+cos(2πf3t)，其中f1=6.5kHz,f2=7kHz,f3=9kHz,以采样频率fs=32kHz对其进行采样，

（1）对xa(t)信号采集16点样本，分别作16点和补零到256点的FFT，并分别绘出对应的幅频特性曲线。

（2）对xa(t)信号采集256点样本，分别作256点和512点的FFT，并分别绘出对应的幅频特性曲线。

（3）比较（1）和（2）中的结果，分析采样点数和傅里叶变换点数对FFT的影响，说明高密度频谱和高分辨率频谱的特点与区别。

实验内容

一、FFT计算线性卷积1、程序：n=1:15;

x=n*ones(1,n);h=ones(1,8);subplot(3,1,1)

stem(n,x);title('x(n)=n*R16(n)');subplot(3,1,2)

stem(0:7,h);title('h(n)=R8(n)');x1=[0:15zeros(1,7)];

h1=[ones(1,8)zeros(1,15)];X=fft(x1);H=fft(h1);y=X.*H;Y=ifft(y);subplot(3,1,3)

stem(0:22,Y);title('卷积结果')

ylabel('对数幅度/db');xlabel('以pi为单位的频率')

2、图像：

第2页

实验结果及分析

3、

结果分析：

两个有限长序列的线性卷积可以用循环卷积来代替，而循环卷积可以用FFT来计算。但是，对于x(n)和h（n），都必须把他们延长到N点（N=N1+N2-1），延长的部分用零补充。相对于直接计算方法，利用FFT计算线性卷积可以得到相同的结果，但是FFT算法却可以大大减少运算量，提高运算速度，二、利用FFT对信号进行谱分析1.程序：（1）N=16;L=256;f1=6500;f2=7000;f3=9000;fs=32000;ws=2*pi*fs;T=1/fs;

第3页

n=0:N-1;

x=cos(2*pi*f1*n*T)+5*cos(2*pi*f2*n*T)+cos(2*pi*f3*n*T);y=x(1:N);z=fft(y,N);

w=((0:N-1)*ws/N)/(2*pi);subplot(2,1,1)plot(w,abs(z));

ylabel('幅度特性曲线');xlabel('采样点为16的16点FFT');y2=[x(1:N)zeros(1,L-N)];z2=fft(y2,L);

w=((0:L-1)*ws/L)/(2*pi);subplot(2,2,3)plot(w,abs(z2));

ylabel('幅度特性曲线');xlabl('采样点为16补零到256后的256点FFT');(2)N=256;L=512;f1=6500;f2=7000;f3=9000;fs=32000;ws=2*pi*fs;T=1/fs;n=0:N-1;

x=cos(2*pi*f1*n*T)+5*cos(2*pi*f2*n*T)+cos(2*pi*f3*n*T);y=x(1:N);z=fft(y,N);

w=((0:N-1)*ws/N)/(2*pi);subplot(2,2,1)plot(w,abs(z));

ylabel('幅度特性曲线');xlabel('采样点为256的256点FFT');y2=[x(1:N)zeros(1,L-N)];z2=fft(y2,L);

w=((0:L-1)*ws/L)/(2*pi);subplot(2,2,3)plot(w,abs(z2));

ylabel('幅度特性曲线');xlabl('采样点为256补零到512后的512点FFT');2.图像：

实验结果及分析

第4页

第5

页

实验名称实验目的

设计性实验一：IIR数字滤波器的设计

1、本实验为设计性实验。

2、掌握用双线性变换法设计IIR数字滤波器的基本原理和设计方法。

3、掌握用双线性变换法设计IIR数字Butterworth滤波器的原理和设计方法。IIR数字滤波器的设计借助模拟滤波器原型,再将模拟滤波器转换成数字滤波器。由双线性变换的S域与Z域间的关系可知:z和s之间可以直接代换,由于S平面与Z平面一一单值对应,S平面的虚轴(整个jΩ)经映射后确已成为z平面上的单位圆,但Ω与为非线性关系,因此,通过双线性变换后两个滤波器的频率特性形状不能保持相同,双线性变换不存在混迭效应。因为s平面的左半平面被映射在单位圆内部,这意味着稳定的模拟滤波器经双线性变换可以映射成稳定的数字滤波器。

设计滤波器的步骤：

1.得到数字指标（Wn等）2.双线性变换为模拟低通指标3.归一化模拟拟指标

4.利用通带衰减与阻带衰减的值求Butterworth的阶数N及归一化的模拟系统函数

5.将这个再经过去归一化得到想要的滤波器类型6.用双线性变换法变为数字滤波器

模拟域频率与数字域频率的关系为=(2/T)tan(wc/2),是一种非线性的关系。这种非线性关系使得模拟滤波器和数字滤波器的频率响应与对应频率的关系上发生了畸变，也造成了相位的非线性变化。为保证各边界频率点为预先指定的频率，在确定模拟低通滤波器的系统函数之前，必须按上式进行频率预畸变。

所用子函数：

（1）非归一化巴特沃斯模拟滤波器设计函数：function[b,a]=afd_buttap(wp,ws,Rp,As)

Wp是指通带截止角频率，ws指阻带截止角频率，Rp、As分别指通带最大衰减、阻带最大衰减。

（2）计算离散系统频率响应的函数

function[db,mag,pha,grd,w]=freqz_m(b,a)（3）设计双线性变换法的函数

function[bz,az]=bilinear(b,a,Fs);

实验原理

第6页

实验内容

IIR数字滤波器的设计

用双线性变换法设计一个IIR数字Butterworth低通滤波器。技术指标为：通带截止频率fp=1kHz，阻带截止频率fs=1.5kHz，通带衰减Rp≤1dB，阻带衰减Rs≥40dB，采样频率Fs=10kHz。绘出滤波器的幅频特性曲线和相频特性曲线，判断设计是否符合要求。1.设计方案

设计滤波器，首先要对模拟频率进行数字与转换、归一化，采用双线性变换法还要进行预畸变。接下来依次采用非归一化巴特沃斯模拟滤波器设计函数、设计双线性变换法的函数、计算离散系统频率响应的函数，最后画出幅频特性、相频特性、群延迟等图像即可完成图像。2、源程序：fp=1000;fs=1500;Fs=10000;Rp=1;As=40;T=1/Fs;

wp=2*pi*fp*T;ws=2*pi*fs*T;

omegap=(2/T)*tan(wp/2);omegas=(2/T)*tan(ws/2);

[cs,ds]=afd_buttap(omegap,omegas,Rp,As);[b,a]=bilinear(cs,ds,Fs);

[db,mag,pha,grd,w]=freqz_m(b,a);

subplot(2,2,1);plot(w/pi,mag);ylabel('幅度');

xlabel('以π为单位的频率');title('幅度响应');axis([0,0.801]);subplot(2,2,3);plot(w/pi,db);title('幅度响应（dB）');grid;

xlabel('以π为单位的频率');ylabel('对数幅度dB');axis([0,0.8-600]);

subplot(2,2,2);plot(w/pi,pha);title('相位响应');grid;

xlabel('以π为单位的频率');ylabel('相位');axis([0,0.8-44]);subplot(2,2,4);plot(w/pi,grd);title('群延迟');grid;

xlabel('以π为单位的频率');ylabel('样本');axis([0,0.8010]);子程序：

afd_buttap函数

function[b,a]=afd_buttap(wp,ws,Rp,As)

N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws)));fprintf('\n****butter=%2.0f\n',N)omega=wp/((10^(Rp/10)-1)^(1/(2*N)));[b,a]=buttap_o(N,omega);

freqz_m函数

function[db,mag,pha,grd,w]=freqz_m(b,a)

第7页

[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);

db=20*log((mag+eps)/max(mag));pha=angle(H);

grd=grpdelay(b,a,w);

实验结果分析

巴特沃斯低通滤波器的幅度响应曲线在通带内具有最平坦的特性，且在通带和阻带内幅度特性都是单调变化的，在阻频带最终下降为零。

由于实验要求IIR的通带截止频率是1KHz，阻带截止频率是1.5KHz。经过归一化后得到过渡带的范围就是0.2pi到0.3pi之间。在通带范围内的响应时1，阻带内是0。满足实验设计要求。

第8页

实验名称实验目的

设计性实验二：FIR数字滤波器的设计

1、本实验为设计性实验。

2、掌握用窗函数法设计IIR数字滤波器的基本原理和设计方法。3、掌握用窗函数法设计线性相位FIR数字低通滤波器的编程实现。理想低通滤波器经加窗处理后，其频率响应在不连续点处出现过渡带，它主要是由窗函数频谱的主瓣引起的，过渡带的宽度取决于窗函数主瓣的宽度。另外滤波器还会在通带和阻带内产生波纹。

窗函数设计的思想是根据给定的滤波器技术指标，选择滤波器的长度和窗函数，使其具有最窄宽度的主瓣和最小的旁瓣。核心是从给定的频率特性通过加窗确定有限长单位脉冲响应序列h(n)。

由于FIR的单位冲激响应是有限长的，因此滤波器是稳定的，任何非因果有限长序列，只要经过一定的延时，都能变成因果的有限长序列，因而总能用因果系统来实现。FIR的系统函数是N-1次多项式，它在Z平面有N-1个零点，原点是N-1阶重极点，因此，H(z)是稳定的。

利用窗函数设计滤波器，是让待设计的滤波器去逼近理想特性，理想低通滤波器的频率特性应该是振幅特性在通带内为1，阻带内为0，在通带内的相位特性与w成线性关系。设计步骤：

1.确定数字滤波器的性能要求，确定滤波器的临界频率，滤波器的长度。2.根据性能要求，合理选择单位脉冲响应，从而确定理想频率响应的幅频特性和相频特性

3.选择适合的窗函数起初所需设计的FIR滤波器的单位脉冲响应。若不满足，可以改变传函数的形式或者是长度N

程序所用子函数：

function[db,mag,pha,grd,w]=freqz_m(b,a)functionhd=ideal_lp(wc,M);FIR数字滤波器的设计

用窗函数法设计一个线性相位FIR数字低通滤波器。技术指标为：通带截止角频率ωp=0.2π，阻带截止角频率ωs=0.3π，通带衰减Rp≤1dB，阻带衰减Rs≥40dB。绘出滤波器的幅频特性曲线和相频特性曲线，判断设计是否符合要求。

根据相同的滤波器要求，选用不同的窗函数进行设计，比较各种窗函数对FIR数字滤波器频率特性的影响。1、设计方案：

本实验采用了哈明窗、凯泽窗、布莱克曼窗这三种窗函数，对FIR滤波器进行设计。根据实验设计目标，本实验采用哈明窗可以最接近要求，其次是布莱克曼窗、凯泽窗。

实验原理

实验内容

第9页

实验结果及分析

2、设计步骤：

（1）给出所要设计的FIR数字滤波器的技术指标，如通带截止频率、阻带截止频率、通带衰减、阻带衰减。

（2）根据允许的过渡带宽度及阻带衰减，初步选择窗函数和N值。（3）若采用理想低通逼近，则求出理想低通的冲击响应hd（n）。（4）将hd（n）与窗函数相乘得FIR数字滤波器的冲激响应h（n）。（5）计算FIR数字滤波器的频率响应，并验证是否达到指标。3.源程序：哈明窗：

wp=0.2*pi;ws=0.3*pi;N=51;

n=[0:1:N-1];wc=(ws+wp)/2;

hd=ideal_lp(wc,N);w_ham=(hamming(N))';h=hd.*w_ham

[db,mag,pha,grd,w]=freqz_m(h,[1]);delta_w=2*pi/1000;

Rp=-(min(db(1:1:wp/delta_w+1)))

As=-round(max(db(wp/delta_w+1:1:501)))

subplot(2,2,1);stem(n,hd);title('理想冲击响应')axis([0N-1-0.10.3]);xlabel('n');ylabel('hd(n)')subplot(2,2,2);stem(n,w_ham);title('Hamming窗')axis([0N-1-0.11.1]);xlabel('n');ylabel('w(n)')subplot(2,2,3);stem(n,h);title('实验冲激响应')axis([0N-1-0.10.3]);xlabel('n');ylabel('h(n)')

subplot(2,2,4);plot(w/pi,db);title('幅度响应（dB）');gridaxis([00.8-1000]);

xlabel('以/pi为单位的频率');ylabel('对数幅度/dB');凯泽窗：

wp=0.2*pi;ws=0.3*pi;As=40;N=65;

n=[0:1:N-1];wc=(ws+wp)/2;

hd=ideal_lp(wc,N);w_bla=(blackman(N))';h=hd.*w_bla;

[db,mag,pha,grd,w]=freqz_m(h,[1]);delta_w=2*pi/1000;

Rp=-min(db(wp/delta_w+1:1:wp/delta_w))As=-round(max(db(ws/delta_w+1:1:501)))

第10页

subplot(2,2,1);plot(w/pi,db);title('幅度相应（dB）');grid;axis([01-15010]);

xlabel('以\pi为单位的频率');ylabel('对数幅度/dB');subplot(2,2,2);plot(w/pi,pha);title('相位响应');grid;

axis([01-44]);xlabel('以\pi为单位的频率');ylabel('相位');布莱克曼窗：

wp=0.2*pi;ws=0.3*pi;As=40;N=65;

n=[0:1:N-1];wc=(ws+wp)/2;

hd=ideal_lp(wc,N);w_bla=(blackman(N))';h=hd.*w_bla;

[db,mag,pha,grd,w]=freqz_m(h,[1]);delta_w=2*pi/1000;

Rp=-min(db(wp/delta_w+1:1:wp/delta_w))As=-round(max(db(ws/delta_w+1:1:501)))subplot(2,2,1);plot(w/pi,db);title('');grid;

axis([01-15010]);

xlabel('');ylabel('/dB');subplot(2,2,2);plot

子函数：

function[db,mag,pha,grd,w]=freqz_m(b,a)[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);

db=20*log10((mag+eps)/max(mag));pha=angle(H);

grd=grpdelay(b,a,w);functionhd=ideal_lp(wc,M);alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha+eps;

hd=sin(wc*m)./(pi*m);一、程序运行图像：1.哈明窗

第11页

2、凯泽窗

3、布莱克曼窗

第12

页

实验结果分析

上面三幅图像依次是哈明窗，凯泽窗，布莱克曼窗。N越大，滤波器的过渡带就越窄，主瓣高度与第一旁瓣高度的差值基本不变。旁瓣的个数越来越多，旁瓣的宽度随N的加大而减小。

由此可以知道，当滤波器阶数N相同时，实验中三个窗函数中哈明窗阻带衰减最小，其次是布莱克曼窗，凯泽窗的阻带最小衰减最大。

对于主瓣宽度方面，布莱克曼比哈明窗的主瓣要宽，主瓣高度与第一旁瓣高度的差值方面哈明窗最小，布莱克曼窗高度差最大。

对于旁瓣衰减的速度，哈明窗比较缓慢，布莱克曼窗最快。

旁瓣幅度方面，矩形窗最大，汉宁窗和哈明窗较大，布莱克曼窗较小。教材上讲，选用窗函数要求有二：旁瓣高度尽可能小，以使能量尽量集中在主瓣；主板宽度尽可能窄，以获得尽可能陡的过渡带。由上面的几个图形可以看出，这两点难以同时满足，当选用主瓣宽度较窄时，虽然得到的幅频特性较陡峭，但通带、阻带波动会明显增加；当选用较低的旁瓣幅度时，虽然得到的幅频特性较平缓匀滑，但过渡带变宽。在实际设计中要选用合适的参数，以使滤波器的效益达到最高。

这次课程设计持续了一周时间。虽然考试早已经结束，但通过这次实验对数字信号处理这门课程又有了更多的收获和更好的理解。

实验一主要研究了FFT的应用。FFT是数字信号处理的重要的一部分，它在能得到正确结果的同时，大大提高了运算速度，提高了工作效率。如实验一所示，FFT不仅可以利用循环卷积来计算线性卷积，也可以对信号进行频谱分析。但是FFT算法过程中还有很多需要注意的地方。由此可以看出，掌握FFT算法是非常重要的。

实验二主要研究了用MATLAB完成对IIR滤波器的设计。通过实验，我掌握了双线性变换法的设计原理和设计步骤，也初步掌握了利用MATLAB设计IIR滤波器的过程，尤其是对子函数的具体含义有了深刻的理解。

实验三主要研究了用MATLAB设计FIR滤波器，通过选用不同的窗函数，得到不同的结果来判决如何选择正确、准确且满足设计目标的窗函数。通过实验所得图像可以知道，滤波器的要求，即旁瓣高度尽可能小，主板宽度尽可能窄，这两点难以同时实现。所以在设计滤波器的时候要综合考虑两个因素，以得到最好、最符合要求、效益最高的滤波器。

本次试验虽然程序在课本、参考书和网上都有，但运行起来还是不可避免的出了一些问题，比如IIR滤波器一开始就无法运行来，但后来经过老师的指导和同学的帮助得到了正确的答案。是自己的粗心大意还有对程序的不理解，所以才出了差错，以后一定要仔细认真，高效率的完成实验。

总的来说，这次实验比较顺利，基本完成了实验目标，理论和实践的结合让我对知识有了更加深刻的理解和认识，受益匪浅。但个别知识点还有弄不清楚的地方，虽然这门课程结束了，但它对以后专业课还是有很大影响的，所以我会继续好好研究这些不懂的地方，争取把这一门课都掌握的清清楚楚。最后十分感谢所有老师的耐心指导！

实验总结

第13页