数字信号处理实验——用FFT对连续信号和时域离散信号进行谱分析

一、实验目的：

1、在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。

2、熟悉应用FFT对典型信号进行频谱分析的方法。熟悉FFT算法原理和FFT子程序的应用。

3、学习用FFT对连续信号和时域离散信号进行谱分析的方法。了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。

二、实验原理：

（一）、在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散傅里叶变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT定义为：

反变换为：

有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。

在信号处理中，DFT的计算具有举足轻重的地位，，信号的相关、滤波、谱估计等都要通过DFT来实现。然而，当N很大的时候，求一个N点的DFT要完成NN次复数乘法和N(N1)次复数加法，其计算量相当大。1965年J.W.Cooley和J.W.Tukey巧妙地利用WN因子的周期性和对称性，构造了一个DFT快速算法，即快速傅立叶变换(FFT)。

（二）、在运用DFT进行频谱分析的过程中可能的产生混叠误差

序列的频谱是被采样信号频谱的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号

1

进行滤波。

（三）、matlab函数应用：

MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，

主要有Fft、Ifft、Fft2 、Ifft2, Fftn、ifftn和Fftshift、Ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。

1、 fft和Ifft函数

（1）调用方式：Y=fft(X) 参数说明

·如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换； ·如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换；

·如果X是（ND)维数组，则是对第一个非单元素的维进行离散傅立叶变换； （2） Y＝fft(X,N) 参数说明

N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。 （3）Y＝fft(X,[],dim) 或Y＝fft(X,N,dim)

参数说明

·在参数dim指定的维上进行离散傅立叶变换；

·当X为矩阵时，dim用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2表明变换按行进行。

函数Ifft的参数应用与函数Fft完全相同。 应用说明

【实例1】fft的应用

X=[2 1 2 8]; Y=fft(X)

 2

运行结果

Y＝ 13.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)

2、 Fftshift和Ifftshift函数

调用方式

Z=fftshift(Y)

此函数可用于将傅立叶变换结果Y（频域数据）中的直流成分（即频率为0处得值）移到频谱的中间位置。

【实例3】 fftshift的应用 X=rand(5,4); y=fft(X);

z=fftshift(y);%只将傅立叶变换结果y中的直流成分移到频谱的中间位置. 运行结果：

y= 3.2250 2.5277 1.4820 1.6314 0.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.0881i 0.3294-0.2368i 0.0768-0.3092i 0.6453-0.4519i -0.7240+0.4116i Z= -0.5515-0.2297i -0.0573+0.0881i -0.2867+0.6435i 0.5657-0.4661i 0.6453-0.4519i -0.7240+0.4116i 0.3294-0.2368i 0.0768-0.3092i 1.4820 1.6314 3.2250 2.5277

3

0.6453+0.4519i -0.7240-0.4116i 0.3294+0.2368i 0.0768+0.3092i -0.5515+0.2297i -0.0573-0.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=1000;%采样频率 N=512;%数据点数 n=0:N-1;

t=0:1/fs:N/fs;%采样时间序列 f0=100;%信号频率 x=sin(2*pi*f0*t); subplot(3,1,1); plot(t,x); xlabel('t');

ylabel('sin(2*pi*100*t)'); title('时域信号');

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,angY);%画出相频响应 xlabel('f'); ylabel('相位') title('N=512'); grid;

4

时域信号sin(2*pi*100*t)10-100.10.20.3tN=5120.40.50.60.7400幅度20000100200300400500fN=5126007008009001000200相位0-2000100200300400500f6007008009001000

结果分析：假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。

这个点的模值，就是该频率值下的幅度特性。假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。

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，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。

前面的程序由：plot(f,magY);%画出幅频响应 三|、实验内容：

1、 对被白噪声污染的信号S23cos100t/61.5cos150t/2进行频谱分析，从中鉴别出有用的信号。

要求：将信号的幅度换算成实际的幅度，信号的频率换算成实际的频率

fs=1000;%采样频率 N=512;%数据点数 n=0:N-1; randn(1,N);

t=0:1/fs:(N-1)/fs;%采样时间序列 f0=100;%信号频率

x=2+3*cos(100*pi*t-pi/6)+(3/2)*cos(150*pi*t+pi/2); subplot(3,1,1); plot(t,x); xlabel('t');

ylabel('sin(2+3*cos(100*pi*t-pi/6)+(3/2)*cos(150*pi*t+pi/2))'); title('时域信号');

5

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,angY);%画出相频响应 xlabel('f'); ylabel('相位') title('N=512'); grid;

2、 对连续的单一频率周期信号, Ssin2t ,信号频率f=1Hz

按采样频率fs=8Hz采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。

并从中得出什么结论？

程序：fs=8;%采样频率 N=16;%数据点数 n=0:N-1;

t=0:1/fs:N/fs;%采样时间序列 f0=1;%信号频率 x=sin(2*pi*t); subplot(3,1,1);

6

stem(t,x); xlabel('t');

ylabel('sin(2*pi*t)'); title('时域信号');

Y=fft(x,N);%对信号进行FFT变换

magY=abs(Y);%求得FFT变换后的幅度

angY=angle(Y)*180/pi;%求得FFT变换后的相位 f=n*fs/N;%频率序列 subplot(3,1,2);

stem(f,magY);%画出幅频响应 xlabel('f'); ylabel('幅度'); title('N=512'); grid;

subplot(3,1,3);

stem(f,angY);%画出相频响应 xlabel('f'); ylabel('相位') title('N=512'); grid;

结果：

n1,0n33、 对三角波序列xc(n)，xc(n)0,else进行频谱分析设n=16。

8-n,4n7程序：clear;

7

fs=16;%采样频率 N=16;%数据点数 n=0:N-1;

t=0:1/fs:(N-1)/fs;%采样时间序列 f0=100;%信号频率

x=[1,2,3,4,4,3,2,1,zeros(1,8)]; subplot(3,1,1); stem(t,x); xlabel('t');

ylabel('[1,2,3,4,4,3,2,1,zeros(1,8)]'); title('时域信号');

Y=fft(x,N);%对信号进行FFT变换

magY=abs(Y);%求得FFT变换后的幅度

angY=angle(Y)*180/pi;%求得FFT变换后的相位 f=n*fs/N;%频率序列 subplot(3,1,2);

stem(f,magY);%画出幅频响应 xlabel('f'); ylabel('幅度'); title('N=512'); grid;

subplot(3,1,3);

stem(f,angY);%画出相频响应 xlabel('f'); ylabel('相位') title('N=512'); grid; 实验结

果8

：

4、 对矩形波进行频谱分析。矩形波为脉冲宽度为2s，持续时间为-5~5s，采用

矩形脉冲函数rectpuls(t,2); 四、实验报告：

通过本实验加深了对FFT的理解，更熟练的应用Matlab软件，提高了自己编写相关程序的能力，熟悉并应用FFT对典型信号进行频谱分析的方法，学习了用FFT对连续信号和时域离散信号进行谱分析的方法，在实验的过程中出现了些许错误，有操作上也有程序上的，通过本实验加强了自己动手纠错的能力。

9