《吴镇杨matlab实验三快速傅里叶变换及其应用》由会员分享，可在线阅读，更多相关《吴镇杨matlab实验三快速傅里叶变换及其应用（14页珍藏版）》请在金锄头文库上搜索。
1、DSP试验04008012 1 实验三快速傅里叶变换及其应用一 :实验目的（ 1） 加 深对 FFT的理解，熟悉matlab 中的有关函数。（ 2） 应 用 FFT对典型信号进行频谱分析。（ 3） 了 解应用 FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT. （ 4） 应 用 FFT实现序列的线性卷积和相关。二 :实验原理 : 在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier 变换 (DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N 时，它的 DFT定义为：反变换
2、为：有限长序列的DFT是其 Z 变换在单位圆上的等距采样，或者说是序列Fourier 变换的等距采样，因此可以用于序列的谱分析。FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。常用的FFT是以 2 为基数的，其长度。它的效率高，程序简单，使用非常方便， 当要变换的序列长度不等于2 的整数次方时， 为了使用以2 为基数的FFT ，可以用末位补零的方法，使其长度延长至2 的整数次方。(一)在运用 DFT进行频谱分析的过程中可能的产生三种误差(1) 混叠序列的频谱是被采样信号的周期延拓，当采样速率不
3、满足Nyquist 定理时， 就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。(2) 泄漏实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT 来对信号进行频谱分析， 这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。泄漏不能与混叠完全分开，因为泄漏导致频谱的扩
4、展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。(3) 栅栏效应DFT是对单位圆上Z 变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看， 用 DFT来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。(二)用 FFT计算线性卷积用 FF
5、T可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性卷积。一般情况，设两个序列的DSP试验04008012 2 长度分别为N1 和 N2，要使圆周卷积等于线性卷积的充要条件是FFT的长度 NN1N2 对于长度不足N 的两个序列，分别将他们补零延长到N。当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。有两种方法：（1）重叠相加法。将长序列分成与短序列相仿的片段，分别用FFT对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。（2） 重叠保留法。 这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输出时只需将各段线性卷积部分直接连接
6、起来，省掉了输出段的直接相加。（三）用FFT计算相关函数两个长为N的实离散时间序列)(nx与)(ny的互相关函数定义为：)()()()()()()(1010mymxnymnxmnynxmrNnNmxy)( nrxy的离散傅里叶变换为：10),()()(NkkYkXkRxy当)()(nynx时，得到)(nx的自相关函数为：km NjNnNkxxekX Nmnxnxmr221010)(1 )()()(利用 FFT求两个有限长序列线性相关的步骤（设)(nx长1N，)( ny长2N） ：（1）为了使两个有限长序列的线性相关可用其圆周相关代替而不产生混淆，选择周期LN2121NN，以便使用 FFT ，将
7、)(nx，)(ny补零至长为N。（2）用 FFT计算10),(),(NkkYkX（3）10),()()(NkkYkXkRxy（4）对)( kRxy作 IFFT ；取后1N项，得11)(mNmrxy；取前N项，得10)(Nmmrxy。三、实验内容及步骤实验中用到的信号序列：a)高斯（ Gaussian ）序列b)衰减正弦序列其它0150)(2)(nenxqpna 其它0150)2sin()(nfnenxanbc)三角波序列d)反三角波序列DSP试验04008012 3 其它074830)(nnnnnxc其它0744304)(nnnnnxd上机实验内容：（1）观察高斯序列的时域和幅频特性，固定信号
8、xa(n)中参数 p=8，改变 q 的值，使q 分别等于2，4，8，观察它们的时域和幅频特性，了解当q 取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变 p，使 p 分别等于 8，13，14，观察参数p 变化对信号序列的时域及幅频特性的影响，观察p 等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。参数 p=8 不变主要代码如下 : i=1:15; p=8;q=2; subplot(3,2,1); x(i)=exp(-(i-p).2)/q); stem(x);xlabel(n); subplot(3,2,2); G=ff
9、t(x); plot(abs(G(1:15);xlabel(k); 运行结果如下：对比三张图可知：P 不变，随着q 值的增大，时域信号幅值变化缓慢时域幅度对应变大,波形变胖，低频分量变多，频域信号频谱泄露程DSP试验04008012 4 度减小。参数 q=8 不变代码和之前相比只要改变相应几个参数而已，故不列出；运行结果如下 : 可见，当q 不变，随着p 的增大，时域信号幅值不变，会在时间轴移位，对应右移,可见 p 决定了波形位置，当实验中q=8,p=13 时，明显出现泄漏。(2)观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状
10、，绘出幅频特性曲线，改变f，使 f 分别等于0.4375 和 0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。其中一段代码如下（其他的都是改变参数而已）：a=0.1;f=0.0625; DSP试验04008012 5 for i=1:16 x(i)=exp(-a*(i-1)*sin(2*pi*f*(i-1); end for i=17:100 x(i)=0; end n=0:15; subplot(3,2,1); plot(n(1:16),x(1:16);xlabel(n); subplot(3,2,2); G=fft(x,16); plot(
11、n(1:16),abs(G(1:16);xlabel(k); 结果：(3)观察三角波和反三角波序列的时域和幅频特性，用 N=8 点 FFT分析信号序列xc(n)和 xd(n)的幅频特性， 观察两者的序DSP试验04008012 6 列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。主要代码如下：for i=1:4 x(i)=i-1; end for i=5:8 x(i)=9-i; end n=0:7; subplot(2,2,1); plot(n,x(1:8);xlabel(n); subplot(2,2,2); G=fft(x,8); plot(n(1:8),abs(G(1:8);xl
12、abel(k); 结果：在 xc(n)和 xd(n)末尾补零，用N=32 点 FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFTDSP试验04008012 7 频谱还有相同之处吗？这些变化说明了什么？代码只需改动几处数据即可，其中一段如下所示：for i=1:4 x(i)=i-1; end for i=5:8 x(i)=9-i; end for i=9:32 x(i)=0; end n=0:31; subplot(2,1,1); G=fft(x,32); plot(n(1:32),abs(G(1:32);xlabel(k); 运行结果：变化：反三角波的低频分量增多，对信
13、号末尾补零加长整数个周期可以对原信号达到细化频谱的作用。(4) 一个连续信号含两个频率分量，经采样得x(n)=sin*2 *0.125n+cos*2 *(0.125+ f)n+n=0,1 ,N-1 已知 N=16，f 分别为 1 /16 和 1/64，观察其频谱；当N=128 时， f 不变，其结果有何不同，为什么？DSP试验04008012 8 其中一段代码如下：实现方法和之前一样N=16; f=1/16; for n=1:N x(n)=sin(2*pi*0.125*(n-1)+cos(2*pi*(0.125+f)*(n-1); end n=0:15; subplot(2,2,1); G=f
14、ft(x,16); plot(n(1:16),abs(G(1:16);xlabel(k); 运行结果：（5）用 FFT分别实现xa(n)（p8， q2）和 xb(n)（a0.1，f0.0625）的 16 点循环卷积和线性卷积。for i=1:16 DSP试验04008012 9 x(i)=exp(-(i-1-8).2)/ 2); y(i)=exp(-0.1*(i-1)*sin(2*pi*0.0625*(i-1); end %for i=17:31 % x(i)=0; % y(i)=0; %end n=0:30 G1=fft(x,31); G2=fft(y,31); z=ifft(G1.*G2,
15、31); subplot(2,1,1); plot(n(1:31),z(1:31); subplot(2,1,2); G1=fft(x,16); G2=fft(y,16); z=ifft(G1.*G2,16); n=0:15; plot(n(1:16),z(1:16); 运行结果：（6）产生一512 点的随机序列xe(n),并用 xc(n)和 xe(n)作线性卷积，观察卷积前后xe(n)频谱的变化。要求将xe(n)分成8 段，分别采用重叠相加法和重叠保留法。DSP试验04008012 10 xe=rand(1,512); n1=0:1:3; xc1=n1; n2=4:7; xc2=8-n2;
16、xc=xc1,xc2; yn=zeros(1,519); for j=0:7 xj=xe(64*j+1:64*(j+1); xak=fft(xj,71); xck=fft(xc,71); yn1=ifft(xak.*xck); temp=zeros(1,519); temp(64*j+1:64*j+71)=yn1; yn=yn+temp; end; n=0:518; figure(1) subplot(2,1,1); plot(n,yn); xlabel(n);ylabel(y(n); subplot(2,1,2); plot(n,abs(fft(yn); xlabel(k);ylabel(Y(k); axis(0,600,0,300); xe=rand(1,512); k=1:7; xe1=k-k; xe_1=xe1,xe; yn_1=zeros(1,519); for j=0:7 xj_1=xe_1(64*j+1:
