数字信号处理课程设计报告实验三.doc

课程设计报告课 程: 数字信号处理课程设计 学 院: 信息工程学院 专 业: 信息工程 学 号: 学生姓名: 教师姓名: 2019年10月13日实验三：FFT频谱分析及应用一、实验目的：（一）通过实验，加深对FFT的理解，熟悉FFT子程序。（二）熟悉用FFT对典型信号进行频谱分析的方法。二、实验原理与方法： 在各种信号序列中，有限长序列占有重要地位。对有限长序列，可以利用离散傅里叶变换（DFT）进行分析。DFT不但可以很好的反应序列频谱特性，而且易于用快速算法（FFT）在计算机上实现。设序列为x(n)，长度为N，其DFT定义为：Xk=n=0N-1xnWNnk，k=0,N-1反变换为xn=1Nk=0N-1X(k)WN-kn，n=0, ,N-1 有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列傅里叶变换的等距采样，因此可以用于序列的谱分析。FFT是DFT的一种快速算法，是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。常用的FFT是以2为基数，其长度为N=2M。它的效率高，程序简单，使用方便。当要变换的序列长度不等于2的整数幂次时，为了使用以2为基数的FFT，可以使用末尾补零的方法，使其长度为2的整数次方。 在MATLAB信号处理工具箱中的函数为fft(x,N)，可用于序列x(n)的N点快速傅里叶变换。经函数fft求得的序列一般是复序列，通常要求其幅值和相位。MATLAB中提供了求复数的幅值和相位函数：abs、angle。三、实验内容：（一）模拟信号xt=2cos4t+5cos(8t)，以0.01n进行采样，其中n=0,N-1：求N=40点FFT的幅度频谱，从图中能否观察出信号的2个频率分量？提高采样点数，如N=128，再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的2个模拟频率和数字频率各为多少？FFT频谱分析结果与理论上是否一致？解：MATLAB程序：N=40;n=0:N-1;t=0.01*n;x=2*cos(4*pi*t)+5*cos(8*pi*t);k=0:N/2;w=2*pi/N*k;X=fft(x,N);magX=abs(X(1:N/2+1);subplot(2,1,1);stem(n,x,.);title(signal x(n);grid on;subplot(2,1,2);plot(w/pi,magX);title(FFT N=40);xlabel(f (unit:pi);ylabel(|X|);grid on;结果截图：能观察出信号的2个频率分量 MATLAB程序：N=128;n=0:N-1;t=0.01*n;x=2*cos(4*pi*t)+5*cos(8*pi*t);k=0:N/2;w=2*pi/N*k;X=fft(x,N);magX=abs(X(1:N/2+1);subplot(2,1,1);stem(n,x,.);title(signal x(n);grid on;subplot(2,1,2);plot(w/pi,magX);title(FFT N=128);xlabel(f (unit:pi);ylabel(|X|);grid on;结果截图：幅度频谱变得更加密集，模拟频率和数字频率各为4hz和100hz 频谱分析结果与理论相一致的。（二）一个连续信号含3个频率分量，经采样得到以下序列：xn=sin2*0.15n+cos20.15+fn+cos20.15+2fn,n=0,N-1。已知N=16，f分别为1/16,1/64，观察其频谱；当N=64、128，f不变，其结果有何不同，为什么？分析参数变化对频率分辨率的影响。解：N=16，f=1/16,1/64：MATLAB程序：N=16;n=0:N-1;df1=1/16;df2=1/64;k=0:N/2;w=2*pi/N*k;x1=sin(2*pi*0.15*n)+cos(2*pi*(0.15+df1)*n)+cos(2*pi*(0.15+2*df1)*n);x2=sin(2*pi*0.15*n)+cos(2*pi*(0.15+df2)*n)+cos(2*pi*(0.15+2*df2)*n);X1=fft(x1,N);magX1=abs(X1(1:N/2+1);X2=fft(x2,N);magX2=abs(X2(1:N/2+1);subplot(2,2,1);plot(n,x1);title(N=16 df=1/16);grid on;subplot(2,2,2);plot(n,x2);title(N=16 df=1/64);grid on;subplot(2,2,3);plot(w/pi,magX1);title(FFT N=16 df=1/16);xlabel(f (unit:pi);ylabel(|X|);grid on;subplot(2,2,4);plot(w/pi,magX2);title(FFT N=16 df=1/64);xlabel(f (unit:pi);ylabel(|X|);grid on;结果截图：N=64和128，f=1/64：MATLAB程序：N1=64;n1=0:N1-1;N2=128;n2=0:N2-1;df=1/64;x1=sin(2*pi*0.15*n1)+cos(2*pi*(0.15+df)*n1)+cos(2*pi*(0.15+2*df)*n1);x2=sin(2*pi*0.15*n2)+cos(2*pi*(0.15+df)*n2)+cos(2*pi*(0.15+2*df)*n2);k1=0:N1/2;w1=2*pi/N1*k1;k2=0:N2/2;w2=2*pi/N2*k2;X1=fft(x1,N1);magX1=abs(X1(1:N1/2+1);X2=fft(x2,N2);magX2=abs(X2(1:N2/2+1);subplot(2,2,1);plot(n1,x1);title(N=64 df=1/64);grid on;subplot(2,2,2);plot(n2,x2);title( N=128 df=1/64);grid on;subplot(2,2,3);plot(w1/pi,magX1);title(FFT N=64 df=1/64);xlabel(f (unit:pi);ylabel(|X|);grid on;subplot(2,2,4);plot(w2/pi,magX2);title(FFT N=128 df=1/64);xlabel(f (unit:pi);ylabel(|X|);grid on;结果截图：由上图可知，当f越小时，频谱分辨率越低。当N越大时，频谱分辨率越高。N从16变化为64,128时，谱峰逐渐清晰，截取长度增加，谱线变得密集，频谱更加接近真实值，泄露和混叠效应、栅栏效应都变小。（三）被噪声污染的信号，比较难看出所包含的频率分量，如一个由50Hz和120Hz正弦信号构成的信号，受零均值随机噪声的干扰，数据采样率为1000Hz。选取合适的采样点数，使用FFT函数来分析其信号频率成分，并绘制信号的时域和频域波形。解：MATLAB程序：t=0:0.001:0.8;x=sin(2*pi*50*t)+cos(2*pi*120*t);y=x+1.5*randn(1,length(t);subplot(3,1,1);plot(t,x);subplot(3,1,2);plot(t,y);Y=fft(y,512);P=Y.*conj(Y)/512;f=1000*(0:255)/512;subplot(3,1,3);plot(f,P(1:256);结果截图：（四）已知AM调幅信号xt=A0+mtcosw0t，其中mt=Acost，以采样频率fs，采样长度N，对信号x(t)进行采样，频率分辨率为1Hz。求采样信号的频谱和功率谱，验证不发生频谱泄露的条件。解：频谱泄露条件f0=m*fs/N=m*F0，此题f0易证为(w/2pi+w0/2pi)，代入上式可得条件为f0是整数时不会产生频谱泄露 当f0为整数时：MATLAB程序：F0=1;N=40;A0=5;A=4;w0=2*pi*2;w=2*pi*8;fs=N*F0;n=0:(N-1);t=1/fs*n;x=(A0+A*cos(w*t).*cos(w0*t);k=0:N/2;w=2*pi/N*k;X=fft(x,N);magX=abs(X(1:N/2+1);subplot(1,2,1);plot(w/pi,magX);title(满足条件的频谱图);xlabel(f (unit:pi);ylabel(|x|);grid ;window=boxcar(length(x);nfft=1024;Pxx,f=periodogram(x,window,nfft,fs);subplot(1,2,2);plot(f,10*log10(Pxx);title(满足条件的功率谱图);grid on;结果截图：当f0不为整数时：MATLAB程序：F0=1;N=40;A0=5;A=4;w0=2*pi*2;w=2*pi*8.4;fs=N*F0;n=0:(N-1);t=1/fs*n;x=(A0+A*cos(w*t).*cos(w0*t);k=0:N/2;w=2*pi/N*k;X=fft(x,N);magX=abs(X(1:N/2+1);subplot(1,2,1);plot(w/pi,magX);title(不满足条件的频谱图);xlabel(f (unit:pi);ylabel(|x|);grid ;window=boxcar(length(x);nfft=1024;Pxx,f=periodogram(x,window,nfft,fs);subplot(1,2,2);plot(f,10*log10(Pxx);title(不满足条件的功率谱图);grid on;结果截图：（五）研究高密度频谱和高分辨率频谱1.高密度频谱：（1）定义：当信号的时间域长度不变时，在频域内对它的频谱进行提高采样频率，得到高密度频谱。（2）特点：由于信号末尾补零没有对源信号增加任何新的信息，因此不能提高频率分辨率，但可以减少栅栏效应，细化当前分辨率下的频谱。（3）应用：减少栅栏效应，细化当前分辨率下的频谱 2.高分辨率频谱：（1）定义：在维持采样频率fs不变的情况下，为提高分辨率智能增加采样点数N，此时所得到的的频谱称为高分辨率频谱。（2）特点：分辨率高（3）应用：提高频谱分辨率（六）考虑一连续信号xat=cos26.5103t+cos27103t+cos(29103t) 以采样率fs=32kHz对信号进行xa(t)采样，分析如下情况的幅频特性： 采集数据长度N=16，做16点的FFT；采集长度N=16点，补零至256点，做256点FFT。 采集数据长度N=64，做64点的FFT；采集数据长度N=16，补零至256点，做256点的FFT。 采集长度N=256点，做256点的FFT。观察以上3种设置的幅频特性曲线，分析和比较它们的特点，并说明原因。解：MATLAB程序：N=16;n=0:N-1;fs=32000;xa=cos(2*pi*6.5*103*n/fs)+cos(2*pi*7*103*n/fs)+cos(2*pi*9*103*n/fs);Xa=fft(xa,N);f=fs*(0:15)/16;magX=abs(Xa);subplot(2,2,1);stem(n,xa,.);title(N=16 16点);grid on;subplot(2,2,3);plot(f,magX);xlabel(Hz);ylabel(幅度);title(FFT N=16 16点);grid on;x=xa(1:16),zeros(1,240);X=fft(x,256);f=fs*(0:255)/256;magX=abs(X);subplot(2,2,2);stem(n,xa,.);title(N=16 256点);grid on;subplot(2,2,4);plot(f,magX);xlabel(Hz);ylabel(幅度);title(FFT N=16 256点);grid on;结果截图：MATLAB程序：N=64;n=0:N-1;fs=32000;xa=cos(2*pi*6.5*103*n/fs)+cos(2*pi*7*103*n/fs)+cos(2*pi*9*103*n/fs);Xa=fft(xa,N);f=fs*(0:63)/64;magX=abs(Xa);subplot(2,2,1);stem(n,xa,.);title(N=64 64点);grid on;subplot(2,2,3);plot(f,magX); grid on;xlabel(Hz);ylabel(幅度);title(FFT N=64 64点);x=xa(1:16),zeros(1,240);X=fft(x,256);f=fs*(0:255)/256;magX=abs(X);subplot(2,2,2);stem(n,xa,.);title(N=16 256点);grid on;subplot(2,2,4);plot(f,magX); grid on;xlabel(Hz);ylabel(幅度);title(FFT N=16 256点);结果截图：MATLAB程序：N=256;n