数字信号处理实验
第五次实验
应用FFT 实现信号的频谱分析
学 号:高 超 姓 名:12081311 指导老师:黄怡、杨萌 选课时间:周一3-5节 实验时间:2014年11月24日
一. 实验目的
(1) 能够熟练掌握快速离散傅里叶变换的原理及应用FFT 进行频谱分析的基本方法。 (2) 对离散傅里叶变换的主要性质及FFT 在数字信号处理的重要作用有进一步的了解。
二. 基本原理
1. 离散傅里叶变换(DFT) 及其主要性质
DFT 表示离散信号的离散频谱,DFT 的主要性质中有奇偶对称特性、虚实特性等。通过实验可以加深理解。
实序列的DFT 具有偶对称的实部和奇对称的虚部, 也即:
X k =X∗(N−k)
实序列DFT 的这个特性,在本实验中可以通过实指数序列及三角序列看出来。
对于单一频率的三角序列来说,他的DFT 谱线也是单一的,这个物理意义可以从实验中得到验证,在理论上可以推到如下:
设
2π
x n =(n) RN(n)
其DFT 为
N−1
X k = x n e
n=0
2π−j2π2π1−j k−1 n−jk+1)n
= (e−e) N−1
n=0
从而,X 0 =
j1 N− en=02j
N−1
1
2π
−e
−j
2π
n =0
4π1NN−jnX 1 = 1−e ==−jn=0
X 2 =0
······
X N−2 =0
2π1NN−j(N−2) nj2πnNX N−1 = (e−e) =−=jN−1
n=0
以上这串式中X(0)反映了x(n) 的直流分量,X(1)是x(n) 的一次谐波,又根据虚实特性
X∗ N−1 =X(1),而其他分量均为零。
当周期减小时显然(
2π
n) RN(n) 的谱只应该在kN 3
=3以及k=N−3才有分量,实验者可
以通过和上述相同的步骤加以理论证明。
由于cos N ·RN(n) 和sin N ·RN(n) 相位差2,所以它的DFT 只包括实部而没有虚部,以上性质在本次实验中可得到验证。
2. 利用DFT对信号进行频谱分析
DFT 的重要应用之一是对时域连续信号的频谱进行分析,称为傅里叶分析,时域连续信号离散傅里叶分析的基本步骤如下图5 – 1所示
2π
2π
π
图5 – 1 时域连续信号离散傅里叶分析步骤
其中消混频低通滤波器LPF (预滤波器)的引入,是为了消除或减少时域连续信号转换成序列时可能出现的频谱混叠的影响。实际工作中,时域离散信号x(n) 的时宽是很长的甚至是无限长。由于DFT 的需要,必须把x(n) 限制在一定的时间间隔之内,即进行数据截断。数据的截断相当于加窗处理。因此,在计算x(n) 的DFT 之前,用一个时域有限的窗函数w(n) 加到x(n) 上是非常有必要的。
在实际应用中,消混叠低通滤波器的阻带不可能无限衰减,故有x(n) 的频谱X(jΩ) 周期延拓得到的X(ejω) 有非零混叠,即出现混叠现象。
加窗对频域的影响,用周期卷积表示
1π
jω
V e = X ejθ W(ej(ω−θ) ) dθ
−π
最后是进行DFT 运算。加窗后的DFT 为
N−1
V k = v n e
m=0
−j
2πnk , 0≤n≤N−1
其中假设窗函数长L小于或等于DFT 长度N。
有限长序列v n =x n ·w(n) 的DFT 相当于v n 傅里叶变换的等间隔取样。
V k =V(ejω)|ω=2πk
N
因为DFT 频率间隔为N,且模拟频率Ω和数字频率ω间的关系为ω=ΩT,所以离散频率点对应的模拟频率为Ω=
2πk
2π
∆f=NT。 NT
2πN
1
利用DFT 计算频谱,只给出频谱ωk=
k或Ωk=NTk的频率分量,即频率的取样值,而不
2π
可能得到连续的频谱函数。
如果在两个离散的谱线之间有一个特别大的频谱分量,就无法检测出来了。
为了在保持原来频谱形状不变的情况下,使频谱加密,即使频域取样点数增加,从而使原来看不到的频谱分量变得可以看到,可以通过在信号数据的末端补加一些零值点,使DFT 计算周期内点数增加,但又不改变原有的记录数据的方法来实现。 3. 快速离散傅里叶变换(FFT)
快速离散傅里叶变换是计算离散傅里叶变换的一种快速算法,为了提高运算速度,FFT 将DFT 的计算逐次分解成较小点数的DFT 。按时间抽取(DIT ) FFT算法把输入序列x(n) 按其n值为偶数或是奇数分解成越来越短的序列。按频域抽取(DIF ) FFT算法把输出序列X(k) 按其k值得奇偶分解成越来越短的序列。
三. 实验内容及实验结果
1. 实验内容
(1) 编写一个调用FFT 函数的通用程序,可以计算下列三种序列的离散频谱。
指数序列:v1(n) =(0.9)nu(n) 。
周期为N的余弦序列:v2 n =cos Nn ,且0≤n≤N−1。 复合函数序列:v3 n =0.9sin Nn +0.6sin N/3n 。
(2) 计算实指数序列v1(n) 的N点离散频谱V1 k ,记录N为不同的2的幂次方时的V1 k 值,
2π
2π
2π
并与理论值V1 ejωk 进行分析比较。
(3) 计算周期为N的余弦序列v2 n 的N点FFT 、2N点FFT 及(N+2) 点FFT ,记录结果并分
析说明。 (4) 已知信号x t =0.15sin 2πf(2πf3t) ,其中 1t +sin 2πf2t −0.1sin
f1=1Hz, f2=2Hz,f3=3Hz,取样频率为32Hz 编程实现:
32点和64点FFT ,画出其幅度谱。比较两者间的差异,思考实际频率与离散频谱图中横坐标k的对应关系。 2. 实验结果
实验源代码:
1). 编写一个调用FFT 函数的通用程序,可以计算题中三种序列的离散频谱 %计算三种序列的通用函数 clear all clc
N = 100; n = 0 : N-1; xn1 = 0.9.^n;
xn2 = cos(2.*pi.*n./N);
xn3 = 0.9.*sin(2.*pi.*n./N) + 0.6.*sin(2.*pi.*n./(3.*N));
XK1 = fft(xn1, N); magXK1 = abs(XK1); phaXK1 = angle(XK1); subplot(1, 2, 1); plot(n, xn1); xlabel('n'); ylabel('v(n)');
title('v1n N=100'); subplot(1, 2, 2);
k = 0 : length(magXK1)-1; stem(k, magXK1, '.'); xlabel('k');
图5 – 2 指数序列时域谱及其离散频ylabel('|V(k)|');
谱 title('V1n N=100');
%指数序列实验结果如图5 – 2所示,
figure(2);
XK2 = fft(xn2, N); magXK2 = abs(XK2); phaXK2 = angle(XK2); subplot(1, 2, 1); plot(n, xn2); xlabel('n'); ylabel('v(n)');
title('v2n N=100'); subplot(1, 2, 2);
k = 0 : length(magXK2)-1; stem(k, magXK2, '.'); xlabel('k');
图5 – 3 余弦序列时域谱及其离散频ylabel('|V(k)|');
谱 title('V2n N=100');
%余弦序列实验结果如图5 – 3所示,
figure(3)
XK3 = fft(xn3, N); magXK3 = abs(XK3); phaXK3 = angle(XK3); subplot(1, 2, 1); plot(n, xn3); xlabel('n'); ylabel('v(n)');
title('v3n N=100'); subplot(1, 2, 2);
k = 0 : length(magXK3)-1; stem(k, magXK3, '.'); xlabel('k');
图5 – 4 复合函数序列时域谱及其离散频谱 ylabel('|V(k)|');
title('V3n N=100');
%复合函数序列实验结果如图5 – 4所示,
2). 计算实指数序列v1(n) 的N点离散频谱V1 k ,记录N为不同的2的幂次方时的V1 k 值,并与理论值V1 ejωk 进行分析比较
%实指数序列频谱的相关计算 clear all clc
N = input('序列长度N='); n = 0 : N-1; xn1 = 0.9.^n; XK1 = fft(xn1, N); magXK1 = abs(XK1); phaXK1 = angle(XK1);
k = 0 : length(magXK1)-1; stem(k, magXK1, '.'); xlabel('k');
ylabel('|X(k)|'); title('指数序列');
%实验结果如下图5 – 5、5 – 6、5 – 7所示
图
5 –
5 v1(n) 的32点离散频谱V1 k 图5 – 6 v1(n) 的64点离散频谱V1 k
图5 – 7 v1(n) 的128点离散频谱V1 k
实验分析:有实验图像可知,对于v1(n) 的N 点离散频谱V1 k ,N 的取值不同,会导致频谱不同:频率最大值不同,各个点的取值情况不同,但是频谱的整体的变化趋势是一致的。
3). 计算周期为N 的余弦序列v_2 (n)的N 点FFT 、2N 点FFT 及(N+2)点FFT ,记录结果并分析说明。
%复合函数序列频谱的相关计算 clear all clc
N = input('序列长度N='); a = 2; b = 2;
n = 0 : N-1;
xn2 = cos(2.*pi.*n./N); XK2 = fft(xn2, N); XK3 = fft(xn2, a*N);
XK4 = fft(xn2, N+b);
figure;
magXK2 = abs(XK2); phaXK2 = angle(XK2); k = 0 : length(magXK2)-1; stem(k, magXK2, '.'); xlabel('k');
ylabel('|X(k)|');
title('周期为N 的余弦序列');
%周期为N 的余弦序列频谱如图5 – 8所示
图5 – 8 v2(n) 的100点离散频谱V2 k
figure(2);
magXK3 = abs(XK3); phaXK3 = angle(XK3); k = 0 : length(magXK3)-1; stem(k, magXK3, '.'); xlabel('k');
ylabel('|X(k)|');
title('周期为2N 的余弦序列');
%周期为2N 的余弦序列频谱如图5 – 9所示
图5 – 9 v2(n) 的200点离散频谱V2 k
figure(3);
magXK4 = abs(XK4); phaXK4 = angle(XK4); k = 0 : length(magXK4)-1; stem(k, magXK4, '.'); xlabel('k');
ylabel('|X(k)|');
title('周期为N+2的余弦序列');
%周期为N 的余弦序列频谱如图5 – 10所示 图5 – 10 v2(n) 的102点离散频谱 V2 k
实验分析:此处取周期N=100。当采样周期M=N
时,在进行一个周期内只取被采样信号的两个点,且两点间距一个周期N ,频谱为单一的谱线,其频谱分辨率F =1Hz 。而对于图5 – 9和图5 – 10采样点分别为200和102,均大于100,因此采样周期小于原信号周期,有原信号失真,造成频谱泄露,从而降低了频谱的分辨率。此外,由于在主频线两边形成很多旁瓣,引起不同频率分量间的干扰
4). 求其32点和64点FFT ,画出其幅度谱。比较两者间的差异,思考实际频率与离散频谱图中横坐标k的对应关系。
%复合函数序列频谱的相关计算 clear all; clc;
f1 = 1;
f2 = 2; f3 = 3;
N = input('序列长度N='); M = input('序列长度M=');
n = 0 : N-1;
x = 0.15*sin(2*pi*f1*n/32) + sin(2*pi*f2*n/32)-0.1*sin(2*pi*f3*n/32); XK = fft(x, N); magXK = abs(XK); phaXK = angle(XK);
k = 0 : length(magXK)-1; subplot(1, 2, 1); stem(k, magXK, '.'); xlabel('k'); ylabel('X|k|'); title('X1(k)');
YK = fft(x, M); magYK = abs(YK); phaYK = angle(YK);
k = 0 : length(magYK)-1; subplot(1, 2, 2); stem(k, magYK, '.'); xlabel('k');
ylabel('X2|k|'); title('X2(k)');
%实验结果如图5–11所示
实验分析:因N=32,故频率分辨率F=fs N=1Hz 。由图X1(k) 可知,k=1,2,3所对应的频谱即为频率f1=1Hz, f2=2Hz,f3=3Hz的正弦波所对应的频谱。而信号的模拟频率为
图5 – 11 x(t) 的32点和 64点离散频谱比较
Ω和数字频率ω之间的关系为ω=ΩT。所以,离散的频率函数第k点对应的模拟频率为
ω2πkkΩk==,fk=
而数字域频率间隔∆ω=2π N对应的模拟域谱线间距为
F=
1fs1==fk−fk−1= p
其中N(tp) 为记录长度(或样本长度) 点数,T为采样时间间隔,F为频谱分辨率。
四. 实验总结
通过本次实验,我掌握了应用FFT 实现信号频谱分析的方法,对于离散傅里叶变换的主要
性质以及FFT 在数字信号处理中的重要作用有了进一步的了解。
我还注意到,在对信号求FFT 时,容易出现频谱混叠失真、栅栏效应以及频谱泄露(谱间干扰) 等误差现象。为避免以上现象出现,需要
1. 对于频谱混叠失真,为兼顾高频容量fℎ与频率分辨率F的唯一办法就是增加记录长度的点数N,即要满足
f2fℎsN=>
这个条件是在为采样任何特殊数据处理(加窗处理)的情况下,为实现基本的FFT 算法所必须满足的最低条件。如果进行了加窗处理,必然加宽频谱分量,频谱分辨率就
可能变差,为了保证频率分辨率不变,则必须增加数据长度tp。
2. 为减小栅栏效应,可以使频域采样更密,即增加频率采样点数N(末尾补零) ,而不改变tp。并且,补零之后的窗函数长度,必须和补零前一致。
3. 在进行FFT 运算时,时域截断是必然的,从而频谱泄露和谱间干扰是不可避免的。为尽量减小泄露和谱间干扰的影响,需增加窗的时域宽度(频域主瓣变窄) 。一般选用各种缓变的窗,例如三角形窗、升余弦窗、改进的升余弦窗等,使得窗谱的旁瓣能力更小,卷积后造成的泄露减小。
数字信号处理实验
第五次实验
应用FFT 实现信号的频谱分析
学 号:高 超 姓 名:12081311 指导老师:黄怡、杨萌 选课时间:周一3-5节 实验时间:2014年11月24日
一. 实验目的
(1) 能够熟练掌握快速离散傅里叶变换的原理及应用FFT 进行频谱分析的基本方法。 (2) 对离散傅里叶变换的主要性质及FFT 在数字信号处理的重要作用有进一步的了解。
二. 基本原理
1. 离散傅里叶变换(DFT) 及其主要性质
DFT 表示离散信号的离散频谱,DFT 的主要性质中有奇偶对称特性、虚实特性等。通过实验可以加深理解。
实序列的DFT 具有偶对称的实部和奇对称的虚部, 也即:
X k =X∗(N−k)
实序列DFT 的这个特性,在本实验中可以通过实指数序列及三角序列看出来。
对于单一频率的三角序列来说,他的DFT 谱线也是单一的,这个物理意义可以从实验中得到验证,在理论上可以推到如下:
设
2π
x n =(n) RN(n)
其DFT 为
N−1
X k = x n e
n=0
2π−j2π2π1−j k−1 n−jk+1)n
= (e−e) N−1
n=0
从而,X 0 =
j1 N− en=02j
N−1
1
2π
−e
−j
2π
n =0
4π1NN−jnX 1 = 1−e ==−jn=0
X 2 =0
······
X N−2 =0
2π1NN−j(N−2) nj2πnNX N−1 = (e−e) =−=jN−1
n=0
以上这串式中X(0)反映了x(n) 的直流分量,X(1)是x(n) 的一次谐波,又根据虚实特性
X∗ N−1 =X(1),而其他分量均为零。
当周期减小时显然(
2π
n) RN(n) 的谱只应该在kN 3
=3以及k=N−3才有分量,实验者可
以通过和上述相同的步骤加以理论证明。
由于cos N ·RN(n) 和sin N ·RN(n) 相位差2,所以它的DFT 只包括实部而没有虚部,以上性质在本次实验中可得到验证。
2. 利用DFT对信号进行频谱分析
DFT 的重要应用之一是对时域连续信号的频谱进行分析,称为傅里叶分析,时域连续信号离散傅里叶分析的基本步骤如下图5 – 1所示
2π
2π
π
图5 – 1 时域连续信号离散傅里叶分析步骤
其中消混频低通滤波器LPF (预滤波器)的引入,是为了消除或减少时域连续信号转换成序列时可能出现的频谱混叠的影响。实际工作中,时域离散信号x(n) 的时宽是很长的甚至是无限长。由于DFT 的需要,必须把x(n) 限制在一定的时间间隔之内,即进行数据截断。数据的截断相当于加窗处理。因此,在计算x(n) 的DFT 之前,用一个时域有限的窗函数w(n) 加到x(n) 上是非常有必要的。
在实际应用中,消混叠低通滤波器的阻带不可能无限衰减,故有x(n) 的频谱X(jΩ) 周期延拓得到的X(ejω) 有非零混叠,即出现混叠现象。
加窗对频域的影响,用周期卷积表示
1π
jω
V e = X ejθ W(ej(ω−θ) ) dθ
−π
最后是进行DFT 运算。加窗后的DFT 为
N−1
V k = v n e
m=0
−j
2πnk , 0≤n≤N−1
其中假设窗函数长L小于或等于DFT 长度N。
有限长序列v n =x n ·w(n) 的DFT 相当于v n 傅里叶变换的等间隔取样。
V k =V(ejω)|ω=2πk
N
因为DFT 频率间隔为N,且模拟频率Ω和数字频率ω间的关系为ω=ΩT,所以离散频率点对应的模拟频率为Ω=
2πk
2π
∆f=NT。 NT
2πN
1
利用DFT 计算频谱,只给出频谱ωk=
k或Ωk=NTk的频率分量,即频率的取样值,而不
2π
可能得到连续的频谱函数。
如果在两个离散的谱线之间有一个特别大的频谱分量,就无法检测出来了。
为了在保持原来频谱形状不变的情况下,使频谱加密,即使频域取样点数增加,从而使原来看不到的频谱分量变得可以看到,可以通过在信号数据的末端补加一些零值点,使DFT 计算周期内点数增加,但又不改变原有的记录数据的方法来实现。 3. 快速离散傅里叶变换(FFT)
快速离散傅里叶变换是计算离散傅里叶变换的一种快速算法,为了提高运算速度,FFT 将DFT 的计算逐次分解成较小点数的DFT 。按时间抽取(DIT ) FFT算法把输入序列x(n) 按其n值为偶数或是奇数分解成越来越短的序列。按频域抽取(DIF ) FFT算法把输出序列X(k) 按其k值得奇偶分解成越来越短的序列。
三. 实验内容及实验结果
1. 实验内容
(1) 编写一个调用FFT 函数的通用程序,可以计算下列三种序列的离散频谱。
指数序列:v1(n) =(0.9)nu(n) 。
周期为N的余弦序列:v2 n =cos Nn ,且0≤n≤N−1。 复合函数序列:v3 n =0.9sin Nn +0.6sin N/3n 。
(2) 计算实指数序列v1(n) 的N点离散频谱V1 k ,记录N为不同的2的幂次方时的V1 k 值,
2π
2π
2π
并与理论值V1 ejωk 进行分析比较。
(3) 计算周期为N的余弦序列v2 n 的N点FFT 、2N点FFT 及(N+2) 点FFT ,记录结果并分
析说明。 (4) 已知信号x t =0.15sin 2πf(2πf3t) ,其中 1t +sin 2πf2t −0.1sin
f1=1Hz, f2=2Hz,f3=3Hz,取样频率为32Hz 编程实现:
32点和64点FFT ,画出其幅度谱。比较两者间的差异,思考实际频率与离散频谱图中横坐标k的对应关系。 2. 实验结果
实验源代码:
1). 编写一个调用FFT 函数的通用程序,可以计算题中三种序列的离散频谱 %计算三种序列的通用函数 clear all clc
N = 100; n = 0 : N-1; xn1 = 0.9.^n;
xn2 = cos(2.*pi.*n./N);
xn3 = 0.9.*sin(2.*pi.*n./N) + 0.6.*sin(2.*pi.*n./(3.*N));
XK1 = fft(xn1, N); magXK1 = abs(XK1); phaXK1 = angle(XK1); subplot(1, 2, 1); plot(n, xn1); xlabel('n'); ylabel('v(n)');
title('v1n N=100'); subplot(1, 2, 2);
k = 0 : length(magXK1)-1; stem(k, magXK1, '.'); xlabel('k');
图5 – 2 指数序列时域谱及其离散频ylabel('|V(k)|');
谱 title('V1n N=100');
%指数序列实验结果如图5 – 2所示,
figure(2);
XK2 = fft(xn2, N); magXK2 = abs(XK2); phaXK2 = angle(XK2); subplot(1, 2, 1); plot(n, xn2); xlabel('n'); ylabel('v(n)');
title('v2n N=100'); subplot(1, 2, 2);
k = 0 : length(magXK2)-1; stem(k, magXK2, '.'); xlabel('k');
图5 – 3 余弦序列时域谱及其离散频ylabel('|V(k)|');
谱 title('V2n N=100');
%余弦序列实验结果如图5 – 3所示,
figure(3)
XK3 = fft(xn3, N); magXK3 = abs(XK3); phaXK3 = angle(XK3); subplot(1, 2, 1); plot(n, xn3); xlabel('n'); ylabel('v(n)');
title('v3n N=100'); subplot(1, 2, 2);
k = 0 : length(magXK3)-1; stem(k, magXK3, '.'); xlabel('k');
图5 – 4 复合函数序列时域谱及其离散频谱 ylabel('|V(k)|');
title('V3n N=100');
%复合函数序列实验结果如图5 – 4所示,
2). 计算实指数序列v1(n) 的N点离散频谱V1 k ,记录N为不同的2的幂次方时的V1 k 值,并与理论值V1 ejωk 进行分析比较
%实指数序列频谱的相关计算 clear all clc
N = input('序列长度N='); n = 0 : N-1; xn1 = 0.9.^n; XK1 = fft(xn1, N); magXK1 = abs(XK1); phaXK1 = angle(XK1);
k = 0 : length(magXK1)-1; stem(k, magXK1, '.'); xlabel('k');
ylabel('|X(k)|'); title('指数序列');
%实验结果如下图5 – 5、5 – 6、5 – 7所示
图
5 –
5 v1(n) 的32点离散频谱V1 k 图5 – 6 v1(n) 的64点离散频谱V1 k
图5 – 7 v1(n) 的128点离散频谱V1 k
实验分析:有实验图像可知,对于v1(n) 的N 点离散频谱V1 k ,N 的取值不同,会导致频谱不同:频率最大值不同,各个点的取值情况不同,但是频谱的整体的变化趋势是一致的。
3). 计算周期为N 的余弦序列v_2 (n)的N 点FFT 、2N 点FFT 及(N+2)点FFT ,记录结果并分析说明。
%复合函数序列频谱的相关计算 clear all clc
N = input('序列长度N='); a = 2; b = 2;
n = 0 : N-1;
xn2 = cos(2.*pi.*n./N); XK2 = fft(xn2, N); XK3 = fft(xn2, a*N);
XK4 = fft(xn2, N+b);
figure;
magXK2 = abs(XK2); phaXK2 = angle(XK2); k = 0 : length(magXK2)-1; stem(k, magXK2, '.'); xlabel('k');
ylabel('|X(k)|');
title('周期为N 的余弦序列');
%周期为N 的余弦序列频谱如图5 – 8所示
图5 – 8 v2(n) 的100点离散频谱V2 k
figure(2);
magXK3 = abs(XK3); phaXK3 = angle(XK3); k = 0 : length(magXK3)-1; stem(k, magXK3, '.'); xlabel('k');
ylabel('|X(k)|');
title('周期为2N 的余弦序列');
%周期为2N 的余弦序列频谱如图5 – 9所示
图5 – 9 v2(n) 的200点离散频谱V2 k
figure(3);
magXK4 = abs(XK4); phaXK4 = angle(XK4); k = 0 : length(magXK4)-1; stem(k, magXK4, '.'); xlabel('k');
ylabel('|X(k)|');
title('周期为N+2的余弦序列');
%周期为N 的余弦序列频谱如图5 – 10所示 图5 – 10 v2(n) 的102点离散频谱 V2 k
实验分析:此处取周期N=100。当采样周期M=N
时,在进行一个周期内只取被采样信号的两个点,且两点间距一个周期N ,频谱为单一的谱线,其频谱分辨率F =1Hz 。而对于图5 – 9和图5 – 10采样点分别为200和102,均大于100,因此采样周期小于原信号周期,有原信号失真,造成频谱泄露,从而降低了频谱的分辨率。此外,由于在主频线两边形成很多旁瓣,引起不同频率分量间的干扰
4). 求其32点和64点FFT ,画出其幅度谱。比较两者间的差异,思考实际频率与离散频谱图中横坐标k的对应关系。
%复合函数序列频谱的相关计算 clear all; clc;
f1 = 1;
f2 = 2; f3 = 3;
N = input('序列长度N='); M = input('序列长度M=');
n = 0 : N-1;
x = 0.15*sin(2*pi*f1*n/32) + sin(2*pi*f2*n/32)-0.1*sin(2*pi*f3*n/32); XK = fft(x, N); magXK = abs(XK); phaXK = angle(XK);
k = 0 : length(magXK)-1; subplot(1, 2, 1); stem(k, magXK, '.'); xlabel('k'); ylabel('X|k|'); title('X1(k)');
YK = fft(x, M); magYK = abs(YK); phaYK = angle(YK);
k = 0 : length(magYK)-1; subplot(1, 2, 2); stem(k, magYK, '.'); xlabel('k');
ylabel('X2|k|'); title('X2(k)');
%实验结果如图5–11所示
实验分析:因N=32,故频率分辨率F=fs N=1Hz 。由图X1(k) 可知,k=1,2,3所对应的频谱即为频率f1=1Hz, f2=2Hz,f3=3Hz的正弦波所对应的频谱。而信号的模拟频率为
图5 – 11 x(t) 的32点和 64点离散频谱比较
Ω和数字频率ω之间的关系为ω=ΩT。所以,离散的频率函数第k点对应的模拟频率为
ω2πkkΩk==,fk=
而数字域频率间隔∆ω=2π N对应的模拟域谱线间距为
F=
1fs1==fk−fk−1= p
其中N(tp) 为记录长度(或样本长度) 点数,T为采样时间间隔,F为频谱分辨率。
四. 实验总结
通过本次实验,我掌握了应用FFT 实现信号频谱分析的方法,对于离散傅里叶变换的主要
性质以及FFT 在数字信号处理中的重要作用有了进一步的了解。
我还注意到,在对信号求FFT 时,容易出现频谱混叠失真、栅栏效应以及频谱泄露(谱间干扰) 等误差现象。为避免以上现象出现,需要
1. 对于频谱混叠失真,为兼顾高频容量fℎ与频率分辨率F的唯一办法就是增加记录长度的点数N,即要满足
f2fℎsN=>
这个条件是在为采样任何特殊数据处理(加窗处理)的情况下,为实现基本的FFT 算法所必须满足的最低条件。如果进行了加窗处理,必然加宽频谱分量,频谱分辨率就
可能变差,为了保证频率分辨率不变,则必须增加数据长度tp。
2. 为减小栅栏效应,可以使频域采样更密,即增加频率采样点数N(末尾补零) ,而不改变tp。并且,补零之后的窗函数长度,必须和补零前一致。
3. 在进行FFT 运算时,时域截断是必然的,从而频谱泄露和谱间干扰是不可避免的。为尽量减小泄露和谱间干扰的影响,需增加窗的时域宽度(频域主瓣变窄) 。一般选用各种缓变的窗,例如三角形窗、升余弦窗、改进的升余弦窗等,使得窗谱的旁瓣能力更小,卷积后造成的泄露减小。