1、哈工大数字信号处理实验报告实验一： 用FFT作谱分析实验目的： (1) 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法， 所以FFT的运算结果必然满足DFT的基本性质)。 (2) 熟悉FFT算法原理和FFT子程序的应用。 (3) 学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。实验原理： DFT的运算量： 一次完整的DFT运算总共需要次复数乘法和复数加法运算，因而直接计算DFT时，乘法次数和加法次数都和成正比，当N很大时，运算量很客观的。例如，当N=8时，DFT运算需64位复数乘法，当N=1024时
2、，DFT运算需1048576次复数乘法。而N的取值可能会很大，因而寻找运算量的途径是很必要的。 FFT算法原理： 大多数减少离散傅里叶变换运算次数的方法都是基于的对称性和周期性。 （1）对称性 （2）周期性 由此可得 这样： 1.利用第三个方程的这些特性，DFT运算中有些项可以合并； 2.利用的对称性和周期性，可以将长序列的DFT分解为短序列的DFT。 前面已经说过，DFT的运算量是与成正比的，所以N越小对计算越有利，因而小点数序列的DFT比大点数序列的DFT运算量要小。 快速傅里叶变换算法正是基于这样的基本思路而发展起来的，她的算法基本上可分成两大类，即按时间抽取法和按频率抽取法。 我们最常
3、用的是的情况，该情况下的变换成为基2快速傅里叶变换。 完成一次完整的FFT计算总共需要次复数乘法运算和次复数加法运算。很明显，N越大，FFT的优点就越突出。实验步骤 (1) 复习DFT的定义、 性质和用DFT作谱分析的有关内容。 (2) 复习FFT算法原理与编程思想， 并对照DIT-FFT运算流图和程序框图， 读懂本实验提供的FFT子程序。 (3) 编制信号产生子程序， 产生以下典型信号供谱分析用： (4) 编写主程序。 (5) 按实验内容要求， 上机实验， 并写出实验报告。实验程序与结果: （1）对x1(n)进行FFT变换（N=8和N=16）clear all;N=4;%读入长度x=1,1,
4、1,1;%调用信号产生子程序产生实验信号n=0:N-1;subplot(3,1,1);stem(n,abs(x);%调用绘图子程序(函数)绘制时间序列波形图title(原时间序列);N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=8);N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;
5、subplot(3,1,3);stem(n,fftshift(abs(y2);hold on;plot(n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=16); 结果图形实验误差分析： 理论FFT频谱为一个抽样信号，由图形看出，随着采样率的提高，得到的FFT频谱分辨率就越高，当N趋于无限时频谱包络接近理论的抽样函数。 （2）对x2(n)进行FFT变换（N=8和N=16）clear all;x=1,2,3,4,4,3,2,1;%调用信号产生子程序产生实验信号N=8;%读入长度n=0:N-1;subplot(3,1,1);stem(n,ab
6、s(x);%调用绘图子程序(函数)绘制时间序列波形图title(原时间序列);N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=8);N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2);hold on;plot(
7、n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=16); 结果图形实验误差分析： 理论FFT频谱为一个抽样函数平方的信号，由图形看出，随着采样率的提高，得到的FFT频谱分辨率就越高，当N趋于无限时频谱包络接近理论的抽样函数。 （3）对x3(n)进行FFT变换（N=8和N=16）clear all;x=4,3,2,1,1,2,3,4;%调用信号产生子程序产生实验信号N=8;%读入长度n=0:N-1;subplot(3,1,1);stem(n,abs(x);%调用绘图子程序(函数)绘制时间序列波形图title(原时间序列);N=8;%读入长
8、度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=8);N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2);hold on;plot(n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线ti
9、tle(N=16); 结果图形实验误差分析： 由图形看出，随着采样率的提高，得到的FFT频谱分辨率就越高，但对于这个函数来说，N仍然过小，因此幅频特性曲线与理论结果有较大差距。 （4）对x4(n)进行FFT变换（N=8和N=16）clear all;n=0:15;%读入长度x=cos(pi/4*n);%调用信号产生子程序产生实验信号subplot(3,1,1);stem(n,abs(x);%调用绘图子程序(函数)绘制时间序列波形图title(cos(pi/4*n);N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);
10、stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=8);N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2);hold on;plot(n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=16); 结果图形：实验误差分析： 这个函数的理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，
11、FFT都可以与理论图形一样的曲线。当N=8,16时，满足抽样定理。 （5）对x5(n)进行FFT变换（N=8和N=16）clear all;n=0:15;%读入长度x=sin(pi/8*n);%调用信号产生子程序产生实验信号subplot(3,1,1);stem(n,abs(x);%调用绘图子程序(函数)绘制时间序列波形图title(sin(pi/8*n);N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1)
12、,r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=8);N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2);hold on;plot(n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=16); 结果图形：实验误差分析： 这个函数的理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=8时，不满足抽样定理，因此产生了混叠失真，当N=8时，满足抽样定理，因
13、此可得到与理论曲线一致的图形。 （6）对x6(n)进行FFT变换（N=16、N=32和N=64）clear all;n=0:64;%读入长度x=cos(pi/8*n)+cos(pi/4*n)+cos(pi/16*5*n);%调用信号产生子程序产生实验信号subplot(4,1,1);stem(n,abs(x);hold on;plot(n,abs(x),r-);%调用绘图子程序(函数)绘制时间序列波形图title(cos(pi/8*n)+cos(pi/4*n)+cos(pi/16*5*n);N=16;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;s
14、ubplot(4,1,2);stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=16);N=32;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(4,1,3);stem(n,fftshift(abs(y2);hold on;plot(n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=32);N=64;%读入长度y3=fft(x,N);%调用FFT子程序（
15、函数）计算信号的DFTn=0:N-1;subplot(4,1,4);stem(n,fftshift(abs(y3);hold on;plot(n,fftshift(abs(y3),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=64); 结果图形：实验误差分析： 这个函数的理论FFT抽样频谱为有三个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=16时，不满足抽样定理，因此产生了混叠失真，当N=32,64时，满足抽样定理，因此可得到与理论曲线一致的图形。 （7）对x(n)=x4(n)+x5(n)进行FFT变换（N=8和N=16）cl
16、ear all;n=0:15;%读入长度x4=cos(pi/4*n);x5=sin(pi/8*n);x=x4+x5;%调用信号产生子程序产生实验信号subplot(3,1,1);stem(n,abs(x);hold on;plot(n,abs(x),r-);%调用绘图子程序(函数)绘制时间序列波形图title(cos(pi/4*n)+sin(pi/8*n);N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1)
17、,r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=8);N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2);hold on;plot(n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=16); 结果图形：实验误差分析： 这个函数的理论FFT抽样频谱为有两个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=8时，不满足抽样定理，因此产生了混叠失真，当N=1
18、6时，满足抽样定理，因此可得到与理论曲线一致的图形。 （8）对x(n)=x4(n)+jx5(n)进行FFT变换（N=8和N=16）clear all;n=0:15;%读入长度x4=cos(pi/4*n);x5=sin(pi/8*n);x=x4+x5*j;%调用信号产生子程序产生实验信号subplot(3,1,1);stem(n,abs(x);hold on;plot(n,abs(x),r-);%调用绘图子程序(函数)绘制时间序列波形图title(cos(pi/4*n)+sin(pi/8*n)*j);N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-
19、1;subplot(3,1,2);stem(n,fftshift(abs(y1);hold on;plot(n,fftshift(abs(y1),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=8);N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2);hold on;plot(n,fftshift(abs(y2),r-);%调用绘图子程序(函数)绘制|X(k)|曲线title(N=16); 结果图形：实验误差分析： 这个函数的理论FFT抽样频谱为有
20、两个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=8时，不满足抽样定理，因此产生了混叠失真，当N=16时，满足抽样定理，因此可得到与理论曲线一致的图形。实验思考题 (1) 在N=8时，x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢? 答：在N=8的时候，x3(n)只是x2(n)以它的长度8为周期，将其延拓成周期 序列然后加以移位，最后截取主值区间（n=0到8）的序列值，因此x3(n)是x2(n) 进行循环位移后的结果。由于序列的循环位移性质 因此，循环位移只影响到DFT后的相频特性，而不影响幅频特性，因此x2(n) 和x3(n)的幅
21、频特性会相同。 若N=16，此时，x2(n)和x3(n)做DFT即为它们分别进行末尾补零后再 进行的DFT，则此时两个经过补零以后的序列就不满足循环位移的性质，因 此x2(n)和x3(n)的幅频特性就会发生变化。 (2) 如果周期信号的周期预先不知道， 如何用FFT进行谱分析? 答：若预先不知道周期信号的周期，应先适当截取M点进行FFT，再将截 取的长度扩大1倍重新截取，比较二者结果，若二者的差别满足分析误差的 要求，就可以近似得到该信号的频谱，假若不满足误差要求则继续加倍截取 的长度进行FFT，直到结果满足误差要求为止。实验二： 用窗函数法设计FIR数字滤波器实验目的： (1)熟悉矩形窗、海
22、宁窗、汉明窗和布莱克曼窗。 (2) 掌握用上述窗函数法设计FIR数字滤波器的原理和方法。 (3) 熟悉线性相位FIR数字滤波器特性。 (4) 了解各种窗函数对滤波特性的影响。 实验原理和方法： 如果所希望的滤波器的理想频率响应函数为Hd(e j)， 则其对应的单位脉冲响应为 用窗函数w(n)将hd(n)截断， 并进行加权处理， 得到: h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列， 其频率响应函数H(e j)为 如果要求线性相位特性， 则h(n)还必须满足： 根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如，
23、 要设计线性相位低通特性， 可选择h(n)=h(N-1-n)一类， 而不能选h(n)=-h(N-1-n)一类。实验内容及步骤： (1) 复习用窗函数法设计FIR数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。 (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
24、;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(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
25、）函数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(幅度曲线和相频曲线（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;
26、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(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
27、;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）函数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)
28、;%画幅度曲线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) 结果图像： （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;%
29、计算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(w,unwrap(phase(H1);%画相频曲线gridxlabel(w/rad)n=0:1:32;wH=0.54
30、-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);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
