实验三 快速傅里叶变换及其应用.docx
《实验三 快速傅里叶变换及其应用.docx》由会员分享，可在线阅读，更多相关《实验三 快速傅里叶变换及其应用.docx（21页珍藏版）》请在冰豆网上搜索。
实验三快速傅里叶变换及其应用
实验三快速傅里叶变换及其应用
一:
实验目的
（1）加深对FFT的理解，熟悉matlab中的有关函数。
（2）应用FFT对典型信号进行频谱分析。
（3）了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT.
（4）应用FFT实现序列的线性卷积和相关。
二:
实验原理:
在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换（DFT）。
这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x（n）的长度为N时，
它的DFT定义为：
反变换为：
有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。
FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。
它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。
常用的FFT是以2为基数的，其长度
。
它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。
（一）在运用DFT进行频谱分析的过程中可能的产生三种误差
（1） 混叠
序列的频谱是被采样信号的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。
避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。
（2） 泄漏
实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。
泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。
为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。
（3） 栅栏效应
DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看，用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。
减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。
（二）用FFT计算线性卷积
用FFT可以实现两个序列的圆周卷积。
在一定的条件下，可以使圆周卷积等于线性卷积。
一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT的长度N≥N1＋N2
对于长度不足N的两个序列，分别将他们补零延长到N。
当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。
有两种方法：
（1）重叠相加法。
将长序列分成与短序列相仿的片段，分别用FFT对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。
（2） 重叠保留法。
这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输出时只需将各段线性卷积部分直接连接起来，省掉了输出段的直接相加。
（三）用FFT计算相关函数
两个长为
的实离散时间序列
与
的互相关函数定义为：
的离散傅里叶变换为：
当
时，得到
的自相关函数为：
利用FFT求两个有限长序列线性相关的步骤（设
长
，
长
）：
（1）为了使两个有限长序列的线性相关可用其圆周相关代替而不产生混淆，选择周期
≥
，以便使用FFT，将
，
补零至长为
。
（2）用FFT计算
（3）
（4）对
作IFFT；取后
项，得
；取前
项，得
。
三、实验内容及步骤
实验中用到的信号序列：
a）高斯（Gaussian）序列b）衰减正弦序列
c）三角波序列d）反三角波序列
上机实验内容：
（1）观察高斯序列的时域和幅频特性，固定信号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=fft（x）;
plot（abs（G（1:
15）））;xlabel（'k'）;
运行结果如下：
对比三张图可知：
P不变，随着q值的增大，时域信号幅值变化缓慢时域幅度对应变大,波形变胖，低频分量变多，频域信号频谱泄露程度减小。
参数q=8不变
代码和之前相比只要改变相应几个参数而已，故不列出；
运行结果如下:
可见，当q不变，随着p的增大，时域信号幅值不变，会在时间轴移位，对应右移,可见p决定了波形位置，
当实验中q=8,p=13时，明显出现泄漏。
（2）观察衰减正弦序列xb（n）的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？
说明产生现象的原因。
其中一段代码如下（其他的都是改变参数而已）：
a=0.1;f=0.0625;
fori=1:
16
x（i）=exp（-a*（i-1））*sin（2*pi*f*（i-1））;
end
fori=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（n（1:
16）,abs（G（1:
16）））;xlabel（'k'）;
结果：
（3）观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc（n）和xd（n）的幅频特性，观察两者的序列形状和频谱曲线有什么异同？
绘出两序列及其幅频特性曲线。
主要代码如下：
fori=1:
4
x（i）=i-1;
end
fori=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）））;xlabel（'k'）;
结果：
在xc（n）和xd（n）末尾补零，用N=32点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？
两情况的FFT频谱还有相同之处吗？
这些变化说明了什么？
代码只需改动几处数据即可，其中一段如下所示：
fori=1:
4
x（i）=i-1;
end
fori=5:
8
x（i）=9-i;
end
fori=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'）;
运行结果：
变化：
反三角波的低频分量增多，对信号末尾补零加长整数个周期可以对原信号达到细化频谱的作用。
（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不变，其结果有何不同，为什么？
其中一段代码如下：
实现方法和之前一样
N=16;
f=1/16;
forn=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=fft（x,16）;
plot（n（1:
16）,abs（G（1:
16）））;xlabel（'k'）;
运行结果：
（5）用FFT分别实现xa（n）（p＝8，q＝2）和 xb（n）（a＝0.1，f＝0.0625）的16点循环卷积和线性卷积。
fori=1:
16
x（i）=exp（（-（i-1-8）.^2）/2）;
y（i）=exp（-0.1*（i-1））*sin（2*pi*0.0625*（i-1））;
end
%fori=17:
31
%x（i）=0;
%y（i）=0;
%end
n=0:
30
G1=fft（x,31）;
G2=fft（y,31）;
z=ifft（G1.*G2,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段，分别采用重叠相加法和重叠保留法。
xe=rand（1,512）;
n1=0:
1:
3;
xc1=n1;
n2=4:
7;
xc2=8-n2;
xc=[xc1,xc2];
yn=zeros（1,519）;
forj=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）;
forj=0:
7
xj_1=xe_1（64*j+1:
64*j+71）;
xak_1=fft（xj_1）;
xck_1=fft（xc,71）;
yn1_1=ifft（xak_1.*xck_1）;
temp_1=zeros（1,519）;
temp_1（64*j+1:
64*j+64）=yn1_1（8:
71）;
yn_1=yn_1+temp_1;
end;
n=0:
518;
figure
（2）
subplot（2,1,1）;
plot（n,yn_1）;
xlabel（'n'）;ylabel（'y（n）'）;
subplot（2,1,2）;
plot（n,abs（fft（yn_1）））;
xlabel（'k'）;ylabel（'Y（k）'）;
axis（[0,600,0,300]）;
运行结果：
（7）用FFT分别计算xa（n）（p=8,q=2）和 xb（n）（a＝0.1，f＝0.0625）的16点循环相关和线性相关，问一共有多少种结果，它们之间有何共同点？
n=0:
1:
15;
xan=exp（-（n-8）.^2/2）;
xbn=exp（-0.1*n）.*sin（2*pi*0.0625*n）;
k=length（xbn）;
xan1=[xanzeros（1,k-1）];
xbn1=[xbnzeros（1,k-1）];
xak=fft（xan1）;
xbk=fft（xbn1）;
rm=real（ifft（conj（xak）.*xbk））;
rm1=[rm（k+1:
2*k-1）rm（1:
k）];
m=（-k+1）:
（k-1）;
subplot（2,1,1）;
stem（m,rm1）;xlabel（'n'）;
subplot（2,1,2）;
xak=fft（xan）;
xbk=fft（xbn）;
rm=real（ifft（conj（xak）.*xbk））;
stem（n,rm）;xlabel（'n'）;
运行结果：
（8）用FFT分别计算xa（n）（p=8,q=2）和 xb（n）（a＝0.1，f＝0.0625）的自相关函数。
n=0:
1:
15;
xan=exp（-（n-8）.^2/2）;
k=length（xan）;
xak=fft（xan,2*k）;
rm=real（ifft（conj（xak）.*xak））;
rm=[rm（k+2:
2*k）rm（1:
k）];
m=（-k+1）:
（k-1）;
subplot（2,1,1）;
stem（m,rm）;xlabel（'m'）;
运行结果：
四，思考题
（1）实验中的信号序列xc（n）和xd（n），在单位圆上的Z变换频谱|Xc（jω）|和|Xd（jω）|会相同吗？
如果不同，你能说出哪一个低频分量更多一些吗？
为什么？
答：
不同。
反三角波的低频分量更多一些。
（2）对一个有限长序列进行DFT等价于将该序列周期延拓后进行DFS展开，因为DFS也只是取其中一个周期来计算，所以FFT在一定条件下也可以用以分析周期信号序列。
如果实正弦信号sin（2πfn）,f=0.1用16点FFT来做DFS运算，得到的频谱时信号本身的真实谱吗？
升级会员 