MATLAB实现FFT及信号的谱分析

一、实验目的 1.通过实验加深对 FFT 理解，熟悉 FFT 程序、结构和编程方法。

２.熟练应用 FFT 分析典型信号的方法。

３.了解应用 FFT 在实践中正确应用信号频域分析可能出现的问题FFT。

4. 理解 FFT 与 IFFT 的关系。

5..熟悉应用 FFT 实现两个序列的线性卷积的方法。

二、实验原理及方法 在各种信号序列中，有限长序列信号处理占有很重要的位置，对有限长序列，我们可以使用离散 Fourier变换（DFT）。这种变化不仅能很好地反映序列的频谱特性，而且在计算机上实现了永快算法。 x(n) 的长度为 N 时，它的 DFT 定义为：

长序列有限 DFT 是其 Z 在单位圆上变换等距采样或序列 Fourier 等距采样可用于序列分析。FFT 并不是与 DFT 另一种不同的变化是为了减少 DFT操作次数的快速算法。它是对变换式进行一次次的分解，是其成为若干小点数的组合，从而减少运算量。常用的 FFT 是以 2 其长度为基数 N = 2L。它的效率高，程序简单，使用非常方便。 2 为了使用以2为基数的整数次方 FFT，末位补零的方法可以延长到 2 整数次方。

(一)在使用 DFT 三个误差可能发生在频谱分析中。 (1) 混叠 为了计算连续信号的频谱，首先需要对连续信号进行取样。如果取样频率过低，即取样周期过大，则频域内会出现混合现象，不可能无失真恢复原始连续信号 数字。当处理模拟信号的最高频率为带限信号时 fh 与抽样频率 fs 满足fs ≥ 2 fh频谱混合不会出现。 然而，时域内有限长信号的频谱宽度是无限的。为了使有限长信号满足抽样定理，可以在抽样前用低通模拟滤波器过滤信号，以确保高于折叠频率的重量不会出现。 （2）泄漏 实际的信号序列通常是长的，甚至是无限的。为了方便起见，我们经常使用短的序列来接近它们。这可以使用较短的序列 DFT 对信号进行频谱分析。 x(n) 截短的过程就是将 原始信号序列乘以矩形窗口函数的过程是频域中两个频谱的卷积。一般来说，由此获得的频谱与原始信号频谱不同，称为泄漏。在实际应用中，可以选择频率 谱主瓣小，旁瓣小，尽量接近(Ω) 窗函数减少泄漏。 泄漏不能与混合物完全分离，因为泄漏导致频谱扩展，从而导致混合物。为了减少泄漏的影响，可以选择合适的窗口函数尽量减少频谱的扩散。 (3)栅栏效应 DFT 是对单位上 z 变换均匀取样，不能将频谱视为连续函数。这就产生了栅栏效应。从某种意义上说，用 DFT看频谱就像通过尖桩围栏看图片，只能在离散点看到真正的频谱。这样，一些频谱的峰点或谷点就可能被尖桩围栏挡住，我们无法观察到。减小栅栏效应的一个方法是借助于原序列的末端增加一些零值，从而变动 DFT 的点数。事实上，这种方法人为地改变了真实频谱采样的点数和位置，相当于移动每个尖桩围栏的位置，从而可以看到原始频谱的峰点或谷点。 IFFT 一般可以通过 FFT 只要是正确的，程序就完成了 X[k]取共轭，进行 FFT 运算，然后取共轭，乘以因子 1/N，就可以完成 IFFT。 实验中使用的信号序列：

三、实验内容 1.用三种不同的 DFT 程序计算 x(n) = R? (n) 傅里叶变换 运行时间。 （1）用 for loop 语句的M 函数文件dft1.m，循环变量逐点计算 X (k ) ； （2）编写用 MATLAB 矩阵运算的 M 函数文件 dft2.m, 完成以下矩阵操作； （3）调用 FFT 直接计算库函数 X (k ) ； (4)以上三种不同的方式分别编写 DFT 程序计算序列 x(n) 傅立叶变换 X (e ?? ) ，并绘制相应的幅频和相频特性，然后比较每个程序的计算机运行时间。 M 函数文件如下：

dft1.m: function[Am,pha]=dft1(x) N=length(x); w=exp(-j*2*pi/N); for k=1:N sum=0; for n=1:N sum=sum x(n)*w^((k-1)*(n-1)); end Am(k)=abs(sum); pha(k)=angle(sum); end dft2.m: function[Am,pha]=dft2(x) N=length(x); n=[0:N-1]; k=[0:N-1]; w=exp(-j*2*pi/N); nk=n'*k; wnk=w.^(nk); Xk=x*wnk; Am=abs(Xk); pha=angle(Xk);

四、思考题 FFT 周期信号序列的频谱也可用于分析？ 假如正弦信号sin(2Πf?k ), f? = 0.1Hz ，用 16 点来做 FFT 获得的频谱是信号本身的真实谱吗？为什么？

五、实验报告要求 1.简要介绍实验原理和目的。 2.给出准备好的实验主程序、实验信号序列的时域和频域图形，并分析获得的图形，说明参数 变化时对时域和频域信号波形的影响。 3.简单回答思考题。

x=ones(1,8); n=0:length(x)-1; t0=clock; [m1,a1]=dft1(x); subplot(321),stem(n,m1,'.'); title('用dft1实现 mag'); subplot(322),stem(n,a1,'.'); title('angle'); t1=etime(clock,t0); t0=clock; [m2,a2]=dft2(x); subplot(323),stem(n,m2,'.'); title('用dft2实现 mag'); subplot(324),stem(n,a2,'.'); title('angle'); t2=etime(clock,t0); t0=clock; xk=fft(x); subplot(325),stem(n,abs(xk),'.'); title('用FFT实现 mag'); subplot(326),stem(n,angle(xk),'.'); title('angle'); t3=etime(clock,t0); figure; subplot(311),stem(t1,'.');ylabel('单位：s');title('用dft1实现'); subplot(312),stem(t2,'.');ylabel('单位：s');title('用dft2实现'); subplot(313),stem(t3,'.');ylabel('单位：s');title('用FFT实现'); 调用函数： ?t1.m function[Am,pha]=dft1(x) N=length(x); w=exp(-j*2*pi/N); for k=1:N sum=0; for n=1:N sum=sum x(n)*w^((k-1)*(n-1)); end Am(k)=abs(sum); pha(k)=angle(sum); end ?t2.m function[a,p]=dft2(x) N=length(x); n=0:N-1; k=n; w=exp(-j*2*pi/N); nk=n'*k; wnk=w.^(nk); xk=x*wnk; a=abs(xk); p=angle(xk);

p=8; n=0:15; figure; for k=1:3 q=2^k; x1=exp(-((n-p).^2)/q); subplot(3,2,2*k-1),stem(n,x1,'.'); title('q变化 x(n)'); [m,a,w]=dtft(x1); subplot(3,2,2*k),plot(w/pi,m); title('q变化 abs[X(w)]'); end q=8; figure; for k=1:3 p=-2*(k^2) 11*k-1; x1=exp(-((n-p).^2)/q); subplot(3,2,2*k-1),stem(n,x1,'.'); title('p变化 x(n)'); [m,a,w]=dtft(x1); subplot(3,2,2*k),plot(w/pi,m); title('p变化 abs[X(w)]'); end figure; f1=[0.0625,0.4375,0.5625]; for k=1:3 f=f1(k); a=0.1; x2=exp(-a)*sin(2*pi*f*n); subplot(3,2,2*k-1),stem(n,x2,'.'); title('f 变化 x(n)'); [m,a,w]=dtft(x2); subplot(3,2,2*k),plot(w/pi,m); title('f 变化 abs[X(w)]'); end 调用函数： %dtft.m function[m,a,w]=dtft(x) N=length(x); n=0:N-1; w=linspace(-2*pi,2*pi,500); y=x*exp(-j*n'*w); m=abs(y); a=angle(y);

nbsp;

3.程序：
f=64;
dt=1/f;
N=[16,32,64];
for k=1:3
n=0:N(k)-1;
nt=n*dt;
x=cos(8*pi*nt)+cos(16*pi*nt)+cos(20*pi*nt);
y=fft(x);
subplot(3,3,3*k-2),stem(n,x,'.');
title('x(n)');
subplot(3,3,3*k-1),stem(n,abs(y),'.');
title('abs(X(k))');
subplot(3,3,3*k),stem(n,angle(y),'.');
title('angle(X(k))');
end

n=0:127;
x=cos(n*pi/36)+cos(1.5*n*pi/36);
y=fft(x);
figure;
subplot(311),stem(n,x,'.');
title('x(n)');
subplot(312),stem(n,abs(y),'.');
title('X(k) 幅度谱');
subplot(313),stem(n,angle(y),'.');
title('X(k) 相位谱');
N=length(n);
n1=0:N/4-1;
x1=cos(4*n1*pi/36)+cos(1.5*4*n1*pi/36);
x1=[x1,zeros(1,N*3/4)];
y1=fft(x1,128);
figure;
subplot(311),stem(n,x1,'.');
title('x(4n)');
subplot(312),stem(n,abs(y1),'.');
title('对应的X(k) 幅度谱');
subplot(313),stem(n,angle(y1),'.');
title('对应的X(k) 相位谱');
n2=0:(N-1)*4;
x2=zeros(1,N*4-3);
for k=0:127
x2(k*4+1)=cos(k*pi/36)+cos(1.5*k*pi/36);
end
x3=zeros(1,N);
x3=x2(1:N);
y2=fft(x3,128);
figure;
subplot(311),stem(n,x3,'.');
title('x(n/4)');
subplot(312),stem(n,abs(y2),'.');
title('对应的X(k) 幅度谱');
subplot(313),stem(n,angle(y2),'.');
title('对应的X(k) 相位谱');