实验三离散傅立叶变换

1、关于实验三离散傅立叶变换第一张，PPT共四十六页，创作于2022年6月一、实验目的 加深对离散傅立叶变换(DFT)的理解。 掌握利用MATLAB语言进行离散傅立叶变换和逆变换的方法。加深对离散傅立叶变换基本性质的理解。掌握离散傅立叶变换快速算法的应用。第二张，PPT共四十六页，创作于2022年6月二、实验原理及方法 建立以时间t为自变量的“信号”与以频率f为自变量的“频率函数”(频谱)之间的某种变换关系。所以“时间”或“频率”取连续还是离散值, 就形成各种不同形式的傅里叶变换对。傅里叶变换第三张，PPT共四十六页，创作于2022年6月四种不同傅里叶变换对傅里叶级数(FS)：连续时间, 离散频率

2、的傅里叶变换。周期连续时间信号傅里叶级数(FS)得到非周期离散频谱密度函数。傅里叶变换(FT)：连续时间, 连续频率的傅里叶变换。非周期连续时间信号通过连续付里叶变换(FT)得到非周期连续频谱密度函数。离散时间的傅里叶变换(DTFT):离散时间,连续频率的傅里叶变换。非周期离散的时间信号(单位园上的Z变换(DTFT)得到周期性连续的频率函数。离散傅里叶变换(DFT):离散时间, 离散频率的傅里叶变换。第四张，PPT共四十六页，创作于2022年6月上面讨论的前三种傅里叶变换对,都不适用在计算机上运算, 因为至少在一个域( 时域或频域)中, 函数是连续的。因为从数字计算角度我们感兴趣的是时域及频域

3、都是离散的情况, 这就是第四种离散傅里叶变换。第五张，PPT共四十六页，创作于2022年6月离散傅里叶级数(DFS)离散时间序列x(n)满足x(n)=x(n+rN)，称为离散周期序列，其中N为周期，x(n)为主值序列。由傅立叶分析知道周期函数可由复指数的线性组合叠加得到。其频率为基本频率的倍数。从离散时间傅立叶变换的频率周期性，我们知道谐波次数是有限的，其频率为周期序列可表示成：第六张，PPT共四十六页，创作于2022年6月其中 叫做离散傅立叶级数系数，也称为周期序列的频谱，可由下式表示注意 也是一个基本周期为N的周期序列。上面两式称为周期序列的傅立叶级数变换对。令 表示复指数，可以得到以下：

4、第七张，PPT共四十六页，创作于2022年6月例：求出下面周期序列的DFSx(n)=,0,1,2,3,0,1,2,3,0,1,2,3, 基本周期为N=4，WN=W4=-j， 因而第八张，PPT共四十六页，创作于2022年6月MATLAB实现矩阵-向量相乘运算来实现。由于 和 均为周期函数，周期为N，可设 和 代表序列 和 的主值区间序列，则前面的两个表达式可写成：式中，矩阵WN为方阵DFS矩阵。第九张，PPT共四十六页，创作于2022年6月利用MATLAB实现傅立叶级数计算编写函数实现DFS计算function xk=dfs(xn,N) n=0:1:N-1; %n的行向量 k=n; %k的行向

5、量 WN=exp(-j*2*pi/N); %WN因子 nk=n*k; %产生一个含nk值的N乘N维矩阵 WNnk=WN.nk; %DFS矩阵 xk=xn* WNnk; %DFS系数行向量第十张，PPT共四十六页，创作于2022年6月例：xn=0,1,2,3，N=4xn=0,1,2,3;N=4;xk=dfs(xn,N)第十一张，PPT共四十六页，创作于2022年6月逆运算IDFSfunction xn=idfs(xk,N) n=0:1:N-1; k=n; WN=exp(-j*2*pi/N); nk=n*k; WNnk=WN.(-nk); xn=（xk*WNnk）/N;第十二张，PPT共四十六页，

6、创作于2022年6月xn=idfs(xk,4)x=xn第十三张，PPT共四十六页，创作于2022年6月周期重复次数对序列频谱的影响理论上讲 ，周期序列不满足绝对可积条件，要对周期序列进行分析，可以先取K个周期进行处理，然后让K无限增大，研究其极限情况。这样可以观察信号序列由非周期到周期变换时 ，频谱由连续谱逐渐向离散谱过渡的过程。第十四张，PPT共四十六页，创作于2022年6月例：已知一个矩形序列的脉冲宽度占整个周期的1/2，一个周期的采样点数为10，用傅立叶级数变换求信号的重复周期数分别为1、4、7、10时的幅度频谱。MATLAB程序：xn=ones(1,5),zeros(1,5);Nx=l

7、ength(xn);Nw=1000;dw=2*pi/Nw;k=floor(-Nw/2+0.5):(Nw/2+0.5);for r=0:3; K=3*r+1; nx=0:(K*Nx-1); x=xn(mod(nx,Nx)+1); Xk=x*(exp(-j*dw*nx*k)/K; subplot(4,2,2*r+1); stem(nx,x) axis(0,K*Nx-1,0,1.1); ylabel(x(n); subplot(4,2,2*r+2); plot(k*dw,abs(Xk) axis(-4,4,0,1.1*max(abs(Xk); ylabel(X(k);end第十五张，PPT共四十六页

8、，创作于2022年6月从上图可以看出，信号序列的周期数越多，则频谱越是向几个频点集中，当信号周期数趋于无穷大时 ，频谱转化为离散谱。第十六张，PPT共四十六页，创作于2022年6月离散傅立叶变换（DFT）有限长序列x(n)表示为x(n)是非周期序列，但可以理解为周期序列 的主值序列。由离散傅立叶级数DFS和IDFS引出有限长序列的离散傅立叶正、逆变换关系式第十七张，PPT共四十六页，创作于2022年6月有限长序列傅立叶变换定义式为： 比较正、逆变换的定义式可以看出，只要把DFT公式中的系数 改为 ，并最后乘以1/N，那么，DFT的计算程序就可以用来计算IDFT。第十八张，PPT共四十六页，创作

9、于2022年6月DFT与DFS的关系比较两者的变换对，可以看出两者的区别仅仅是将周期序列换成了有限长序列。有限长序列x(n)可以看作是周期序列 的一个周期；反之周期序列 可以看作是有限长序列x(n)以N为周期的周期延拓。由于公式非常相似，在程序编写上也基本一致。第十九张，PPT共四十六页，创作于2022年6月例：已知序列x(n)=0,1,2,3,4,5,6,7，求x(n)的DFT和IDFT，画出序列傅立叶变换的幅度和相位图，并将原图像与逆变换图像进行比较。N=8;xn=0:N-1; n=0:N-1;xk=dft(xn,N);x=idft(xk,N);subplot(2,2,1); stem(n

10、,xn)subplot(2,2,2); stem(n,abs(x)subplot(2,2,3); stem(n,abs(xk)subplot(2,2,4); stem(n,angle(xk)第二十张，PPT共四十六页，创作于2022年6月第二十一张，PPT共四十六页，创作于2022年6月三、快速傅立叶变换有限长序列通过离散傅里叶变换(DFT)将其频域离散化成有限长序列.但其计算量太大(与N的平方成正比）, 很难实时地处理问题, 因此引出了快速傅里叶变换(FFT)。FFT并不是一种新的变换形式,它只是DFT的一种快速算法.并且根据对序列分解与选取方法的不同而产生了FFT的多种算法.第二十二张，P

11、PT共四十六页，创作于2022年6月DFT的快速算法FFT是数字信号处理的基本方法和基本技术，是必须牢牢掌握的。 时间抽选FFT算法的理论推导和流图详见数字信号处理教材。该算法遵循两条准则： （1）对时间奇偶分；（2）对频率前后分。 这种算法的流图特点是： （1）基本运算单元都是蝶形 任何一个长度为N=2M的序列，总可通过M次分解最后成为2点的DFT计算。如图所示： 第二十三张，PPT共四十六页，创作于2022年6月WNk称为旋转因子 计算方程如下： Xm+1(p)=Xm(p)+WNkXm(q) Xm+1(q)=Xm(p)-WNkXm(q)第二十四张，PPT共四十六页，创作于2022年6月 （

12、2）同址（原位）计算 这是由蝶形运算带来的好处，每一级蝶形运算的结果Xm+1(p)无须另外存储，只要再存入Xm(p)中即可，Xm+1(q)亦然。这样将大大节省存储单元。 （3）变址计算 输入为“混序”（码位倒置）排列，输出按自然序排列，因而对输入要进行“变址”计算（即码位倒置计算）。“变址”实际上是一种“整序”的行为，目的是保证“同址”。 第二十五张，PPT共四十六页，创作于2022年6月FFT的应用凡是利用付里叶变换来进行分析、综合、变换的地方，都可以利用FFT算法来减少其计算量。FFT主要应用在 1、快速卷积 2、快速相关 3、频谱分析第二十六张，PPT共四十六页，创作于2022年6月快速

13、傅立叶变换的MATLAB实现提供fft函数计算DFT格式 X=fft(x) X=fft(x,N) 如果x的长度小于N，则在其后填零使其成为N点序列，反之对x进行截断，若省略变量N，则DFT的长度即为x的长度。如果N为2的幂，则得到高速的基-2FFT算法；若N不是2的乘方，则为较慢的混合算法。如果x是矩阵，则X是对矩阵的每一列向量作FFT。第二十七张，PPT共四十六页，创作于2022年6月快速傅立叶逆变换（IFFT）函数调用格式 y=ifft(x) y=ifft(x,N)当N小于x长度时，对x进行截断，当N大于x长度时，对x进行补零。第二十八张，PPT共四十六页，创作于2022年6月fftshi

14、ft函数功能： 对fft的输出进行重新排列，将零频分量移到频谱的中心。调用格式 y=fftshift(x) 当x为向量时，fftshift(x)直接将x中左右两半交换而产生y。 当x为矩阵时，fftshift(x)直接将x中左右、上下进行交换而产生y。第二十九张，PPT共四十六页，创作于2022年6月由题目可得x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)fs=100N=128/1024例：已知信号由15Hz幅值0.5的正弦信号和40Hz幅值2的正弦信号组成，数据采样频率为100Hz,试绘制N=128点DFT的幅频图。第三十张，PPT共四十六页，创作于2022年6月

15、fs=100;N=128;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,N);f=(0:length(y)-1)*fs/length(y);mag=abs(y);stem(f,mag);title(N=128点)第三十一张，PPT共四十六页，创作于2022年6月第三十二张，PPT共四十六页，创作于2022年6月利用FFT进行功率谱的噪声分析已知带有测量噪声信号 其中f1=50Hz,f2=120Hz, 为均值为零、方差为1的随机信号，采样频率为1000Hz，数据点数N=512。试绘制信号的功率谱图。第三十三张，PPT共四

16、十六页，创作于2022年6月t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(1,length(t);Y=fft(y,512);P=Y.*conj(Y)/512; %求功率f=1000*(0:255)/512;subplot(2,1,1);plot(y);subplot(2,1,2);plot(f,P(1:256);第三十四张，PPT共四十六页，创作于2022年6月第三十五张，PPT共四十六页，创作于2022年6月序列长度和FFT的长度对信号频谱的影响。已知信号 其中f1=15Hz,f2=40Hz,采样频率为100Hz. 在下

17、列情况下绘制其幅频谱。 Ndata=32,Nfft=32; Ndata=32,Nfft=128;第三十六张，PPT共四十六页，创作于2022年6月fs=100;Ndata=32; Nfft=32;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(2,1,1)plot(f(1:Nfft/2),mag(1:Nfft/2)title(Ndata=32,Nfft=32)第三十七张，PPT共四十六页，创作于2022

18、年6月Nfft=128;n=0:Ndata-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);y=fft(x,Nfft);mag=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(2,1,2)plot(f(1:Nfft/2),mag(1:Nfft/2)title(Ndata=32,Nfft=128)第三十八张，PPT共四十六页，创作于2022年6月第三十九张，PPT共四十六页，创作于2022年6月线性卷积的FFT算法在MATLAB实现卷积的函数为CONV，对于N值较小的向量，这是十分有效的。对于N值较大的向量卷积可用FFT加快计算速度。由DFT性质可知，若DFTx1(n)=X1(k),DFTx2(n)=X2(k)则 若DFT和IDFT均采用FFT和IFFT算法，可提高卷积速度。第四十张，PPT共四十六页，创作于2022年6月计算x1(n)和x2(n)的线性卷积的FFT算法可由下面步骤实现计算X1(k)=FFTx1(n);计算X2(