用MATLAB进行FFT频谱分析

用MATLAB 进行FFT 频谱分析

假设一信号：

()()292.7/2cos 1.0996.2/2sin 1.06.0+++=t t R ππ

画出其频谱图。 分析：

首先，连续周期信号截断对频谱的影响。

DFT 变换频谱泄漏的根本原因是信号的截断。即时域加窗，对应为频域卷积，因此，窗函数的主瓣宽度等就会影响到频谱。

实验表明，连续周期信号截断时持续时间与信号周期呈整数倍关系时，利用DFT 变换可以得到精确的模拟信号频谱。举一个简单的例子：

()ππ2.0100cos +=t Y

其周期为0.02。截断时不同的持续时间影响如图一.1：（对应程序shiyan1ex1.m ）

图 错误!文档中没有指定样式的文字。.1

140.0160.0180.02

截断时，时间间期为周期整数倍，频谱图

0.0250.03

20

40

60

80

100

截断时，时间间期不为周期整数倍，频谱图

其次，采样频率的确定。

根据Shannon 采样定理，采样带限信号采样频率为截止频率的两倍以上，给定信号的采样频率应>1/7.92，取16。

再次，DFT 算法包括时域采样和频域采样两步，频域采样长度M 和时域采样长度N 的关系要符合M ≧N 时，从频谱X(k)才可完全重建原信号。

实验中信号R 经采样后的离散信号不是周期信号，但是它又是一个无限长的信号，因此处理时时域窗函数尽量取得宽一些已接近实际信号。

实验结果如图一.2：其中，0点位置的冲激项为直流分量0.6造成（对应程序为shiyan1.m ）

图 错误!文档中没有指定样式的文字。.2

♣ARMA （Auto Recursive Moving Average ）模型：

将平稳随机信号x(n)看作是零均值，方差为σu 2的白噪声u(n)经过线性非移变系统H(z)后的输出，模型的传递函数为

020406080100120140160180200

0.4

0.50.60.7

0.800.050.10.150.20.250.30.350.40.450.5

50100

150

()()()

∑∑=-=-+==P

k k

k Q

r r r

z

a z

b z A z B z H 1

11

用差分方程表示为

()()()∑∑==-+--=Q

r r P k k r n u b k n x a n x 0

1

AR （Auto Recursive ）自回归模型，即ARMA 模型中系数b 只有在r=0的情况下为1，其余都是零，获得一个全极点模型：

()()()

∑=-+==P

k k

k z

a z A z B z H 111

差分方程表示为：

()()()n u k n x a n x P

k k +--=∑=1

AR 模型的功率谱估计为：

()()(

)Ω

-Ω

Ω

=j j u

j x e A e A e

S 1

2σ

程序：%%------------------------------------------------------------------------

%%功能：利用MATLAB 的FFT 函数做双正弦信号频谱分析

%%------------------------------------------------------------------------ fs=16; t=0:1/fs:200;

x6=0.6+sin(2*pi*t/2.996)*0.1+cos(2*pi*t/7.92+2)*0.1;

subplot(2,1,1); plot(t,x6);

N=length(t); subplot(212);

plot((-N/2:N/2-1)*fs/N,abs(fftshift(fft(x6,N)))) %绘制信号的频谱，横轴对应实际频率axis([0 0.5 0 160]);

例子：

%%------------------------------------------------------------------------

%%功能：连续周期信号截断对频谱的影响

%%------------------------------------------------------------------------

fs=8000;

n1=0.02;

n=0:1/fs:n1;

n=n(1,1:end-1);

N=length(n);

y=cos(100*pi*n+0.2*pi);

subplot(2,2,1);

plot(n,y);

title('函数y=cos(100{\pi}t+0.2{\pi})');

subplot(2,2,2);

stem((-N/2:N/2-1)*fs/N,abs(fftshift(fft(y,N))));

axis([0 1000 0 100]);

grid on;

title('截断时，时间间期为周期整数倍，频谱图');

n1=0.03;

n=0:1/fs:n1;

n=n(1,1:end-1);

N=length(n);

y=cos(100*pi*n+0.2*pi);

subplot(2,2,3);

plot(n,y);

title('函数y=cos(100{\pi}t+0.2{\pi})');

subplot(2,2,4);

stem((-N/2:N/2-1)*fs/N,abs(fftshift(fft(y,N))));

axis([0 1000 0 100]);

grid on;

title('截断时，时间间期不为周期整数倍，频谱图');