1、1实验名称：数字信号处理实验 实验项目：用 FFT 作谱分析 指导老师： 李秋菊 班级 ： 10 电科 姓名： 张卫娟 学号：201000804084 成绩： 一、实验目的：1、在理论学习的基础上，通过本实验，加深对 FFT 的理解，熟悉 MATLAB 中的有关函数。2、熟悉应用 FFT 对典型信号进行频谱分析的方法。熟悉 FFT 算法原理和 FFT 子程序的应用。3、学习用 FFT 对连续信号和时域离散信号进行谱分析的方法。了解应用 FFT 进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用 FFT。2、实验原理：（一）、在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序
2、列，我们可以使用离散傅里叶变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列 x(n)的长度为 N 时，它的 DFT 定义为：反变换为：有限长序列的 DFT 是其 Z 变换在单位圆上的等距采样，或者说是序列 Fourier变换的等距采样，因此可以用于序列的谱分析。在信号处理中，DFT 的计算具有举足轻重的地位， ，信号的相关、滤波、谱估计等都要通过 DFT 来实现。然而，当 很大的时候，求一个 点的 DFTNN要完成 次复数乘法和 次复数加法，其计算量相当大。1965 年N)1(J.W.Cooley 和 J.W.Tukey 巧妙地利用 因子的周期
3、性和对称性，构造了一个NWDFT 快速算法，即快速傅立叶变换(FFT)。FFT 并不是与 DFT 不同的另一种变换，而是为了减少 DFT 运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减2少运算量。常用的 FFT 是以 2 为基数的，其长度 。它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于 2 的整数次方时，为了使用以 2 为基数的 FFT，可以用末位补零的方法，使其长度延长至 2 的整数次方。（二）、在运用 DFT 进行频谱分析的过程中可能的产生混叠误差序列的频谱是被采样信号频谱的周期延拓，当采样速率不满足 Nyquist 定理时，就会发生频
4、谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。（三）、matlab 函数应用： MATLAB 为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有 Fft、Ifft、Fft2 、Ifft2, Fftn、ifftn 和 Fftshift、Ifftshift 等。当所处理的数据的长度为 2 的幂次时，采用基-2 算法进行计算，计算速度会显著增加。所以，要尽可能使所
5、要处理的数据长度为 2 的幂次或者用添零的方式来添补数据使之成为 2 的幂次。1、fft 和 Ifft 函数（1）调用方式：Y=fft(X)参数说明如果 X 是向量，则采用傅立叶变换来求解 X 的离散傅立叶变换；如果 X 是矩阵，则计算该矩阵每一列的离散傅立叶变换；如果 X 是（N 维数组，则是对第一个非单元素的维进行离散傅立叶变换；)D（2） Yfft(X,N)3参数说明N 是进行离散傅立叶变换的 X 的数据长度，可以通过对 X 进行补零或截取来实现。（3）Yfft(X,dim) 或 Yfft(X,N,dim)参数说明在参数 dim 指定的维上进行离散傅立叶变换；当 X 为矩阵时，dim 用
6、来指定变换的实施方向： dim=1,表明变换按列进行；dim=2 表明变换按行进行。函数 Ifft 的参数应用与函数 Fft 完全相同。应用说明【实例 1】fft 的应用X=2 1 2 8;Y=fft(X，4)运行结果Y13.0000 0+7.0000i -5.0000 0-7.0000【实例 2】fft(X,N,dim)的应用A=2 5 7 8;1 4 0 5;3 8 5 1;9 1 2 7;Z=fft(A,1)42、 Fftshift 和 Ifftshift 函数调用方式Z=fftshift(Y)此函数可用于将傅立叶变换结果 Y（频域数据）中的直流成分（即频率为0 处得值）移到频谱的中间位
7、置。【实例 3】 fftshift 的应用X=rand(5,4);y=fft(X);z=fftshift(y);%只将傅立叶变换结果 y 中的直流成分移到频谱的中间位置.运行结果：y=3.2250 2.5277 1.4820 1.63140.3294+0.2368i 0.0768+0.3092i 0.6453+0.4519i -0.7240-0.4116i-0.2867-0.6435i 0.5657+0.4661i -0.5515+0.2297i -0.0573-0.0881i-0.2867+0.6435i 0.5657-0.4661i -0.5515-0.2297i -0.0573+0.08
8、81i0.3294-0.2368i 0.0768-0.3092i 0.6453-0.4519i -0.7240+0.4116iZ=-0.5515-0.2297i -0.0573+0.0881i -0.2867+0.6435i 0.5657-0.4661i0.6453-0.4519i -0.7240+0.4116i 0.3294-0.2368i 0.0768-0.3092i51.4820 1.6314 3.2250 2.52770.6453+0.4519i -0.7240-0.4116i 0.3294+0.2368i 0.0768+0.3092i-0.5515+0.2297i -0.0573-0
9、.0881i -0.2867-0.6435i 0.5657+0.4661i【实例 4】fft 在信号分析中的应用使用频率分析方法对模拟信号 x(t) x=sin(2*pi*100*t)进行频谱分析。采样频率fs=1000，采样点数 N=512。并画出信号的时域波形及 FFT 变换后的幅频响应与相频响应。程序：fs=512;%采样频率N=512;%数据点数n=0:N-1;t=0:1/fs:(N-1)/fs;%采样时间序列f0=100;%信号频率x=sin(2*pi*f0*t);subplot(3,1,1);plot(t,x);xlabel(t);ylabel(sin(2*pi*100*t);ti
10、tle(时域信号 );Y=fft(x,N);%对信号进行 FFT 变换magY=abs(Y);%求得 FFT 变换后的幅度angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);plot(f,magY);%画出幅频响应xlabel(f);ylabel(幅度);title(N=512);grid;subplot(3,1,3);plot(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=512);grid;60 0.1 0.2 0.3 0.4 0.5 0.6 0.7
11、-101tsin(2*pi*100*t)位 位 位 位0 100 200 300 400 500 600 700 800 900 10000200400f位位N=5120 100 200 300 400 500 600 700 800 900 1000-2000200f位位N=512结果分析：假设采样频率为 Fs，信号频率 F，采样点数为 N。那么 FFT 之后结果就是一个为 N 点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。假设原始信号的峰值为 A，那么 FFT 的结果的每个点（除了第一个点直流分量之外）的模值就是 A 的 N/2 倍。（例如 f=100Hz
12、点处，A=1*N/2=1*512/2=256,而第一个点就是直流分量，它的模值就是直流分量的 N 倍)。magY=abs(Y);%求得 FFT 变换后的幅度如果要换算为实际的幅度需要改为：magY=magY/(N/2);magY(1)=magY(1)/2;每个点的相位,就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz） ，而最后一个点 N 的再下一个点（假设的第 N+1 个点）则表示采样频率Fs，这中间被 N-1 个点平均分成 N 等份，每个点的频率依次增加。 第 n 个点所表示的频率为：Fn=(n-1)*Fs/N 。由上面的公式可以看出， Fn 所能分辨到频率为为 Fs/N，如果采
13、样频率 Fs 为 1024Hz，采样点数为 1024 点，则可以分辨到1Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。由于 FFT 结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。前面的程序由：plot(f(1:N/2),magY(1:N/2);%画出幅频响应三、实验内容：1、 对被白噪声污染的信号 进行频2/150cos.6/10cos32Stt谱分析，从中鉴别出有用的信号。7要求：将信号的幅度换算成实际的幅度，信号的频率换算成实际的频率程序： clear;fs=512;%采样频率N=512;%数据点数n=0:N-1;t=0:1/fs:(N-1)/fs;
14、%采样时间序列f0=100;%信号频率x=2+3*cos(100*pi*t-pi/6)+1.5*cos(150*pi*t+pi/2)+randn(1,N);subplot(3,1,1);plot(t,x);xlabel(t);ylabel(sin(2*pi*100*t);title(时域信号 );Y=fft(x,N);%对信号进行 FFT 变换magY=abs(Y)/N*2;%求得 FFT 变换后的幅度magY(1)=magY(1)/2;angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);plot(f(1:N/2),m
15、agY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title(N=512);grid;subplot(3,1,3);plot(f(1:N/2),angY(1:N/2);%画出相频响应 xlabel(f);ylabel(相位)title(N=512);grid;图形：2、 对连续的单一频率周期信号, ,信号频率 f=1Hzt2sinS按采样频率 fs=8Hz 采样，截取长度 N 分别选 N =20 和 N =16，观察其 DFT 结果的幅度谱。幅频N=20程序：clear;fs=8;%采样频率N=20;%数据点数n=0:N-1;t=0:1/fs:(N-1)/fs;%采
16、样时间序列f0=100;%信号频率x=sin(2*pi*t);subplot(3,1,1);stem(t,x);xlabel(t);ylabel(sin(2*pi*100*t);title(时域信号 );Y=fft(x,N);%对信号进行 FFT 变换magY=abs(Y)/(N/2);%求得 FFT 变换后的幅度magY(1)=magY(1)/2;8angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title
17、(N=20);grid;subplot(3,1,3);stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=20);grid;N=16程序：clear;fs=16;%采样频率N=16;%数据点数n=0:N-1;x=1,2,3,4,4,3,2,1,zeros(1,8);subplot(3,1,1);stem(x);xlabel(n);ylabel(sin(2*pi*100*t);title(时域信号 );Y=fft(x);%对信号进行 FFT 变换magY=abs(Y)/N*2;%求得 FFT 变换后的幅度magY(1)=mag
18、Y(1)/2;angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title(N=16);grid;subplot(3,1,3);stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=16);grid;图形：3、 对三角波序列 ， 进行频谱分析74,n-8301)(xc elsnxc,0)(设 n=16。程序：clear;fs=16;%采样频率N
19、=16;%数据点数n=0:N-1;x=1,2,3,4,4,3,2,1,zeros(1,8);9subplot(3,1,1);stem(x);xlabel(n);ylabel(sin(2*pi*100*t);title(时域信号 );Y=fft(x);%对信号进行 FFT 变换magY=abs(Y)/N*2;%求得 FFT 变换后的幅度magY(1)=magY(1)/2;angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(
20、幅度);title(N=16);grid;subplot(3,1,3);stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f);ylabel(相位)title(N=16);grid;图形：4、 对矩形波进行频谱分析。矩形波为脉冲宽度为 2s，持续时间为-55s，采用矩形脉冲函数 rectpuls(t,2);程序：clear;fs=64;%采样频率N=64;%数据点数n=0:N-1;t=-5:1/fs:5;x=rectpuls(t,2);subplot(3,1,1);plot(t,x);xlabel(n);ylabel(sin(2*pi*100*t);title(时
21、域信号 );Y=fft(x);%对信号进行 FFT 变换magY=abs(Y)/(N/2);%求得 FFT 变换后的幅度magY(1)=magY(1)/2;angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);plot(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f);ylabel(幅度);title(N=64);grid;10图形：四 实验结论通过此次试验，加深对 FFT 的理解，在各种信号序列中，对有限长序列，我们可以使用离散傅里叶变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，当序列 x(n)的长度为 N 时，它的 DFT 定义为：，反变换为：在信号处理中，信号的相关、滤波、谱估计等都要通过 DFT 来实现。但是，当很大的时候，计算很困难，我们可以巧妙地利用 因子的周期性和对称性，N NW构造了一个 DFT 快速算法，即快速傅立叶变换(FFT)来计算。当然，在进行频谱分析的过程中，当采样速率不满足 Nyquist 定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。为避免混叠现象，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。
