FFT算法的应用--频率估计.docx

11信工2班 黄锌 07 实验2 FFT算法的应用---频率估计 实验目的：加深对离散信号的DFT的理解及其FFT算法的运用。 实验原理：N点序列的DFT和IDFT变换定义式如下： ， 利用旋转因子具有周期性，可以得到快速算法（FFT）。 在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列的DFT正、反变换。 实验内容： 1． 对连续的单一频率周期信号 按采样频率 采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。 程序如下： n1=[0:1:15]; x1=sin(2*pi*n1/8); figure(1) grid; subplot(2,1,1); plot(n1,x1); x11=fft(x1); x11=abs(x11); subplot(2,1,2); stem(n1,x11); n2=[0:1:19]; x2=sin(2*pi*n2/8); figure(2) subplot(2,1,1); plot(n2,x2); x21=fft(x2); x21=abs(x21); subplot(2,1,2); stem(n2,x21); 运行结果如下： N=16 N=20 结果分析：所用的信号是一个周期为8的正弦信号，分别截取长度16和长度20来进行DFT变换。前者得到的频率幅度谱只有两个???率分量，后者有多个频率分量，原因是16是周期的整数倍，20不是周期的整数倍。 2 . 2N点实数序列 N=64。用一个64点的复数FFT程序，一次算出，并绘出。 程序如下： %方法一（复数FFT程序） N=64; n=[0:1:N-1]; n1=2*n; n2=2*n+1; k=[0:1:N-1]; xn1=cos(2*pi/N*7*n1)+1/2*cos(2*pi/N*19*n1); %偶序列 xn2=cos(2*pi/N*7*n2)+1/2*cos(2*pi/N*19*n2); %奇序列 XK1=fft(xn1); XK2=fft(xn2); X1=XK1+exp(-j*pi*k/N).*XK2; X2=XK1-exp(-j*pi*k/N).*XK2; X1=[X1 zeros(1,N)]; X2=[zeros(1,N) X2]; XK=X1+X2; k=[0:1:2*N-1]; XK=abs(XK); figure(1) stem(k,XK); xlabel(k);ylabel(|X(k)|); title(X(k)=DFT[x(n)]2N) %方法2(直接FFT) N1=64; n0=0:2*N1-1; x1=cos(2*pi*7*n0/N1)+1/2*cos(2*pi*19*n0/N1); Xk1=fft(x1,128); k1=n0; figure(2) stem(k1,abs(Xk1)) xlabel(k1);ylabel(|X(k1)|); title(直接FFT X(k1)) 运行结果： 结果分析：第二种方法是直接对信号进行FFT变换，实际上是对第一种方法的一种检测，二者运行结果完全一致，说明结果正确。 3. 频率估计 1）产生一个单频实信号，加上一定信噪比的噪声。 2）对含噪声的信号进行频率估计。估计方法可以查阅相关文献。 3）统计估计出来的频率和真实频率之间的误差。 4）验证该频率估计算法在不同信噪比、不同数据长度下、不同频率时候的性能。 程序如下： Function(a) x=(0:1000); y1=sin(200*pi*x/1000);%原始频率100Hz figure(1) plot(x,y1); axis([0,50,-1,1]); title(未受干扰时信号的波形); ylabel(y); xlabel(x/1000); grid; y2=awgn(y1,a); y=y1+y2; figure(2); plot(x,y); axis([0,200,-3,3]); title(叠加了高斯白噪声的信号波形); ylabel(y); xlabel(x/1000); grid; Y=fft(y); Y=abs(Y); f=x; figure(3); plot(f,Y) axis([0,200,-100,2000]); title(含噪声信号y(t)的频谱) xlabel(频率(Hz)) ylabel(|Y(f)|) 运行结果： （1）信噪比为1时 （2）信噪比为5时 （3）信噪比为10时 （4）信噪比为20时