哈工大数字信号处理实验报告

HarbinInstituteofTechnology实验报告课程名称：数字信号处理实验题目：用FFT作谱分析、用窗函数法设计数字滤波器院系：电子与信息工程学院班级：哈尔滨工业大学一、实验目的(1)进一步加深DFT算法原理和根本性质的理解。(2)熟悉FFT算法原理和FFT子程序的应用。(3)学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。二、实验步骤(1)复习DFT的定义、性质和用DFT作谱分析的有关内容。(2)复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图。(3)编制信号产生程序，产生以下典型信号供谱分析用：(4)按实验内容要求，上机实验，并写出实验报告。三、实验原理1、FFT产生背景。离散傅立叶变换从理论上解决了傅立叶变换应用于实际的可能性，但假设直接按DFT公式计算，运算量太大〔与N2成比例〕。快速傅立叶变换〔FFT〕是离散傅立叶变换的快速算法，它大大减少了离散傅立叶变换的运算次量，一般可缩短一、二故数量级，且N越大改善越明显，从而使DFT的运算在实际工作中才真正得到广泛应用。2、FFT原理。减少离散傅立叶变换运算次数的方法基于的对称性和周期性。〔1〕对称性〔2〕周期性由此可得这样可以将DFT运算中有些项进行合并，利用的对称性和周期性，将长序列的DFT分解为短序列的DFT。3、FFT运算量。在采用按时间抽取法的情况下，当N=2M时共有M级蝶形，每级都有N/2个蝶形运算。计算每个蝶形需要一次复数乘法和两次复数加法，因而每级运算需N/2次复数乘法和N次复数加法，这样M级运算总共需要：复数乘法次数：mF=N/2.M=N/2.log2N复数加法次数：aF=NM=Nlog2N在一般情况下，复数乘法所需时间比复数加法多，因此以复数乘法为例将DFT运算量与FFT运算量进行比照。直接进行DFT运算复数乘法次数为N2次，利用FFT运算复数乘法为N/2.log2N次。计算量之比为：四、实验结果、源程序及误差分析1、实验程序：N1=8;N2=16;x1=[1,1,1,1];y11=fft(x1,8);y12=fft(x1,16)subplot(2,2,1);stem(0:3,x1);title('函数x1时域图像');boxoffsubplot(2,2,2);stem(0:(N1-1),abs(y11));title('N=8，x1的DFT图像');boxoffsubplot(2,2,3);stem(0:(N2-1),abs(y12));title('N=16，x1的DFT图像');boxoff实验误差分析：理论的DFT函数应为抽样函数，由实验结果可以看出，随着FFT点数的增加，FFT图像包络函数越趋近于抽样函数。2、实验结果：实验程序：N1=8;N2=16;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);boxoffstem(0:(N1-1),abs(y21));title('N=8，x2的DFT图像');boxoffsubplot(2,2,4);stem(0:(N2-1),abs(y22));title('N=16，x2的DFT图像');boxoff实验误差分析：理论的DFT函数应为形如抽样函数的平方的函数，由实验结果可以看出，随着FFT点数的增加，FFT图像包络函数接近抽样函数的平方。3、实验结果：实验程序：N1=8;N2=16;xa=4:-1:1;xb=1:1:4;x3=[xa,xb];y31=fft(x3,8);y32=fft(x3,16);subplot(2,2,1);stem(0:7,x3);title('函数x3时域图像');subplot(2,2,2);boxoffstem(0:(N1-1),abs(y31));title('N=8时，x3的DFT图像');boxoffsubplot(2,2,3);stem(0:(N2-1),abs(y32));title('N=16时，x3的DFT图像');boxoff实验误差分析：随着N的增加，曲线接近理论值，具体的分析和2中的信号相同4、实验结果：实验程序：N1=8;N2=16;n=0:1:19x4=cos(0.25*pi*n);n=0:1:7;x411=cos(0.25*pi*n);n=8:1:15;x412=cos(0.25*pi*n);n=16:1:19;x413=cos(0.25*pi*n);y411=fft(x411,8);y412=fft(x412,8);y413=fft(x413,8);y41=y411+y412+y413;n=0:1:15;x421=cos(0.25*pi*n);n=16:1:19;x422=cos(0.25*pi*n);y421=fft(x421,16);y422=fft(x422,16);y42=y421+y422;subplot(2,2,1);stem(0:19,x4);title('函数x4时域图像');boxoffsubplot(2,2,2);stem(0:(N1-1),abs(y41));title('N=8,x4的DFT图像');boxoffsubplot(2,2,3);stem(0:(N2-1),abs(y42));title('N=16,x4的DFT图像');boxoff实验误差分析：理论的DFT函数为冲激函数，从实验结果来看，有两个明显的峰值可以看到，与理论值接近。5、实验结果：实验程序：n=0:1:19;N1=8;N2=16;x5=sin(0.125*pi*n);n=0:1:7;x511=sin(0.125*pi*n);n=8:1:15;x512=sin(0.125*pi*n);n=16:1:19;x513=sin(0.125*pi*n);y511=fft(x511,8);y512=fft(x512,8);y513=fft(x513,8);y51=y511+y512+y513;n=0:1:15;x521=sin(0.125*pi*n);n=16:1:19;x522=sin(0.125*pi*n);y521=fft(x521,16);y522=fft(x522,16);y52=y521+y522;subplot(2,2,1);stem(0:19,x5);title('函数x5时域图像');boxoffsubplot(2,2,2);stem(0:(N1-1),abs(y51));title('N=8,x5的DFT图像');boxoffsubplot(2,2,3);stem(0:(N2-1),abs(y52));title('N=16,x5的DFT图像');boxoff实验误差分析:理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=16时，不满足抽样定理，因此产生了混叠失真，当N=8时，满足抽样定理，因此可得到与理论曲线一致的图形。6、x4(n)+x5(n)实验结果：实验程序：n=0:1:19;N1=8;N2=16;x4=cos(0.25*pi*n);x5=sin(0.125*pi*n);x45=x4+x5;n=0:1:7;x4511=cos(0.25*pi*n)+sin(0.125*pi*n);n=8:1:15;x4512=cos(0.25*pi*n)+sin(0.125*pi*n);n=16:1:19;x4513=cos(0.25*pi*n)+sin(0.125*pi*n);y4511=fft(x4511,8);y4512=fft(x4512,8);y4513=fft(x4513,8);y451=y4511+y4512+y4513;n=0:1:15;x4521=cos(0.25*pi*n)+sin(0.125*pi*n);n=16:1:19;x4522=cos(0.25*pi*n)+sin(0.125*pi*n);y4521=fft(x4521,16);y4522=fft(x4522,16);y452=y4521+y4522;subplot(2,2,1);stem(0:19,x45);title('函数x45时域图像');boxoffsubplot(2,2,2);stem(0:(N1-1),abs(y451));title('N=8,x45的DFT图像');boxoffsubplot(2,2,3);stem(0:(N2-1),abs(y452));title('N=16,x45的DFT图像');boxoff实验误差分析：这个函数的理论FFT抽样频谱为有两个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=16时，不满足抽样定理，因此产生了混叠失真。7、x4(n)+j*x5(n)实验结果：实验程序：n=0:1:19;N1=8;N2=16;x4=cos(0.25*pi*n);x5=sin(0.125*pi*n);x45=x4+x5;n=0:1:7;x4511=cos(0.25*pi*n)+j*sin(0.125*pi*n);n=8:1:15;x4512=cos(0.25*pi*n)+j*sin(0.125*pi*n);n=16:1:19;x4513=cos(0.25*pi*n)+j*sin(0.125*pi*n);y4511=fft(x4511,8);y4512=fft(x4512,8);y4513=fft(x4513,8);y451=y4511+y4512+y4513;n=0:1:15;x4521=cos(0.25*pi*n)+j*sin(0.125*pi*n);n=16:1:19;x4522=cos(0.25*pi*n)+j*sin(0.125*pi*n);y4521=fft(x4521,16);y4522=fft(x4522,16);y452=y4521+y4522;subplot(2,2,1);stem(0:19,x45);title('函数x45时域图像');boxoffsubplot(2,2,2);stem(0:(N1-1),abs(y451));title('N=8,x45的DFT图像');boxoffsubplot(2,2,3);stem(0:(N2-1),abs(y452));title('N=16,x45的DFT图像');boxoff实验误差分析：理论FFT抽样频谱为两个单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=8时，不满足抽样定理，因此产生了混叠失真，当N=16时，满足抽样定理，因此可得到与理论曲线一致的图形。8、实验结果：实验程序：n=0:1:69;N1=16;N2=32;N3=64;fs=64;x6=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=0:1:15;x611=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=16:1:31;x612=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=32:1:47;x613=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=48:1:63;x614=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=64:1:69;x615=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);y611=fft(x611,N1);y612=fft(x612,N1);y613=fft(x613,N1);y614=fft(x614,N1);y615=fft(x615,N1);y61=y611+y612+y613+y614+y615;n=0:1:31;x621=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=32:1:63;x622=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=64:1:69;x623=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);y621=fft(x621,N2);y622=fft(x622,N2);y623=fft(x623,N2);y62=y621+y622;n=0:1:63;x631=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);n=64:1:69;x632=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);y631=fft(x631,N3);y632=fft(x632,N3);y63=y631+y632;subplot(2,2,1);stem(0:69,x6);title('函数x6时域图像');boxoffsubplot(2,2,2);stem(0:(N1-1),abs(y61));title('N=16,x6的DFT图像');boxoffsubplot(2,2,3);stem(0:(N2-1),abs(y62));title('N=32,x6的DFT图像');boxoffsubplot(2,2,4);stem(0:(N3-1),abs(y63));title('N=64,x6的DFT图像');boxoff实验误差分析：理论FFT抽样频谱为两个单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。当N=16时，不满足抽样定理，因此产生了混叠失真，当N=32,64时，满足抽样定理，因此可得到与理论曲线一致的图形。五、思考题(1)在N=8时，x2(n)和x3(n)的幅频特性会相同吗?为什么?N=16呢?答：从实验结果看出，N=8时，x2(n)和x3(n)的幅频特性相同。而N=16时，两者不同。原因是N=8时，x3(n)只是以长度8为周期，将其延拓称周期序列然后加以移位，最后截取主值区间的序列值。因此x3(n)是x2(n)进行循环位移后的结果。移位后只影响DFT相频特性而不影响幅频特性。而当N=16时，两个序列不满足循环位移的性质，因此幅频特性发生变化。(2)如果周期信号的周期预先不知道，如何用FFT进行谱分析?答：假设该信号是模拟信号，用FFT进行频谱分析时，首先对信号进行采样，使之变成离散信号，然后用FFT对离散信号进行频谱分析。由采样定理，采样频率fs应当大于信号最高频率的2倍。如果信号的周期预先不知道，可以选择一个较大的采样频率，用试验法逐步增大采样频率，观察频谱特性。一、实验目的(1)熟悉矩形窗、汉宁窗、海明窗和布莱克曼窗。(2)掌握用上述窗函数法设计FIR数字滤波器的原理和方法。(3)熟悉线性相位FIR数字滤波器特性。(4)了解各种窗函数对滤波特性的影响。二、实验原理与方法如果所希望的滤波器的理想频率响应函数为Hd(ejω)，那么其对应的单位脉冲响应为用窗函数w(n)将hd(n)截断，并进行加权处理，得到:h(n)作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数H(ejω)为如果要求线性相位特性，那么h(n)还必须满足：根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性，可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。三、实验内容与步骤(1)复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。(2)编写程序，其中幅度特性要求用dB表示。备注：不许用freqz()这个函数设设画图时，用20log|H(k)|打印幅度特性。第k点对应的频率为为使曲线包络更接近H(ejw)的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h(n)末尾补零的方法，使长度变为64，再进行64点DFT，那么可得到更精确的幅度衰减特性曲线。四、上机实验内容1、用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率Wc=rad。窗口N=15，33。要求在两种窗口长度情况下，分别求出h(n)，打印出相应的幅频特性和相频曲线。观察3db带宽和20db带宽，总结窗口N对滤波特性的影响。设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数即五、实验结果1、矩形窗，N=152、汉宁窗，N=153、海明窗，N=154、布莱克曼窗，N=155、矩形窗，N=336、汉宁窗，N=337、海明窗，N=338、布莱克曼窗，N=33分析：通过对上述四种窗函数对同一低通滤波器的幅频特性分析可以看出，在N一定的条件下矩形窗、汉宁窗、哈明窗、布莱克曼窗的过渡带愈来愈宽，阻带衰减越来愈大，他们都是牺牲过渡带的宽度来减少阻带的波动，增大阻带的衰减，减小吉布斯现象。实验程序：clearall;N=input('请输入窗口长度N\nN=');while(N)k=input('请选择窗函数类型：1矩形窗，2汉宁窗，3海明窗，4布莱克曼窗\nk=');ifk==1W=boxcar(N);elseifk==2W=hanning(N);elseifk==3W=hamming(N);elseifk==4W=blackman(N);endn=0:(N-1);wc=pi/4;a=n-(N-1)/2+eps;hd=sin(wc*a)./(pi*a);h=hd.*(W)';N2=1024;H=fft(h,N2);mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=unwrap(angle(H));%得到相位w=0:2*pi/N2:2*pi/N2*(N2-1);figuresubplot(2,2,1);n=0:(N-1);stem(n,h);axis([0,N-1,-0.1,0.5]);xlabel('n');ylabel('h(n)');ifN==15&&k==1title('矩形窗N=15时的h(n)');elseifN==15&&k==2title('汉宁窗N=15时的h(n)');elseifN==15&&k==3