哈工大数字信号处理报告.doc

1、数字信号处理实验报告课程名称： 数字信号处理 院 系： 班 级： 姓 名： 学 号： 指导教师： 实验时间： 2013 年11月 实验一： 用FFT作谱分析一、 实验目的1进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质)。2熟悉FFT算法原理和FFT子程序的应用。3学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。二、 实验步骤1复习DFT的定义、性质和用DFT作谱分析的有关内容。复习FFT算法原理与编程思想，并对照DITFFT运算流图和程序框图，2读懂本

2、实验提供的FFT子程序。3编制信号产生子程序，产生以下典型信号供谱分析： ， ， ， 其它 ， ， 0 ， 其它 应当注意，如果给出的是连续信号，则首先要根据其最高频率确定采样速率以及由频率分辨率选择采样点数N，然后对其进行软件采样(即计算，)，产生对应序列。对信号，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。4 编写主程序下图给出了主程序框图，供参考。本实验提供FFT子程序和通用绘图子程序。开始读入长度N调用信号产生子程序产生实验信号调用绘图子程序(函数)绘制时间序列波形图 调用FFT子程序(函数)计算

3、信号的DFT调用绘图子程序(函数)绘制曲线结束主程序框图三、 实验结果直接运行程序，按照实验内容及程序提示键入18，分别对及、进行谱分析。输出的波形及其8点DFT和16点DFT，的16点、32点和64点采样序列及其DFT。4.实验程序word文档 可自由复制编辑X1 的实验程序x1=1 1 1 1 ;y11=fft(x1,8);y12=fft(x1,16);subplot(2,3,1);stem(0:3,x1);title(函数x1的图像);box onsubplot(2,2,2);stem(0:7,abs(y11);title(N=8时，x1的DFT图像);box onsubplot(2,2

4、,4);stem(0:15,abs(y12);title(N=16时，x1的DFT图像);box onX2 的实验程序xa=1:1:4;xb=4:-1:1;x2=xa xb;y21=fft(x2,8);y22=fft(x2,16);subplot(2,2,1);stem(0:7,x2);title(函数x2的图像);subplot(2,2,2);box onstem(0:7,abs(y21);title(N=8时，x2的DFT图像);box onsubplot(2,2,4);stem(0:15,abs(y22);title(N=16时，x1的DFT图像);box onX3的实验程序xa=4:-

5、1:1;xb=1:1:4;x3=xa xb;y31=fft(x2,8);y32=fft(x2,16);subplot(2,2,1);stem(0:7,x3);title(函数x3的图像);subplot(2,2,2);stem(0:7,abs(y31);title(N=8时，函数x3的DFT图像);subplot(2,2,3);stem(0:15,abs(y32);title(N=16时，函数x3的DFT图像);title(N=16时，x2的DFT图像);box onX4 X5 X7 X8的实验程序subplot(441)n=0:1:7x4=cos(0.25*pi*n)stem(n,x4)xl

6、abel(n)ylabel(x4(n)title(x4(n)subplot(442)n=0:1:15x44=cos(0.25*pi*n)stem(n,x44)xlabel(n)ylabel(x44(n)title(x44(n)subplot(443)n=0:1:7stem(n,abs(fft(x4)xlabel(n)ylabel(|X4(k)|)title(8 point FFT)subplot(444)n=0:1:15stem(n,abs(fft(x44)xlabel(n)ylabel(|X44(k)|)title(16 point FFT)subplot(445)n=0:1:7x5=sin

7、(0.125*pi*n)stem(n,x5)xlabel(n)ylabel(x5(n)title(x5(n)subplot(446)n=0:1:15x55=sin(0.125*pi*n)stem(n,x55)xlabel(n)ylabel(x55(n)title(x55(n)subplot(447)n=0:1:7stem(n,abs(fft(x5)xlabel(n)ylabel(|X5(k)|)title(8 point FFT)subplot(448)n=0:1:15stem(n,abs(fft(x55) subplot(449)n=0:1:7x7=cos(0.25*pi*n)+sin(0.

8、125*pi*n)stem(n,x7)xlabel(n)ylabel(x7(n)title(x7(n)subplot(4,4,10)n=0:1:15x77=cos(0.25*pi*n)+sin(0.125*pi*n)stem(n,x77)xlabel(n)ylabel(x77(n)title(x77(n)subplot(4,4,11)n=0:1:7stem(n,abs(fft(x7)xlabel(n)ylabel(|X7(k)|)title(8 point FFT)subplot(4,4,12)n=0:1:15stem(n,abs(fft(x77)xlabel(n)ylabel(|X77(k)

9、|)title(16 point FFT) subplot(4,4,13)n=0:1:7x8=cos(0.25*pi*n)+j*sin(0.125*pi*n)stem(n,x8)xlabel(n)ylabel(x8(n)title(x8(n)subplot(4,4,14)n=0:1:15x88=cos(0.25*pi*n)+j*sin(0.125*pi*n)stem(n,x88)xlabel(n)ylabel(x88(n)title(x88(n)subplot(4,4,15)n=0:1:7stem(n,abs(fft(x8)xlabel(n)ylabel(|X8(k)|)title(8 poi

10、nt FFT)subplot(4,4,16)n=0:1:15stem(n,abs(fft(x88)xlabel(n)ylabel(|X88(k)|)title(16 point FFT) X6的实验程序subplot(231)n=0:16fs=64x6=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs)stem(n,x6)xlabel(n)ylabel(x6(n)title(x6(n)subplot(232)n=0:32fs=64x66=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs)stem(n,x66)h

11、old onxlabel(n)ylabel(x66(n)title(x66(n)subplot(233)n=0:64fs=64x666=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs)stem(n,x666)hold onxlabel(n)ylabel(x666(n)title(x666(n)subplot(234)stem(abs(fft(x6)hold onxlabel(n)ylabel(|X6(k)|)title(16 point FFT)subplot(235)stem(abs(fft(x66)hold onxlabel(n)ylabel(|

12、X66(k)|)title(32 point FFT)subplot(236)stem(abs(fft(x666)hold onxlabel(n)ylabel(|X666(k)|)title(64 point FFT)1及其8点和16点DFT2及其8点和16点DFT3及其8点和16点DFT4的8点和16点波形及其DFT5的8点和16点波形及其DFT6的16点、32点和64点采样序列波形及其DFT7计算并图示和及其DFT。程序自动计算并绘图验证DFT的共轭对称性。当N=16时，。即为的共轭对称分量，而是的共轭反对称分量。根据DFT的共轭对称性，应有以下结果：的8点和16点波形及其DFT绘出和的模

13、。它们正是图中16点的和。8计算并图示和及其DFT。程序自动计算并绘图验证DFT的共轭对称性的第二种形式：如果，则， 。其中，。的8点和16点DFT程序计算结果如下：及，正好与图中的16点及相同。及，正好与图中16点的及相同。四、 实验总结本实验主要是求、的DFT变换。其中是直接给出了离散序列，而、则是经过、运算得到的、。离散傅立叶变换可以看作是在时的Z变换，即表明的N点DFT是的Z变换在单位圆上的N点等间隔采样。离散傅立叶变换也可以看作在时的傅立叶变换，即表明可以看作的傅立叶变换在区间上的N点等间隔采样。五、 思考题1在N=8时，和的幅频特性会相同吗？为什么？N=16呢？答：N=8时两个的幅

14、频特性相同，因为其不为0的区间长度正好是8。N=16时两个的幅频特性不相同。2、如果周期信号的周期预先不知道，如何用FFT进行谱分析？ 答：可以先求出它的离散傅立叶变换，通过其DFT的图形来确定出原始序列的周期。实验二： 用窗函数法设计FIR数字滤波器实验目的：(1)熟悉矩形窗、海宁窗、汉明窗和布莱克曼窗。(2) 掌握用上述窗函数法设计FIR数字滤波器的原理和方法。(3) 熟悉线性相位FIR数字滤波器特性。(4) 了解各种窗函数对滤波特性的影响。 实验原理和方法：如果所希望的滤波器的理想频率响应函数为Hd(e j)， 则其对应的单位脉冲响应为用窗函数w(n)将hd(n)截断， 并进行加权处理，

15、 得到: h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列， 其频率响应函数H(e j)为如果要求线性相位特性， 则h(n)还必须满足：根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性， 可选择h(n)=h(N-1-n)一类， 而不能选h(n)=-h(N-1-n)一类。实验内容及步骤：(1) 复习用窗函数法设计FIR数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。(2) 编写程序。 编写能产生四种窗函数的子程序。 编写主程序。 实验程序及结果：（1）矩形窗clear all;n=0

16、:1:14;wR=ones(1,15);% 编写矩形窗hd=sin(0.25*pi*(n-7+eps)./(pi*(n-7+eps);%读入hd（n）函数h1=hd.*wR;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);plot(w,fftshift(20*log10(abs(H1);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB)title(幅度曲线和相频曲线（n=15）);n=0:N-1;w=2*pi/64*n;subplot(2,2,3);plot(

17、w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)clear all;n=0:1:32;wR=ones(1,33);% 编写矩形窗hd=sin(0.25*pi*(n-16+eps)./(pi*(n-16+eps);%读入hd（n）函数h1=hd.*wR;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,2);plot(w,fftshift(20*log10(abs(H1);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB)title(

18、幅度曲线和相频曲线（n=33）);n=0:N-1;w=2*pi/64*n;subplot(2,2,4);plot(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)结果图像：（2）汉宁窗clear all;n=0:1:14;wH=0.5*(1-cos(2*pi/14*n);% 编写汉宁窗hd=sin(0.25*pi*(n-7+eps)./(pi*(n-7+eps);%读入hd（n）函数h1=hd.*wH;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot

19、(2,2,1);plot(w,fftshift(20*log10(abs(H1);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB);title(幅度曲线和相频曲线（n=15）);n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,3);plot(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)n=0:1:32;wH=0.5*(1-cos(2*pi/32*n);% 编写汉宁窗hd=sin(0.25*pi*(n-16+eps)./(pi*(n-16+eps);%读入hd（n）函数

20、h1=hd.*wH;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,2);plot(w,fftshift(20*log10(abs(H1);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB)title(幅度曲线和相频曲线（n=33）);n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,4);plot(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)结果图像：

21、（3）海明窗：clear all;n=0:1:14;wH=0.54-0.46*cos(2*pi*n/(14+eps);% 编写海明窗hd=sin(0.25*pi*(n-7+eps)./(pi*(n-7+eps);%读入hd（n）函数h1=hd.*wH;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,1);plot(w,fftshift(20*log10(abs(H1);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB)title(幅度曲

22、线和相频曲线（n=15）);n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,3);plot(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)n=0:1:32;wH=0.54-0.46*cos(2*pi*n/(32+eps);% 编写海明窗hd=sin(0.25*pi*(n-16+eps)./(pi*(n-16+eps);%读入hd（n）函数h1=hd.*wH;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplo

23、t(2,2,2);plot(w,fftshift(20*log10(abs(H1);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB)title(幅度曲线和相频曲线（n=33）);n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,4);plot(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad) 结果图像 (4)布莱克曼窗n=0:1:14;wB=0.42-0.5*cos(2*pi/(14+eps)*n)+0.08*cos(4*pi/(14+eps)*n);% 编写布莱克曼窗hd

24、=sin(0.25*pi*(n-7+eps)./(pi*(n-7+eps);%读入hd（n）函数h1=hd.*wB;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,1);plot(w,fftshift(20*log10(abs(H1);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB)title(幅度曲线和相频曲线（n=15）);n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,3);plot

25、(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)n=0:1:32;wB=0.42-0.5*cos(2*pi/(32+eps)*n)+0.08*cos(4*pi/(32+eps)*n);% 编写布莱克曼窗hd=sin(0.25*pi*(n-16+eps)./(pi*(n-16+eps);%读入hd（n）函数h1=hd.*wB;%计算h（n）N=64;H1=fft(h1,N);%调用子程序计算H（k）n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,2);plot(w,fftshift(20*log10(abs(H1

26、);%画幅度曲线gridxlabel(w/rad)ylabel(20lg|H(jw)|/dB)title(幅度曲线和相频曲线（n=33）);n=0:N-1;w=2*pi/64*n;subplot(2,2,1);subplot(2,2,4);plot(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)结果图像：实验分析： 窗函数法特点：窗口法设计的主要优点是简单，使用方便。窗口函数大多有封闭的公式可循，性能、参数都已有表格、资料可供参考， 计算程序简便， 所以很实用。缺点是通带和阻带的截止频率不易控制。 窗口长度N和窗函数类型对滤波特性的影响：增加窗的长度，可以

27、减少窗的主瓣宽度，从而减少过渡带的带宽，但增加N并不能减少带内波动以及加大阻带。这两个指标只能从窗函数的形状上找解决方法。常用的6种窗口的基本参数及性能窗类型窗谱性能指标加窗后滤波器的性能指标最大旁瓣(dB)主瓣带宽阻带衰减(dB)过渡带宽矩形窗-132*2/N-210.9*2/N三角窗-254*2/N-252.1*2/Nhanning窗-314*2/N-443.1*2/Nhamming窗-414*2/N-533.3*2/NBlackman窗-576*2/N-745.5*2/NKaiser窗（=7.865）-57.07-805.0*2/N通过对窗函数选择，可以获得符合要求的阻带衰减指标。实验总

28、结及心得用窗函数法设计FIR滤波器的主要特点：设FT为希望逼近的频响特性函数,H() =FTh(n)为用窗函数法设计的实际滤波器的频响函数。通常取H()相应的理想频响特性作为。知识要点如下：1希望逼近的理想滤波器频响函数的表达式。因为数字滤波器一般要求设计成线性相位特性，所以必须满足上述线性相位FIR滤波器的频域特点。2熟悉各种常用窗函数的技术指标和加窗后对滤波特性的影响，这样才能根据设计指标正确选择窗函数类型及其长度N。3检验设计结果：窗函数法的设计结果单位脉冲响应。而检验一般在频域进行。所以要计算检验3dB截止频率和阻带最小衰减，其计算量相当大，必须用计算机进行。4熟悉窗函数设计法的特点：

29、设计过程简单.方便实用。但边界频率不易精确控制。所以设计完以后，必须检验结果。思考题(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器?写出设计步骤。答：a) 求出对应的数字频率。 通带截止频率 阻带截止频率b）求。设为理想线性相位滤波器 其中为线性相位所必须的位移，且满足。 c）求窗函数。由阻带衰减确定窗形状，由过度带宽确定N。 d）求。由窗函数表达式确定FIR滤波器的。 e）由求，检验各项指标是否满足要求。如不满足要求，则要改变N，或改变窗形状，然后重新计算。(2) 如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为和，试求理想带通的单位脉冲响应hd(n) .答：只需将一个高通滤波器器和一个低通滤波器相加，即可得到带通滤波器，具体设计过程与低通滤波器相同。具体步骤如下： a）理想带通滤波器的频率响应是 其中 b）求此滤波器的单位冲激响应，即 c）根据对通带、阻带衰减的要求以及对过渡带宽的要求，可选择窗函数及窗的点数N。由此可求得所需的线性相位带通滤波器的单位冲激相应 d）由求，检验各项指标是否满足要求。如不满足要求，则要改变N，或改变窗形状，然后重新计算。word文档 可自由复制编辑