数字信号处理课程设计正余弦信号的谱分析

正余弦信号的谱分析

中文摘要：使用MATLAB 软件，通过编写程序，对正余弦信号进行傅里叶变换。用DFT 和FFT 实现对正余弦信号的谱分析，并且分析DFT 长度对频谱的影响。

关键词：matlab ； 傅里叶变换； DFT; FFT;

一、概述

数字信号处理方法的一个重要用途是在离散时间域中确定一个连续时间信号的频谱，通常称为频谱分析，更具体的说，它也包括确定能量谱或功率谱。数字频谱分析可以应用在很广阔的领域。

二、设计目的

1. 用DFT 实现对正余弦信号的谱分析；

2. 观察DFT 长度和窗函数长度对频谱的影响；

3. 对DFT 进行谱分析中的误差现象获得感性认识。

三、设计原理

1、谱分析原理

频谱分析方法是基于以下的观测：如果连续时间信号g a (t ) 是频带有限的，那么对其离散时间等效信号g (n ) 的DFT 进行谱分析。它的离散时间等效物g(n)应当能给出g a (t)频谱的一个很近似的估计两者之间只差一个带数因子T 。然而，在大多数情况下，g a (t ) 是在-∞

假定表征正余弦信号的基本参数，如振幅、频率和相位不随时间改变，则此信号的傅立叶变换G (e j ω) 可以用计算它的DTFT 得到

G (e j ω) =∑g (n ) e

n =-∞∞-j ωn （1） 实际上无限长序列g (n ) 首先乘以一个长度为M 的窗函数w (n ) ，使它变成一个长为M 的有限长序列，

g 1(n ) =g (n ) w (n ) ，对g 1(n ) 求出的DTFT G 1(e j ω) 应该可以作为原连续模拟信号g a (t ) 的频谱估计，然后求出G 1(e j ω) 在0≤ω≤2π区间等分为N 点的离散傅立叶变换DFT 。为保证足够的分辨率，DFT 的长度N 选的比窗长度M 大，其方法是在截断了的序列后面补上N －M 个零。计算采用FFT 算法。 更详细地考察一下上面的方法。这样才能了解它的限制，并正确利用它所得出的结果。特别要分析加窗的效果，以及和由DFT 样本来估计DTFT 频率采样值的问题。

在讨论由G 1(k ) 来估计频谱G 1(e ) 和G (e ) 时，需要重新探讨一下这些变换和它们所对应的频率之间的关系，R 点的DFT G 1(k ) 与它的DTFT G 1(e ) 的关系为： jw jw jw

G 1(k ) =G 1(e jw ) |w =2πk /R 0≤k ≤R -1 (2) 归一化的数字角频率ωk 和DFT 样本序号k 的关系为

ωk =2πk (3) R

同样的，模拟角频率和DFT 样本序号k 的关系为

Ωk =2πk (4) RT

其中T 是采样周期。为了能正确地解释一DFT 为基础的频谱分析结果，首先来考虑单频率正余弦序列的频谱分析。设一个具有数字角频率ω0的余弦信号为：

g (n ) =cos(ω0n +φ) (5)

把这个序列表为： g (n ) =

查表得知它的DTFT 为：

jw 1j (ω0n +φ) (e +e -j (ω0n +φ) ) (6) 2j φ-j φG (e ) =π[e δ(ω-ω+2πl ) +e δ(ω+ω0+2πl )] (7) ∑0∞

因此，它是一个以2π为周期的ω的周期信号，每个周期中包含两个冲击信号。在-π≤ω≤π的频率范围内，ω=ω0处的冲击具有复数幅特性πe j ϕt =-∞，而在ω=-ω0处具有πe -j ϕ的复数幅特性。

为了用DFT 分析个g(n)的频域特性，取g(n)的一个有限长序列

g 1(n ) =cos(ω0n +φ) 0≤n ≤N -1 (8)

2、MATLAB 函数介绍

1. 输入函数input( )

格式：R=input(string)

功能：在屏幕上显示input 括号后的’string ’内容，提示用户从键盘输入某值，并将输入的值赋给R 。

例如，在命令窗口输入R=input(‘How many apples’)

会显示How many apples

从键盘输入 3

会显示 R=3

2. 一维快速傅里叶变换函数fft( )

格式： y=fft(x)

y=fft(x,n)

说明：fft 函数用于计算矢量或矩阵的傅里叶变换。

３．幅度函数abs( )

格式：abs （X ）

功能：对X 取绝对值，当X 是复数时，得到X 的复模值。

四、设计过程

1. 理解题目，编写程序；

2. 运行程序；

3. 输入相关参数；

4. 观察图像。

设计过程框图：

五、设计内容

1. 对一个频率为10Hz ，采样频率为64Hz 的32点余弦序列进行谱分析，画出其频谱图；若将频率改为11Hz ，其他参数不变，重新画出该序列的频谱图，观察频率泄漏现象，分析原因。

（1）. 频率为10Hz ，采样频率为64Hz 的32点余弦序列

clear all;close all;clc;

N=32;n=0:N-1;Fs=64;T=1/Fs;

x1=cos(10*2*pi*n*T); %定义10Hz 序列x1

k=0:N-1;

X1=abs(fft(x1,N)); %求余弦序列x1的32点FFT

subplot(1,2,1);stem(n,x1); %绘制10Hz 序列x1的波形

xlabel('n');ylabel('x1(n)'); title('余弦序列(f=10Hz)');

subplot(1,2,2);stem(k,X1); %绘制序列x1的幅频特性曲线

xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线

(f=10Hz)');

图1. 频率为10Hz ，采样频率为64Hz 的32点余弦序列频谱图

（2）频率为11Hz ，采样频率为64Hz 的32点余弦序列

clear all;close all;clc;

N=32;n=0:N-1;Fs=64;T=1/Fs;

x2=cos(11*2*pi*n*T); %定义11Hz 序列x2

k=0:N-1;

X2=abs(fft(x2,N)); %求余弦序列x2的32点FFT

subplot(1,2,1);stem(n,x2); %绘制11Hz 序列x2的波形

xlabel('n');ylabel('x2(n)'); title('余弦序列(f=11Hz)');

subplot(1,2,2);stem(k,X2); %绘制序列x2的幅频特性曲线

xlabel('k');ylabel('X2(k)');title('32点FFT 幅频曲线

(f=11Hz)');

图2. 频率为11Hz ，采样频率为64Hz 的32点余弦序列频谱图

分析：

由图1可见，它的DFT 只有两个点不等于零，这样DFT 确实正确地分辨了余弦信号的频率。但是这样理想的结果是碰巧得到的，因为我们恰好截取了五个完整的余弦周期(f*N/Fs=5)。如果截取的不是整数周期，情况就不同了。把频率改为11Hz ，而采样频率仍为64Hz 和窗长度仍为32点，用同样的程序计算此余弦信号的频谱，可以得到，频谱图上有两个较大的峰值，而其他的点上幅度也不再为零。，因此可以判断，信号的频谱峰值确实位于两者之间。本来是单一的11Hz 频率的能量会分布到许多DFT 频率上的现象称为频率泄漏，来源于截断效应。

2. 考察DFT 的长度对双频率信号频谱分析的影响。设待分析的信号为

x (n ) =0. 5sin(2πf 1n ) +sin(2πf 2n ) 0≤n ≤15

令两个长度为16的正余弦序列的数字频率为f 1=0. 22及f 2=0. 34。取N 为四个不同值16，32，64，128。画出四个DFT 幅频图，分析DFT 长度对频谱分辨率的影响。

n=0:15; N1=16;N2=32;N3=64;N4=128;f1=0.22;f2=0.34;

x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n); %待分析信号

%谱分析

X1=abs(fft(x,N1));

X2=abs(fft(x,N2));

X3=abs(fft(x,N3));

X4=abs(fft(x,N4));

k1=0:N1-1;k2=0:N2-1;k3=0:N3-1;k4=0:N4-1;

subplot(4,2,1),stem(n,x); %绘制序列的波形

xlabel('n');ylabel('x1(n)');title('余弦序列');

subplot(4,2,2);stem(k1,X1); %绘制N1的幅频特性

xlabel('n');ylabel('X(n)'); title('N=16的幅频特性曲线图');

subplot(4,2,3),stem(n,x); %绘制序列的波形

xlabel('n');ylabel('x1(n)');title('余弦序列');

subplot(4,2,4);stem(k2,X2); %绘制N2的幅频特性

xlabel('n');ylabel('X(n)'); title('N=32的幅频特性曲线图');

subplot(4,2,5),stem(n,x); %绘制序列的波形

xlabel('n');ylabel('x1(n)');title('余弦序列');

subplot(4,2,6);stem(k3,X3); %绘制N3的幅频特性

xlabel('n');ylabel('X(n)'); title('N=64的幅频特性曲线图');

subplot(4,2,7),stem(n,x); %绘制序列的波形

xlabel('n');ylabel('x1(n)');title('余弦序列');

subplot(4,2,8);stem(k4,X4); %绘制N4的幅频特性

xlabel('n');ylabel('X(n)'); title('N=128的幅频特性曲线图

')

图3. N=16的DFT 幅频图

图4. N=32的DFT 幅频图

图5. N=64的DFT 幅频图

图6. N=128的DFT 幅频图

分析：

采样速率Fs>2fc，谱分辨率F=Fs/N，如果保持采样点数N 不变，要提高频谱分辨率（减小F ），就必须降低采样频率，采样频率的降低会引起谱分析范围变窄和谱频率混叠失真。如维持Fs 不变，为提高频谱分辨率可以增加采样点数N 。更长的时域信号能够提供更高的频域分辨率，因为一个N 点的时域信号能被分解为N/2+1个余弦信号和N/2+1个正弦信号，N 增大则(N/2+1)也增大，频域间隔(1/2的时域采样频率)/(N/2+1)减小，所以频域分辨率提高了。所以利用DFT 计算频谱时增加取样点的长度范围可以提高分辨率。

3. 在上题中若把两个正弦波的频率取得较近，令f 1=0. 22，f 2=0. 25，试问怎样选择FFT 参数才能在频谱分析中分辨出这两个分量？

N=input('输入谱分析的长度');

n=1:N-1;

f1=0.22;f2=0.25;

x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n); %待分析信号

subplot(1,2,1),stem(n,x); %绘制序列的波形

xlabel('n');ylabel('x1(n)');

title('余弦序列');

X=abs(fft(x,N)); %求余弦序列的N 点FFT

subplot(1,2,2);

k=0:N-1;

stem(k,X); %绘制序列的幅频特性曲线

xlabel('k');ylabel('X(k)');

string=[num2str(N),'点FFT 幅频曲线'];

title(string);

图7. N=16的FFT 幅频曲线

图8. N=32的FFT 幅频曲线

图9. N=64的FFT 幅频曲线

图10. N=128的FFT 幅频曲线

分析：

从题可以看出当Np=1/|f1-f2|=1/(0.25-0.22)约为33，而N=32最接近33，因为32点的分辨率最好。只有当Np=2*pi/|f1-f2|时才能很好的恢复原信号的频谱。因而，对于双频信号而言，要增大其频

谱的分辨率不仅与其窗函数的宽度N 有关还与两频率差有关。因此，随着N 的增大，其频谱的分辨率增大，但是当N 增加到一定时，其频谱分辨率反而下降。因而不能通过无线的增大N 来扩大频谱的分辨率。其第2题中Np=1/(0.34-0.22)约等于8，因而四个中16点的频谱的分辨率最好。

4. 扩展练习：

DFT 长度对单频率信号频谱分析的影响。

对连续的单一频率周期信号按采样频率f =8f 1采样，截取长度N ，分别选N=20和N=16，观察其DFT 结果的幅度谱。此时离散序列x (n ) =sin(2πnf 1/f ) =sin(2πn /8) ，即k=8。

N=input('输入谱分析的长度');

n=1:N-1;

f1=1/8;

x=sin(2*pi*f1*n); %待分析信号

subplot(1,2,1),stem(n,x); %绘制序列的波形

xlabel('n');ylabel('x1(n)');

title('余弦序列');

X=abs(fft(x,N)); %求余弦序列的N 点DFT

subplot(1,2,2);

k=0:N-1;

stem(k,X); %绘制序列的幅频特性曲线

xlabel('k');ylabel('X(k)');

string=[num2str(N),'点DFT 幅频曲线'];

title(string);

图11. N=20的单频DFT 频谱曲线

图12. N=16的单频DFT 频谱曲线

分析：

计算结果示于图11图12，图11分别是N=20时的截取信号和DFT 结果，由于截取了两个半周期，频谱出现泄漏；图12分别是N=16时的截取信号和DFT 结果，由于截取了两个整周期得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏 。

结束语： 通过此次课程设计，使我更加扎实的掌握了有关数字信号处理以及MATLAB 方面的知识，在设计过程中虽然遇到了一些问题，但经过一次又一次的思考，一遍又一遍的检查以及与同学老师的交流终于找出了原因所在，也暴露出了前期我在这方面的知识欠缺和经验不足。实践出真知，通过亲自动手制作，使我们掌握的知识不再是纸上谈兵。

在课程设计过程中，我们不断发现错误，不断改正，不断领悟，不断获取。最终的检测调试环节，本身就是在践行“过而能改，善莫大焉”的知行观。这次课程设计终于顺利完成了，在设计中遇到了很多问题，最后在老师的指导下，终于游逆而解。在今后社会的发展和学习实践过程中，一定要不懈努力，不能遇到问题就想到要退缩，一定要不厌其烦的发现问题所在，然后一一进行解决，只有这样，才能成功的做成想做的事，才能在今后的道路上劈荆斩棘，而不是知难而退，那样永远不可能收获成功，收获喜悦，也永远不可能得到社会及他人对你的认可。

参考文献：

[1] 张威 . MATLAB 基础与编程入门 . 西安电子科技大学出版社，20008.1

[2] 程佩青 数字信号处理教程 清华大学出版社，2010.5

[3] 曹志刚、钱亚生 现代通信原理 清华大学出版社 2010.6

[4] 王福昌 通信原理学习辅导 华中科技大学出版社 2008.8

[5] 高成 matlab 图像处理与应用 国防工业出版社 2007.4

[6] 宋建华．基于Matlab 的FFT 及其应用[J]．牡丹江师范学院学报(自然科学版) ，2007，2，3—4．