实验三离散傅立叶变换PPT课件.ppt

关于实验三离散傅立叶变换第一张，PPT共四十六页，创作于2022年6月一、实验目的一、实验目的 加深对离散傅立叶变换加深对离散傅立叶变换(DFT)(DFT)的理解。的理解。掌握利用掌握利用MATLABMATLAB语言进行离散傅立叶变换和逆变换的方法。语言进行离散傅立叶变换和逆变换的方法。加深对离散傅立叶变换基本性质的理解。加深对离散傅立叶变换基本性质的理解。掌握离散傅立叶变换快速算法的应用。掌握离散傅立叶变换快速算法的应用。第二张，PPT共四十六页，创作于2022年6月二、实验原理及方法二、实验原理及方法 建立以时间建立以时间t t为自变量的为自变量的“信号信号”与以频率与以频率f f为为自变量的自变量的“频率函数频率函数”(”(频谱频谱)之间的某种变换关系。之间的某种变换关系。所以所以“时间时间”或或“频率频率”取连续还是离散值取连续还是离散值,就形成就形成各种不同形式的傅里叶变换对。各种不同形式的傅里叶变换对。傅里叶变换傅里叶变换第三张，PPT共四十六页，创作于2022年6月四种不同傅里叶变换对四种不同傅里叶变换对傅里叶级数傅里叶级数(FS)FS)：连续时间连续时间,离散频率的傅里叶变换。周期离散频率的傅里叶变换。周期连续时间信号傅里叶级数连续时间信号傅里叶级数(FS)FS)得到非周期离散频谱密度函数。得到非周期离散频谱密度函数。傅里叶变换傅里叶变换(FT)FT)：连续时间连续时间,连续频率的傅里叶变换。非周期连续频率的傅里叶变换。非周期连续时间信号通过连续付里叶变换连续时间信号通过连续付里叶变换(FT)FT)得到非周期连续频谱密度得到非周期连续频谱密度函数。函数。离散时间的离散时间的傅里叶变换傅里叶变换(DTFT):DTFT):离散时间离散时间,连续频率的傅里连续频率的傅里叶变换。非周期离散的时间信号叶变换。非周期离散的时间信号(单位园上的单位园上的Z Z变换变换(DTFT)DTFT)得到周期性连续的频率函数。得到周期性连续的频率函数。离散傅里叶变换离散傅里叶变换(DFT):DFT):离散时间离散时间,离散频率的傅里叶变换。离散频率的傅里叶变换。第四张，PPT共四十六页，创作于2022年6月上面讨论的前三种傅里叶变换对上面讨论的前三种傅里叶变换对,都不适用在计算机上都不适用在计算机上运算运算,因为至少在一个域因为至少在一个域(时域或频域时域或频域)中中,函数是连函数是连续的。因为从数字计算角度我们感兴趣的是时域及频续的。因为从数字计算角度我们感兴趣的是时域及频域都是离散的情况域都是离散的情况,这就是第四种离散傅里叶变换。这就是第四种离散傅里叶变换。第五张，PPT共四十六页，创作于2022年6月离散傅里叶级数离散傅里叶级数(DFS)DFS)离散时间序列离散时间序列x(n)x(n)满足满足x(n)=x(n+rN)x(n)=x(n+rN)，称为离散周期序列，称为离散周期序列，其中其中N N为周期，为周期，x(n)x(n)为主值序列。为主值序列。由傅立叶分析知道周期函数可由复指数的线性组合叠加得由傅立叶分析知道周期函数可由复指数的线性组合叠加得到。其频率为基本频率的倍数。从离散时间傅立叶变换的到。其频率为基本频率的倍数。从离散时间傅立叶变换的频率周期性，我们知道谐波次数是有限的，其频率为频率周期性，我们知道谐波次数是有限的，其频率为周期序列可表示成：周期序列可表示成：第六张，PPT共四十六页，创作于2022年6月其中其中 叫做离散傅立叶级数系数，也称为周叫做离散傅立叶级数系数，也称为周期序列的频谱，可由下式表示期序列的频谱，可由下式表示注意注意 也是一个基本周期为也是一个基本周期为N的周期序列。的周期序列。上面两式称为周期序列的傅立叶级数变换对。上面两式称为周期序列的傅立叶级数变换对。令令 表示复指数，可以得到以下：表示复指数，可以得到以下：第七张，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月利用利用MATLABMATLAB实现傅立叶级数计算实现傅立叶级数计算编写函数实现编写函数实现DFSDFS计算计算function xk=dfs(xn,N)function xk=dfs(xn,N)n=0:1:N-1;%n n=0:1:N-1;%n的行向量的行向量 k=n;%k k=n;%k的行向量的行向量 WN=exp(-j*2*pi/N);%W WN=exp(-j*2*pi/N);%WN N因子因子 nk=n*k;%nk=n*k;%产生一个含产生一个含nknk值的值的N N乘乘N N维矩阵维矩阵 WNnk=WN.nk;%DFS WNnk=WN.nk;%DFS矩阵矩阵 xk=xn*WNnk;%DFS xk=xn*WNnk;%DFS系数行向量系数行向量第十张，PPT共四十六页，创作于2022年6月例：例：xn=0,1,2,3xn=0,1,2,3，N=4N=4xn=0,1,2,3;xn=0,1,2,3;N=4;N=4;xk=dfs(xn,N)xk=dfs(xn,N)第十一张，PPT共四十六页，创作于2022年6月逆运算逆运算IDFSIDFSfunction xn=idfs(xk,N)function xn=idfs(xk,N)n=0:1:N-1;n=0:1:N-1;k=n;k=n;WN=exp(-j*2*pi/N);WN=exp(-j*2*pi/N);nk=n*k;nk=n*k;WNnk=WN.(-nk);WNnk=WN.(-nk);xn=xn=（xk*WNnkxk*WNnk）/N;/N;第十二张，PPT共四十六页，创作于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=length(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共四十六页，创作于2022年6月从上图可以看出，信号序列的周期数越多，则频谱越是向几个从上图可以看出，信号序列的周期数越多，则频谱越是向几个频点集中，当信号周期数趋于无穷大时频点集中，当信号周期数趋于无穷大时，频谱转化为离散谱。，频谱转化为离散谱。第十六张，PPT共四十六页，创作于2022年6月离散傅立叶变换（离散傅立叶变换（DFT）有限长序列有限长序列x(n)表示为表示为x(n)是非周期序列，但可以理解为周期序列是非周期序列，但可以理解为周期序列 的主值的主值序列。由离散傅立叶级数序列。由离散傅立叶级数DFS和和IDFS引出有限长序列的离引出有限长序列的离散傅立叶正、逆变换关系式散傅立叶正、逆变换关系式第十七张，PPT共四十六页，创作于2022年6月有限长序列傅立叶变换定义式为：有限长序列傅立叶变换定义式为：比较正、逆变换的定义式可以看出，只要把比较正、逆变换的定义式可以看出，只要把DFTDFT公式中的系公式中的系数数 改为改为 ，并最后乘以，并最后乘以1/1/N N，那么，那么，DFTDFT的计算的计算程序就可以用来计算程序就可以用来计算IDFTIDFT。第十八张，PPT共四十六页，创作于2022年6月DFT与与DFS的关系的关系比较两者的变换对，可以看出两者的区别仅仅是将周期序列换比较两者的变换对，可以看出两者的区别仅仅是将周期序列换成了有限长序列。成了有限长序列。有限长序列有限长序列x(n)x(n)可以看作是周期序列可以看作是周期序列 的一个周期；反的一个周期；反之周期序列之周期序列 可以看作是有限长序列可以看作是有限长序列x(n)x(n)以以N N为周期的周为周期的周期延拓。期延拓。由于公式非常相似，在程序编写上也基本一致。由于公式非常相似，在程序编写上也基本一致。第十九张，PPT共四十六页，创作于2022年6月例：已知序列例：已知序列x(n)=0,1,2,3,4,5,6,7，求，求x(n)的的DFT和和IDFT，画出序列傅立叶变换的幅度和相位图，并，画出序列傅立叶变换的幅度和相位图，并将原图像与逆变换图像进行比较。将原图像与逆变换图像进行比较。N=8;N=8;xn=0:N-1;n=0:N-1;xn=0:N-1;n=0:N-1;xk=dft(xn,N);xk=dft(xn,N);x=idft(xk,N);x=idft(xk,N);subplot(2,2,1);stem(n,xn)subplot(2,2,1);stem(n,xn)subplot(2,2,2);stem(n,abs(x)subplot(2,2,2);stem(n,abs(x)subplot(2,2,3);stem(n,abs(xk)subplot(2,2,3);stem(n,abs(xk)subplot(2,2,4);stem(n,angle(xk)subplot(2,2,4);stem(n,angle(xk)第二十张，PPT共四十六页，创作于2022年6月第二十一张，PPT共四十六页，创作于2022年6月三、快速傅立叶变换三、快速傅立叶变换有限长序列通过离散傅里叶变换有限长序列通过离散傅里叶变换(DFT)DFT)将其频域离散化将其频域离散化成有限长序列成有限长序列.但其计算量太大但其计算量太大(与与N N的平方成正比）的平方成正比）,很难实时地处理问题很难实时地处理问题,因此引出了快速傅里叶变换因此引出了快速傅里叶变换(FFT)FFT)。FFTFFT并不是一种新的变换形式并不是一种新的变换形式,它只是它只是DFTDFT的一种快速算的一种快速算法法.并且根据对序列分解与选取方法的不同而产生了并且根据对序列分解与选取方法的不同而产生了FFTFFT的多种算法的多种算法.第二十二张，PPT共四十六页，创作于2022年6月DFTDFT的快速算法的快速算法FFTFFT是数字信号处理的基本方法和基本技术，是数字信号处理的基本方法和基本技术，是必须牢牢掌握的。是必须牢牢掌握的。时间抽选时间抽选FFTFFT算法的理论推导和流图详见数字信号处理教算法的理论推导和流图详见数字信号处理教材。该算法遵循两条准则：材。该算法遵循两条准则：（1 1）对时间奇偶分；（）对时间奇偶分；（2 2）对频率前后分。）对频率前后分。这种算法的流图特点是：这种算法的流图特点是：（1 1）基本运算单元都是蝶形）基本运算单元都是蝶形 任何一个长度为任何一个长度为N=2N=2M M的序列，总可通过的序列，总可通过M M次分解最后成为次分解最后成为2 2点的点的DFTDFT计算。如图所示：计算。如图所示：第二十三张，PPT共四十六页，创作于2022年6月W WN Nk k称为旋转因子称为旋转因子 计算方程如下：计算方程如下：X Xm+1m+1(p)=X(p)=Xm m(p)+W(p)+WN Nk kX Xm m(q)(q)X Xm+1m+1(q)=X(q)=Xm m(p)-W(p)-WN Nk kX Xm m(q)(q)第二十四张，PPT共四十六页，创作于2022年6月 （2 2）同址（原位）计算）同址（原位）计算 这是由蝶形运算带来的好处，每一级蝶形运算的结果这是由蝶形运算带来的好处，每一级蝶形运算的结果X Xm+1m+1(p)(p)无须另外存储，只要再存入无须另外存储，只要再存入X Xm m(p)(p)中即可，中即可，X Xm+1m+1(q)(q)亦然。这样将大大节省存储单元。亦然。这样将大大节省存储单元。（3 3）变址计算）变址计算 输入为输入为“混序混序”（码位倒置）排列，输出按自然序排（码位倒置）排列，输出按自然序排列，因而对输入要进行列，因而对输入要进行“变址变址”计算（即码位倒置计算）。计算（即码位倒置计算）。“变址变址”实际上是一种实际上是一种“整序整序”的行为，目的是保证的行为，目的是保证“同址同址”。第二十五张，PPT共四十六页，创作于2022年6月FFTFFT的应用的应用凡是利用付里叶变换来进行分析、综合、变换的地凡是利用付里叶变换来进行分析、综合、变换的地方，都可以利用方，都可以利用FFTFFT算法来减少其计算量。算法来减少其计算量。FFTFFT主要应用在主要应用在 1 1、快速卷积、快速卷积 2 2、快速相关、快速相关 3 3、频谱分析、频谱分析第二十六张，PPT共四十六页，创作于2022年6月快速傅立叶变换的快速傅立叶变换的MATLABMATLAB实现实现提供提供fftfft函数计算函数计算DFTDFT格式格式 X=fft(x)X=fft(x)X=fft(x,N)X=fft(x,N)如果如果x x的长度小于的长度小于N N，则在其后填零使其成为，则在其后填零使其成为N N点序列，反点序列，反之对之对x x进行截断，若省略变量进行截断，若省略变量N N，则，则DFTDFT的长度即为的长度即为x x的长度。的长度。如果如果N N为为2 2的幂，则得到高速的基的幂，则得到高速的基-2FFT-2FFT算法；若算法；若N N不是不是2 2的乘方，的乘方，则为较慢的混合算法。则为较慢的混合算法。如果如果x x是矩阵，则是矩阵，则X X是对矩阵的每一列向量作是对矩阵的每一列向量作FFTFFT。第二十七张，PPT共四十六页，创作于2022年6月快速傅立叶逆变换（快速傅立叶逆变换（IFFTIFFT）函数调用格式函数调用格式 y=ifft(x)y=ifft(x)y=ifft(x,N)y=ifft(x,N)当当N N小于小于x x长度时，对长度时，对x x进行截断，当进行截断，当N N大于大于x x长度时，长度时，对对x x进行补零。进行补零。第二十八张，PPT共四十六页，创作于2022年6月fftshiftfftshift函数函数功能：功能：对对fftfft的输出进行重新排列，将零频分量移到频谱的中心。的输出进行重新排列，将零频分量移到频谱的中心。调用格式调用格式 y=fftshift(x)y=fftshift(x)当当x x为向量时，为向量时，fftshift(x)fftshift(x)直接将直接将x x中左右两半交换而产生中左右两半交换而产生y y。当当x x为矩阵时，为矩阵时，fftshift(x)fftshift(x)直接将直接将x x中左右、上下进行交换而中左右、上下进行交换而产生产生y y。第二十九张，PPT共四十六页，创作于2022年6月由题目可得由题目可得x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)fs=100N=128/1024例：已知信号由例：已知信号由15Hz15Hz幅值幅值0.50.5的正弦信号和的正弦信号和40Hz40Hz幅值幅值2 2的正弦信号组成，数据采样频率为的正弦信号组成，数据采样频率为100Hz,100Hz,试绘制试绘制N=128N=128点点DFTDFT的幅频图。的幅频图。第三十张，PPT共四十六页，创作于2022年6月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月利用利用FFTFFT进行功率谱的噪声分析进行功率谱的噪声分析已知带有测量噪声信号已知带有测量噪声信号 其中其中f1=50Hz,f2=120Hz,为均值为零、方差为为均值为零、方差为1的的随机信号，采样频率为随机信号，采样频率为1000Hz，数据点数，数据点数N=512。试绘制信号的功率谱图。试绘制信号的功率谱图。第三十三张，PPT共四十六页，创作于2022年6月t=0:0.001:0.6;t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(1,length(t);y=x+2*randn(1,length(t);Y=fft(y,512);Y=fft(y,512);P=Y.*conj(Y)/512;%P=Y.*conj(Y)/512;%求功率求功率f=1000*(0:255)/512;f=1000*(0:255)/512;subplot(2,1,1);subplot(2,1,1);plot(y);plot(y);subplot(2,1,2);subplot(2,1,2);plot(f,P(1:256);plot(f,P(1:256);第三十四张，PPT共四十六页，创作于2022年6月第三十五张，PPT共四十六页，创作于2022年6月序列长度和序列长度和FFTFFT的长度对信号频谱的影响。的长度对信号频谱的影响。已知信号已知信号 其中其中f1=15Hz,f2=40Hz,f1=15Hz,f2=40Hz,采样频率为采样频率为100100Hz.Hz.在下列情况下绘制其幅频谱。在下列情况下绘制其幅频谱。Ndata=32,Nfft=32;Ndata=32,Nfft=32;Ndata=32,Nfft=128;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年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月线性卷积的线性卷积的FFTFFT算法算法在在MATLABMATLAB实现卷积的函数为实现卷积的函数为CONVCONV，对于，对于N N值较小的向量，这值较小的向量，这是十分有效的。对于是十分有效的。对于N N值较大的向量卷积可用值较大的向量卷积可用FFTFFT加快计算速加快计算速度。度。由由DFTDFT性质可知，若性质可知，若DFTxDFTx1 1(n)=X(n)=X1 1(k),DFTx(k),DFTx2 2(n)=X(n)=X2 2(k)(k)则则 若若DFTDFT和和IDFTIDFT均采用均采用FFTFFT和和IFFTIFFT算法，可提高卷积速度。算法，可提高卷积速度。第四十张，PPT共四十六页，创作于2022年6月计算计算x1(n)x1(n)和和x2(n)x2(n)的线性卷积的的线性卷积的FFTFFT算法可由算法可由下面步骤实现下面步骤实现计算计算X1(k)=FFTx1(n);计算计算X2(k)=FFTx2(n);计算计算Y(k)=X1(k)X2(k);计算计算x1(n)*x2(n)=IFFTY(k).第四十一张，PPT共四十六页，创作于2022年6月用函数用函数convconv和和FFTFFT计算同一序列的卷积，比较其计算同一序列的卷积，比较其计算时间。计算时间。clockclock函数读取瞬时时钟函数读取瞬时时钟etime(t1,t2)etime(t1,t2)函数计算时刻函数计算时刻t1,t2t1,t2间所经历的时间。间所经历的时间。第四十二张，PPT共四十六页，创作于2022年6月四、练习题四、练习题验证线性性质验证线性性质如果有两个有限长序列分别为如果有两个有限长序列分别为x x1 1(n)(n)和和x x2 2(n)(n)，长度分别为，长度分别为N N1 1和和N N2 2，且，且y(n)=axy(n)=ax1 1(n)+bx(n)+bx2 2(n)(n)，(a(a，b b均为常数均为常数)则该则该y(n)y(n)的的N N点点DFTDFT为为 Y(k)=DFTy(n)=aXY(k)=DFTy(n)=aX1 1(k)+bX(k)+bX2 2(k)(0=k=N-1)(k)(0=k=N-1)其中：其中：N=max(NN=max(N1 1,N,N2 2)，X X1 1(k)(k)和和X X2 2(k)(k)分别为分别为x x1 1(n)(n)和和x x2 2(n)(n)的的N N点点DFTDFT。例：已知序列：例：已知序列：x x1 1(n)=0,1,2,4(n)=0,1,2,4 x x2 2(n)=1,0,1,0,1(n)=1,0,1,0,1第四十三张，PPT共四十六页，创作于2022年6月调试例题程序，掌握基本指令、函数的使用。调试例题程序，掌握基本指令、函数的使用。针对快速算法，比较下列几种情况针对快速算法，比较下列几种情况(设置不同序列长度设置不同序列长度N)N)DFT DFT和和FFTFFT运算速度运算速度 FFTFFT和卷积运算速度和卷积运算速度第四十四张，PPT共四十六页，创作于2022年6月五、实验报告要求五、实验报告要求 简述实验目的和实验原理。简述实验目的和实验原理。列写练习题的代码并绘制程序产生的图形。列写练习题的代码并绘制程序产生的图形。总结实验中你的收获和体会。总结实验中你的收获和体会。第四十五张，PPT共四十六页，创作于2022年6月感感谢谢大大家家观观看看第四十六张，PPT共四十六页，创作于2022年6月