实验三离散傅立叶变换

实验三离散傅立叶变换第一页，共四十三页，2022年，8月28日一、实验目的

加深对离散傅立叶变换(DFT)的理解。掌握利用MATLAB语言进行离散傅立叶变换和逆变换的方法。加深对离散傅立叶变换基本性质的理解。掌握离散傅立叶变换快速算法的应用。第二页，共四十三页，2022年，8月28日二、实验原理及方法

建立以时间t为自变量的“信号”与以频率f为自变量的“频率函数”(频谱)之间的某种变换关系。所以“时间”或“频率”取连续还是离散值,就形成各种不同形式的傅里叶变换对。傅里叶变换第三页，共四十三页，2022年，8月28日四种不同傅里叶变换对傅里叶级数(FS)：连续时间,离散频率的傅里叶变换。周期连续时间信号傅里叶级数(FS)得到非周期离散频谱密度函数。傅里叶变换(FT)：连续时间,连续频率的傅里叶变换。非周期连续时间信号通过连续付里叶变换(FT)得到非周期连续频谱密度函数。离散时间的傅里叶变换(DTFT):离散时间,连续频率的傅里叶变换。非周期离散的时间信号(单位园上的Z变换(DTFT))得到周期性连续的频率函数。离散傅里叶变换(DFT):离散时间,离散频率的傅里叶变换。第四页，共四十三页，2022年，8月28日上面讨论的前三种傅里叶变换对,都不适用在计算机上运算,因为至少在一个域(时域或频域)中,函数是连续的。因为从数字计算角度我们感兴趣的是时域及频域都是离散的情况,这就是第四种离散傅里叶变换。第五页，共四十三页，2022年，8月28日离散傅里叶级数(DFS)离散时间序列x(n)满足x(n)=x(n+rN)，称为离散周期序列，其中N为周期，x(n)为主值序列。由傅立叶分析知道周期函数可由复指数的线性组合叠加得到。其频率为基本频率的倍数。从离散时间傅立叶变换的频率周期性，我们知道谐波次数是有限的，其频率为周期序列可表示成：第六页，共四十三页，2022年，8月28日其中叫做离散傅立叶级数系数，也称为周期序列的频谱，可由下式表示注意也是一个基本周期为N的周期序列。上面两式称为周期序列的傅立叶级数变换对。令表示复指数，可以得到以下：第七页，共四十三页，2022年，8月28日例：求出下面周期序列的DFSx(n)={……,0,1,2,3,0,1,2,3,0,1,2,3,……}基本周期为N=4，WN=W4=-j，

因而第八页，共四十三页，2022年，8月28日MATLAB实现矩阵-向量相乘运算来实现。由于和均为周期函数，周期为N，可设和代表序列和的主值区间序列，则前面的两个表达式可写成：式中，矩阵WN为方阵——DFS矩阵。第九页，共四十三页，2022年，8月28日利用MATLAB实现傅立叶级数计算编写函数实现DFS计算functionxk=dfs(xn,N)n=[0:1:N-1];%n的行向量k=n;%k的行向量WN=exp(-j*2*pi/N);%WN因子nk=n’*k;%产生一个含nk值的N乘N维矩阵WNnk=WN.^nk;%DFS矩阵xk=xn*WNnk;%DFS系数行向量第十页，共四十三页，2022年，8月28日例：xn=[0,1,2,3]，N=4xn=[0,1,2,3];N=4;xk=dfs(xn,N)’第十一页，共四十三页，2022年，8月28日逆运算IDFSfunctionxn=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;第十二页，共四十三页，2022年，8月28日xn=idfs(xk',4)x=xn'第十三页，共四十三页，2022年，8月28日周期重复次数对序列频谱的影响理论上讲，周期序列不满足绝对可积条件，要对周期序列进行分析，可以先取K个周期进行处理，然后让K无限增大，研究其极限情况。这样可以观察信号序列由非周期到周期变换时，频谱由连续谱逐渐向离散谱过渡的过程。第十四页，共四十三页，2022年，8月28日例：已知一个矩形序列的脉冲宽度占整个周期的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));forr=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第十五页，共四十三页，2022年，8月28日从上图可以看出，信号序列的周期数越多，则频谱越是向几个频点集中，当信号周期数趋于无穷大时，频谱转化为离散谱。第十六页，共四十三页，2022年，8月28日离散傅立叶变换（DFT）有限长序列x(n)表示为x(n)是非周期序列，但可以理解为周期序列的主值序列。由离散傅立叶级数DFS和IDFS引出有限长序列的离散傅立叶正、逆变换关系式第十七页，共四十三页，2022年，8月28日有限长序列傅立叶变换定义式为：

比较正、逆变换的定义式可以看出，只要把DFT公式中的系数改为，并最后乘以1/N，那么，DFT的计算程序就可以用来计算IDFT。第十八页，共四十三页，2022年，8月28日DFT与DFS的关系比较两者的变换对，可以看出两者的区别仅仅是将周期序列换成了有限长序列。有限长序列x(n)可以看作是周期序列的一个周期；反之周期序列可以看作是有限长序列x(n)以N为周期的周期延拓。由于公式非常相似，在程序编写上也基本一致。第十九页，共四十三页，2022年，8月28日例：已知序列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,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))第二十页，共四十三页，2022年，8月28日第二十一页，共四十三页，2022年，8月28日三、快速傅立叶变换有限长序列通过离散傅里叶变换(DFT)将其频域离散化成有限长序列.但其计算量太大(与N的平方成正比）,很难实时地处理问题,因此引出了快速傅里叶变换(FFT)。FFT并不是一种新的变换形式,它只是DFT的一种快速算法.并且根据对序列分解与选取方法的不同而产生了FFT的多种算法.第二十二页，共四十三页，2022年，8月28日DFT的快速算法—FFT是数字信号处理的基本方法和基本技术，是必须牢牢掌握的。时间抽选FFT算法的理论推导和流图详见《数字信号处理》教材。该算法遵循两条准则：（1）对时间奇偶分；（2）对频率前后分。这种算法的流图特点是：（1）基本运算单元都是蝶形

任何一个长度为N=2M的序列，总可通过M次分解最后成为2点的DFT计算。如图所示：第二十三页，共四十三页，2022年，8月28日WNk称为旋转因子计算方程如下：Xm+1(p)=Xm(p)+WNkXm(q)Xm+1(q)=Xm(p)-WNkXm(q)第二十四页，共四十三页，2022年，8月28日（2）同址（原位）计算这是由蝶形运算带来的好处，每一级蝶形运算的结果Xm+1(p)无须另外存储，只要再存入Xm(p)中即可，Xm+1(q)亦然。这样将大大节省存储单元。（3）变址计算输入为“混序”（码位倒置）排列，输出按自然序排列，因而对输入要进行“变址”计算（即码位倒置计算）。“变址”实际上是一种“整序”的行为，目的是保证“同址”。第二十五页，共四十三页，2022年，8月28日FFT的应用凡是利用付里叶变换来进行分析、综合、变换的地方，都可以利用FFT算法来减少其计算量。FFT主要应用在1、快速卷积2、快速相关3、频谱分析第二十六页，共四十三页，2022年，8月28日快速傅立叶变换的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。第二十七页，共四十三页，2022年，8月28日快速傅立叶逆变换（IFFT）函数调用格式

y=ifft(x)y=ifft(x,N)当N小于x长度时，对x进行截断，当N大于x长度时，对x进行补零。第二十八页，共四十三页，2022年，8月28日fftshift函数功能：对fft的输出进行重新排列，将零频分量移到频谱的中心。调用格式

y=fftshift(x)

当x为向量时，fftshift(x)直接将x中左右两半交换而产生y。当x为矩阵时，fftshift(x)直接将x中左右、上下进行交换而产生y。第二十九页，共四十三页，2022年，8月28日由题目可得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的幅频图。第三十页，共四十三页，2022年，8月28日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点’)第三十一页，共四十三页，2022年，8月28日第三十二页，共四十三页，2022年，8月28日利用FFT进行功率谱的噪声分析已知带有测量噪声信号其中f1=50Hz,f2=120Hz,为均值为零、方差为1的随机信号，采样频率为1000Hz，数据点数N=512。试绘制信号的功率谱图。第三十三页，共四十三页，2022年，8月28日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));第三十四页，共四十三页，2022年，8月28日第三十五页，共四十三页，2022年，8月28日序列长度和FFT的长度对信号频谱的影响。已知信号其中f1=15Hz,f2=40Hz,采样频率为100Hz.

在下列情况下绘制其幅频谱。

Ndata=32,Nfft=32;Ndata=32,Nfft=128;第三十六页，共四十三页，2022年，8月28日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’)第三十七页，共四十三页，2022年，8月28日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’)第三十八页，共四十三页，2022年，8月28日第三十九页，共四十三页，2022年，8月28日线性卷积的FFT算法在MATLAB实现卷积的函数为CONV，对于N值较小的向量，这是十分有效的。对于N值较大的向量卷积可用FFT加快计算速度。由DFT性质可知，若DFT[x1(n)]=X1