DSP实验报告二.docx
《DSP实验报告二.docx》由会员分享，可在线阅读，更多相关《DSP实验报告二.docx（28页珍藏版）》请在冰豆网上搜索。
DSP实验报告二
实验二应用FFT对信号进行频谱分析
一、实验目的
1、在理论学习的基础上，通过本次实验，加深对快速傅里叶变换的理解，熟悉FFT算法及其程序的编写。
2、熟悉应用FFT对典型信号进行频谱分析的方法。
3、了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。
二、实验原理与方法
①
一个连续信号的频谱可以用它的傅立叶变换表示为
（1）
如果对该信号进行理想采样，可以得到采样序列
（2）
同样可以对该序列进行z变换，其中T为采样周期
（3）
令z为
，则序列的傅立叶变换
（4）
其中ω为数字频率，它和模拟域频率的关系为
（5）
式中的是采样频率。
上式说明数字频率是模拟频率对采样率的归一化。
同模拟域的情况相似，数字频率代表了序列值变化的速率，而序列的傅立叶变换称为序列的频谱。
序列的傅立叶变换和对应的采样信号频谱具有下式的对应关系。
（6）
即序列的频谱是采样信号频谱的周期延拓。
从式（6）可以看出，只要分析采样序列的谱，就可以得到相应的连续信号的频谱。
注意：
这里的信号必须是带限信号，采样也必须满足Nyquist定理。
在各种信号序列中，有限长序列在数字信号处理中占有很重要的地位。
无限长的序列也往往可以用有限长序列来逼近。
有限长的序列可以使用离散傅立叶变换（DFT）。
当序列的长度是N时，定义离散傅立叶变换为：
（7）
其中
，它的反变换定义为：
（8）
根据式（3）和（7），则有
（9）
可以得到
W是z平面单位圆上幅角为
的点，就是将单位圆进行N等分以后第k个点。
所以，X（k）是z变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样。
时域采样在满足Nyquist定理时，就不会发生频谱混淆；同样地，在频率域进行采样的时候，只要采样间隔足够小，也不会发生时域序列的混淆。
②
DFT是对序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。
在运用DFT进行频谱分析的时候可能有三种误差，分析如下：
（1）混淆现象
从式（6）中可以看出，序列的频谱是采样信号频谱的周期延拓，周期是2π/T，因此当采样速率不满足Nyquist定理，即采样频率fs=1/T小于两倍的信号（这里指的是实信号）频率时，经过采样就会发生频谱混淆。
这导致采样后的信号序列频谱不能真实地反映原信号的频谱。
所以，在利用DFT分析连续信号频谱的时候，必须注意这一问题。
避免混淆现象的唯一方法是保证采样的速率足够高，使频谱交叠的现象不出现。
这就告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。
在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。
（2）泄漏现象
实际中的信号序列往往很长，甚至是无限长序列。
为了方便，往往用截短的序列来近似它们。
这样可以使用较短的DFT来对信号进行频谱分析。
这种截短等价于给原信号序列乘以一个矩形窗函数。
而矩形窗函数的频谱不是有限带宽的，从而它和原信号的频谱进行卷积以后会扩展原信号的频谱。
值得一提的是，泄漏是不能和混淆完全分离开的，因为泄露导致频谱的扩展，从而造成混淆。
为了减小泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。
（3）栅栏效应
因为DFT是对单位圆上z变换的均匀采样，所以它不可能将频谱视为一个连续函数，这样就产生了栅栏效应。
减小栅栏效应的一个方法是在源序列的末端补一些零值，从而变动DFT的点数。
这种方法的实质是认为地改变了对真实频谱采样的点数和位置，相当于搬动了“栅栏”的位置，从而使得原来被挡住的一些频谱的峰点或谷点显露出来。
注意，这时候每根谱线多对应的频率和原来的已经不相同了。
从上面的分析过程可以看出，DFT可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减小和消除这些误差的影响。
③
快速傅立叶变换FFT并不是与DFT不相同的另一种变换，而是为了减少DFT运算次数的一种快速算法。
它是对变换式（2-7）进行一次次的分解，使其成为若干小点数DFT的组合，从而减小运算量。
常用的FFT是以2为基数，其长度
。
它的运算效率高，程序比较简单，使用也十分地方便。
当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT，可以用末尾补零的方法，使其长度延长至2的整数次方。
IFFT一般可以通过FFT程序来完成，比较式（2-7）和（2-8），只要对X（k）取共轭，进行FFT运算，然后再取共轭，并乘以因子1/N，就可以完成IFFT。
三、实验步骤
（一）编制实验用主程序及相应子程序
1、在实验之前，认真复习DFT和FFT有关的知识，阅读本实验原理与方法和实验附录部分中和本实验有关的子程序，掌握子程序的原理并学习调用方法。
2、编制信号产生子程序及本实验的频掊分析主程序。
实验中需要用到的基本信号包括：
（1）高斯序列：
（2）衰减正弦序列：
（3）三角波序列：
（4）反三角序列：
四、实验内容及结果分析：
1、观察高斯序列的时域和频域特性
①固定信号中的参数p=8，改变q的值，使q分别等于2，4，8。
观察它们的时域和幅频特性，了解q取不同值的时候，对信号时域特性和幅频特性的影响。
②固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列时域及幅频特性的影响。
注意p等于多少时，会发生明显的泄漏现象，混淆现象是否也随之出现？
记录实验中观察到的现象，绘制相应的时域序列和幅频特性曲线。
结果分析：
①固定p=8时，结果的时域特性：
在时域上，p=8表明最大值1在p=8处取得，q值越大，序列以p为中心衰减速度就越慢。
序列的形状就越平缓。
幅频特性：
上图表明频域的情况恰好相反，p值越大，频谱的带宽就越小，波形越陡峭。
②当固定q=8时，高斯序列的时域特性：
当固定q=8时，序列的形状固定了，p值的不同只是将序列平移了。
但由于序列长度有限，造成信号的截断。
由于q值一样，频谱幅值理应一样（因为序列的形状没变化），但由于信号被截断，当p=14时，产生了明显的频谱泄漏（即截断等价的窗函数引起的频谱搬移）。
与此同时，也产生了混淆（高频分量高于1/2的采样率）。
2、观察衰减正弦序列的时域和幅频特性
①令α=0.1并且f=0.0625，检查谱峰出现的位置是否正确，注意频谱的形状，绘制幅频特性曲线。
②改变f=0.4375，再变化f=0.5625，观察这两种情况下，频谱的形状和谱峰出现的位置，有无混淆和泄漏现象发生？
说明产生现象的原因。
结果分析：
①α=0.1并且f=0.0625时，衰减正弦序列的时域和频域特性如下：
从图中可以看到谱峰出现在n=1和15处，和理论计算一致。
②当α=0.1，f=0.4375和0.5625时的时域图形：
从上面看来，它们时域只相差正负号。
则频域应相差
。
相频如下：
可以看到它们相频确实相差
。
同时从幅频图可以看到产生了混淆和泄漏。
因为f比较大，频谱的主分量比较大，naquist频率超过了25Hz,故产生了混淆和泄漏。
3、观察三角波序列和反三角波序列的时域和幅频特性
①用8点FFT分析信号
和
的幅频特性，观察两者的序列形状和频谱曲线有什么异同？
（注意：
这时候的可以看作是经过圆周移位以后得到的）绘制两者的序列和幅频特性曲线。
②在和的尾补零，用16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？
两个信号之间的FFT频谱还有没有相同之处？
这些变化说明了什么？
结果分析：
1
和
的序列形状和幅频特性分别如下：
可以看到尽管三角序列和反三角序列在时域上时是互为圆周移位形成的，但在幅度频谱上是完全一样的。
区别只是相位的不同：
在和的末尾补零产生的频谱如下：
可以看到此时频谱除了在一些点上仍然一致但大部分不一样了。
我们知道在序列后面补零产生的频谱只是增加了分辨率，本身信息并没有增加。
由于本例中补了和序列长度一样多的零，故在偶数点（图中是n=1，3，5……15）仍然和8点频谱值一致，两者在偶数点的值相等。
但在奇数点的值不一样。
可以认为是两个不同的序列对应的插值方式不一样。
五、*MatLab上机内容
1、略。
2、将信号xb（n）的长度N设为63，用MatLab中randn（1,N）函数产生一个噪声信号w（n），计算将这个噪声信号叠加到xb（n）上以后新信号y（n）的频谱，观察发生的变化并记录。
当信噪比为0dB时，
可以看到由于叠加了随机噪声淹没了原来的信号（时域上如此），频谱除了个别点，所有值几乎是一样的（由于信噪比过低）。
改善了信噪比之后（约为20dB），同上作图：
此时比低信噪比有了明显改善，不但时域形状相差不大，频谱也仍大致保持了原样。
在频谱分量几乎为零的地方出现较大的偏差，这是有白噪声均匀分布引起的。
可以想见，当信噪比更高时，会有进一步的改善。
3、在步骤2的基础上，改变参数α和f，观察在出现混淆现象和泄漏现象的时候有噪声的y（n）信号的频谱有什么变化，是否明显？
下面是一个有明显混淆和泄漏现象的例子：
（α=0.01和f=0.025）
可见由于信号本来就混淆和混叠，只要信噪比足够，效应并不明显。
六、实验总结
本次实验是应用FFT对序列进行频谱分析。
通过本次实验，加深了对FFT的理解，
熟悉了FFT算法及其程序的编写，并应用FFT对典型信号进行了频谱分析。
了解了应用FFT进行信号频谱分析过程中可能出现的问题。
FFT并不是一种新变换，它是DFT的快速算法。
它是利用DFT系数的对称性和周期性将较大的序列的FFT分解为若干个较小的FFT并不断持续进行的递归过程，使得复杂度大大降低。
分为按时间抽取和按频率抽取两种。
利用FFT可以方便地进行序列的频谱分析和卷积运算。
但是，FFT可能会产生误差，即混淆、泄漏和栅栏效应（实际上也是DFT的问题）。
混淆、泄漏可以通过加大抽样频率解决，栅栏效应可以在序列后面补零解决。
七、思考题
1、实验中的信号序列号
和
，在单位圆上的z变换频谱
和）
会相同吗？
如果不同，你能说出哪一个低频分量更多一些吗？
为什么？
答：
对于连续的
，
、
不一样。
但是，由DFT的相移性质知在2πi/M离散点处，它们的幅频是一样的。
三角序列的低频分量更多些。
如果给序列前后补零（并不影响dtft曲线），观察序列形状，反三角序列的变化更为剧烈，又从0到N的跳变。
下面是两者dtft的比较（N=5），可以印证这一点：
三角序列的DTFT：
反三角序列的DTFT
2、对一个有限长序列进行离散傅里叶变换（DFT），等价于将该序列周期延拓后进行傅里叶级数（DFS）展开。
因为DFS也只是取其中一个周期来运算，所以FFT在一定条件下也可以用以分析周期信号序列。
如果实正弦信号sin（2πfn），f=0.1，用16点的FFT来做DFS运算，得到的频谱是信号本身的真实谱吗？
答：
不是。
给出的正弦信号周期是1/f=10,所以用FFT分析周期序列时，截取信号必须是整数个周期，否则会因为截断产生混淆和混叠。
八、源程序
主程序：
%产生高斯序列固定了p
fork=1:
3
p=pro_gauss（9,2^k）;subplot（3,1,k）;stem（abs（p））;if（k==1）title（'p=8,q分别等于2,4,8高斯序列的时域特性'）;
end
end
%高斯序列的幅频
fork=1:
3
p=pro_gauss（9,2^k）;subplot（3,1,k）;stem（abs（my_FFT（p,4）））;if（k==1）title（'p=8,q分别等于2,4,8高斯序列的时域特性'）;end
end
%高斯序列的固定了q
temp=[81314];
fork=1:
3
p=pro_gauss（temp（k）+1,8）;subplot（3,1,k）;stem（abs（p））;if（k==1）title（'q=8,p分别等于8,13,14高斯序列的时域特性'）;
end
end
%高斯序列的幅频
fork=1:
3
p=pro_gauss（temp（k）+1,8）;subplot（3,1,k）;stem（abs（my_FFT（p,4）））;if（k==1）title（'q=8,p分别等于8,13,14高斯序列的频域特性'）;
end
end
%产生衰减正弦信号
%
（1）
p=pro_desin（.1,0.0625）;subplot（2,1,1）;stem（p）;title（'α=0.1，f=0.0625时的衰减正弦序列'）
subplot（2,1,2）;
stem（abs（my_FFT（p,4）））;
title（'α=0.1，f=0.0625时的衰减正弦序列的幅频'）;
%
（2）
temp=[.4375,.5625]
forn=1:
2
p=pro_desin（.5,temp（n））;subplot（3,1,n）;stem（p）;
end
forn=1:
2
p=pro_desin（.5,temp（n））;subplot（3,1,n）;stem（abs（my_FFT（p,50）））;
end
p=pro_tringle（8）;
subplot（2,1,1）;stem（p）;title（'三角序列的时域及幅度频谱'）;
x=my_FFT（p,3）;subplot（2,1,2）;
stem（abs（x））;
stem（abs（x））;
p=pro_transtr（8）;
subplot（2,1,1）;stem（p）;title（'反三角序列的时域及幅度频谱'）;
x=my_FFT（p,3）;subplot（2,1,2）;
stem（abs（x））;
p=pro_tringle（16）;
x=my_FFT（p,4）;subplot（2,1,1）;title（'三角序列补零后的幅度频谱'）;
stem（abs（x））;p=pro_transtr（16）;
x=my_FFT（p,4）;subplot（2,1,2）;title（'反三角序列补零后的幅度频谱'）;
stem（abs（x））;
%上机实验
m=pro_desin（.1,.125）;subplot（2,1,1）;stem（m）;title（'一个衰减正弦序列'）
p=0.1*rand（1,64）;y=p+m;subplot（2,1,2）;stem（y）;title（'叠加了随机噪声的衰减正弦序列'）
subplot（2,1,1）;stem（abs（my_FFT（m,6）））;title（'一个衰减正弦序列的频谱'）
subplot（2,1,2）;stem（abs（my_FFT（y,6）））;title（'叠加了随机噪声的衰减正弦序列的频谱'）
子程序：
function[xb]=pro_desin（alpha,f）;
%产生衰减正弦序列
xb=zeros（1,64）;
forn=1:
64;
xb（n）=exp（-alpha*（n-1））.*sin（2*pi*（n-1）*f）;
end
function[gauss]=pro_gauss（p,q）
%产生高斯序列
gauss=zeros（1,16）;
forn=1:
16
gauss（n）=exp（-（n-p）^2/q）;
end
function[xb]=pro_desin（alpha,f）;
%产生衰减正弦序列
xb=zeros（1,64）;
forn=1:
64;
xb（n）=exp（-alpha*（n-1））.*sin（2*pi*（n-1）*f）;
end
function[xd]=pro_transtr（N）
%产生反三角序列
xd=zeros（1,N）;
forn=1:
8
ifn<5
xd（n）=5-n;
else
xd（n）=n-4;
end
end
function[xc]=pro_tringle（N）
%产生三角序列，长度为N¨
xc=zeros（1,N）;
forn=1:
8
ifn<5
xc（n）=n;
else
xc（n）=9-n;
end
end
function[X]=my_FFT（b,N）
%基2按时间抽取乱序输入顺序输出算法
%输入参数为b（n）,序列长度为2^N-1
%输出为序列X（k）.
%%%求反序，实现反序输入，正序输出%%%%
%采用递推的方法其中c（n）+1表示反序输入的第n个值的序列号%
c=zeros（1,2^N）;%c表示乱序和顺序的对应关系
c
（1）=0;c
（2）=2^（N-1）;
fori=1:
N-1
c（（2^i+1）:
2^（i+1））=c（1:
2^i）+2^（N-i-1）;
end
X=zeros（1,2^N）;
d=zeros（1,2^N）;
fori=1:
2^N
d（i）=b（c（i）+1）;
end%完成初始值的反序%
%%%求出所有的蝶算指数因子%%%%
%a是蝶形因子
x=rand（2^N,N+1）;a=zeros（1,2^（N-1））;q=1;w=zeros（1,2^（N-1））;
foro=1:
2^（N-1）
w（o）=exp（-j*2*pi*（o-1）/2^N）;
end
x（:
1）=d';
fork=2:
N+1
i=1;p=1;
form=1:
2^（k-2）%求出每级所需的蝶算指数因子%
q=（m-1）*2^N/（2^（k-1））;
a（p）=w（q+1）;p=p+1;
end
forn=1:
（2^N/（2^（k-1）））%对偶节点个数的确定及变换计算%
p=1;
forh=i:
（i+（2^（k-2）-1））
x（h,k）=x（h,k-1）+a（p）*x（h+2^（k-2）,k-1）;
x（h+2^（k-2）,k）=x（h,k-1）-a（p）*x（h+2^（k-2）,k-1）;
p=p+1;
end
i=（2^（k-1））*n+1;
end
end
X=（x（:
N+1））';
fori=1:
2^N
X（i）=conj（X（i））;
end
升级会员 