哈工大数字信号处理报告

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

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

3、验信号调用绘图子程序(函数)绘制时间序列波形图 调用FFT子程序(函数)计算信号的DFT调用绘图子程序(函数)绘制曲线结束主程序框图实验结果直接运行程序，按照实验内容及程序提示键入18，分别对及、进行谱分析。输出的波形及其8点DFT和16点DFT，的16点、32点和64点采样序列及其DFT。4.实验程序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

4、图像);box onsubplot(2,2,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图像

5、);box onX3的实验程序xa=4:-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

6、.25*pi*n)stem(n,x4)xlabel(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)subp

7、lot(445)n=0:1:7x5=sin(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=

8、cos(0.25*pi*n)+sin(0.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)xl

9、abel(n)ylabel(|X77(k)|)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)ylabe

10、l(|X8(k)|)title(8 point FFT)subplot(4,4,16)n=0:1:15stem(n,abs(fft(x88)xlabel(n)ylabel(|X88(k)|)title(16 point FFT) X6的实验程序PAGE PAGE 25subplot(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*

11、n/fs)+cos(20*pi*n/fs)stem(n,x66)hold 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

12、(fft(x66)hold onxlabel(n)ylabel(|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 = 1 * GB3 计算并图示和及其DFT。 = 2 * GB3 程序自动计算并绘图验证DFT的共轭对称性。当N=16时，。即为的共

13、轭对称分量，而是的共轭反对称分量。根据DFT的共轭对称性，应有以下结果：的8点和16点波形及其DFT绘出和的模。它们正是图中16点的和。8 = 1 * GB3 计算并图示和及其DFT。 = 2 * GB3 程序自动计算并绘图验证DFT的共轭对称性的第二种形式：如果，则， 。其中，。的8点和16点DFT程序计算结果如下：及，正好与图中的16点及相同。及，正好与图中16点的及相同。实验总结本实验主要是求、的DFT变换。其中是直接给出了离散序列，而、则是经过、运算得到的、。离散傅立叶变换可以看作是在时的Z变换，即表明的N点DFT是的Z变换在单位圆上的N点等间隔采样。离散傅立叶变换也可以看作在时的傅立

14、叶变换，即表明可以看作的傅立叶变换在区间上的N点等间隔采样。五、 思考题1在N=8时，和的幅频特性会相同吗？为什么？N=16呢？答：N=8时两个的幅频特性相同，因为其不为0的区间长度正好是8。N=16时两个的幅频特性不相同。2、如果周期信号的周期预先不知道，如何用FFT进行谱分析？ 答：可以先求出它的离散傅立叶变换，通过其DFT的图形来确定出原始序列的周期。实验二： 用窗函数法设计FIR数字滤波器实验目的：(1)熟悉矩形窗、海宁窗、汉明窗和布莱克曼窗。(2) 掌握用上述窗函数法设计FIR数字滤波器的原理和方法。(3) 熟悉线性相位FIR数字滤波器特性。(4) 了解各种窗函数对滤波特性的影响。

15、实验原理和方法：如果所希望的滤波器的理想频率响应函数为Hd(e j)， 则其对应的单位脉冲响应为用窗函数w(n)将hd(n)截断， 并进行加权处理， 得到: h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列， 其频率响应函数H(e j)为如果要求线性相位特性， 则h(n)还必须满足：根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性， 可选择h(n)=h(N-1-n)一类， 而不能选h(n)=-h(N-1-n)一类。实验内容及步骤：(1) 复习用窗函数法设计FIR数字滤波器一节内容， 阅读

16、本实验原理， 掌握设计步骤。(2) 编写程序。 编写能产生四种窗函数的子程序。 编写主程序。 实验程序及结果：（1）矩形窗clear all;n=0: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

17、|H(jw)|/dB)title(幅度曲线和相频曲线（n=15）);n=0:N-1;w=2*pi/64*n;subplot(2,2,3);plot(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,fft

18、shift(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,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（

19、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(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)n=0:1:32;wH=0.5*(1-

20、cos(2*pi/32*n);% 编写汉宁窗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);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(

21、2,2,1);subplot(2,2,4);plot(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)结果图像：（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,fftshi

22、ft(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.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

23、（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) 结果图像 (4)布莱克曼窗n=0:1

24、:14;wB=0.42-0.5*cos(2*pi/(14+eps)*n)+0.08*cos(4*pi/(14+eps)*n);% 编写布莱克曼窗hd=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)tit

25、le(幅度曲线和相频曲线（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;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

26、;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)结果图像：实验分析： 窗函数法特点：窗口法设计的主要优点是简单，使用方便。窗口函数大多有封闭的公式可循，性能、参数都已有表格、资料可供参考， 计算程序简便， 所以很实用。缺点是通带和阻带的截止频率不易控制。 窗口长度N和窗函数类型对滤波特性的影响：增加窗的长度，可以减少窗的主瓣宽度，从而减少过渡带的带宽，但增加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/NBla