用Matlab对信号进行傅里叶变换实例

1目录用Matlab对信号进行傅里叶变换.........................................................2Matlab的傅里叶变换实例.....................................................................5Matlab方波傅立叶变换画出频谱图......................................................72用Matlab对信号进行傅里叶变换1.离散序列的傅里叶变换DTFT(DiscreteTimeFourierTransform)代码：1N=8;%原离散信号有8点2n=[0:1:N-1]%原信号是1行8列的矩阵3xn=0.5.^n;%构建原始信号，为指数信号45w=[-800:1:800]*4*pi/800;%频域共-800----+800的长度（本应是无穷，高频分量很少，故省去）6X=xn*exp(-j*(n'*w));%求dtft变换，采用原始定义的方法，对复指数分量求和而得7subplot(311)8stem(n,xn);9title('原始信号(指数信号)');10subplot(312);11plot(w/pi,abs(X));12title('DTFT变换')结果：分析：可见，离散序列的dtft变换是周期的，这也符合Nyquist采样定理的描述，连续时间信号经周期采样之后，所得的离散信号的频谱是原连续信号频谱的周期延拓。2.离散傅里叶变换DFT(DiscreteFourierTransform)3与1中DTFT不一样的是，DTFT的求和区间是整个频域，这对计算机的计算来说是不可以实现的，DFT就是序列的有限傅里叶变换。实际上，1中代码也只是对频域的-800----+800中间的1601结果图：1234567891011121314151617N=8;%原离散信号有8点n=[0:1:N-1]%原信号是1行8列的矩阵xn=0.5.^n;%构建原始信号，为指数信号w=[-8:1:8]*4*pi/8;%频域共-800----+800的长度（本应是无穷，高频分量很少，故省去）X=xn*exp(-j*(n'*w));%求dtft变换，采用原始定义的方法，对复指数分量求和而得subplot(311)stem(n,xn);w1=[-4:1:4]*4*pi/4;X1=xn*exp(-j*(n'*w1));title('原始信号(指数信号)');subplot(312);stem(w/pi,abs(X));title('原信号的16点DFT变换')subplot(313)stem(w1/pi,abs(X1));title('原信号的8点DFT变换')4分析：DFT只是DTFT的现实版本，因为DTFT要求求和区间无穷，而DFT只在有限点内求和。3.快速傅里叶变换FFT（FastFourierTransform）虽然DFT相比DTFT缩减了很大的复杂度，但是任然有相当大的计算量，不利于信息的实时有效处理，1965年发现的DFT解决了这一问题。实现代码：1N=64;%原离散信号有8点2n=[0:1:N-1]%原信号是1行8列的矩阵3xn=0.5.^n;%构建原始信号，为指数信号4Xk=fft(xn,N);5subplot(221);6stem(n,xn);7title('原信号');8subplot(212);9stem(n,abs(Xk));10title('FFT变换')效果图：分析：由图可见，fft变换的频率中心不在0点，这是fft算法造成的，把fft改为fftshift可以将频率中心移到0点。5Matlab的傅里叶变换实例1.傅里叶变换的定义傅里叶变换从数学上的定义，为F(w)=int(x*exp(-jwt),t=-inf...inf)其中，int表示积分，t是时间，x是时域信号，inf表示无穷，exp表示指数运算。其含义说的是给一个无限长的时域信号和一个频点w，可以唯一确定一个复数F。于是，F和w就有了这种对应关系，考虑到F是个复数。F的绝对值和w关系叫幅频，F的幅角和w关系叫相频。2.matlab的fft命令helpfft可以知道这个和数学上的傅里叶不一样，因为计算机是离散的！因为计算机的时域信号存储量是有限的！比如等时采样得到的信号，高频分量是不可能获得的，对于比较大的w将无法计算。于是，fft这样计算傅里叶变换：把时域信号进行周期延拓，取一组w，就是时域信号的周期及该周期的二分之一，三分之一，直到n分之一，其中n是一个周期内的采样点。这样做的结果，就是对一段有限长的时域信号，将其长度作为基频率，分析基频和高频含量。当然，能分析到的最高频为n次谐波，再高次谐波由于香农定理而无法体现。3.写一个数学定义傅里叶变换的程序将有限长时域信号不延拓，时域信号外的时间内，认为信号为零。于是获得无限长时域信号，取频点若干，分析其傅里叶变换。考虑到matlab对于由描点法定义的函数，数值积分时常用的方法有：矩形法，6梯形法。一下代码采用梯形法，算例如下：clearclc%%输入信号t=0:1e-3:20;%时域信号的时间范围x=sin(t)+sin(1.5*t+1)+5*cos(0.5*t)+2*randn(size(t));%时域信号xw=[0:1e-2:2];%想要观察的频率范围%%预定义y=w;a=w;j=sqrt(-1);%先定义变量维度，提高运算速度%%计算频点fori=1:length(w)f=trapz(t,x.*exp(-j*w(i)*t));y(i)=abs(f);a(i)=angle(f);end%%输出subplot(3,1,1),plot(t,x)subplot(3,1,2),plot(w,y)subplot(3,1,3),plot(w,a)算例中的时域信号，里有三个正余弦分量，一个干扰分量。等时采样，并认为采样频率足够高，即得到的信号是连续信号。7Matlab方波傅立叶变换画出频谱图代码如下：clc;clear;closeall;fs=30;%采样频率T=1/fs;t=0:T:2*pi;A=2;P=4;y=A*square(P*t);subplot(2,1,1),plot(t,y)title('方波信号')Fy=abs(fft(y,512));f2=fs*(0:256)/512;subplot(2,1,2),plot(f2,Fy(1:257))title('频谱图');set(gcf,'unit','normalized','position',[0011]);set(gca,'xtick',0:0.6:8);axis([0,8,0300]);方波大概0.6左右。角频率W=4，所以频率f=4/(2*pi)=2/pi没错。求信号的512点FFT，点数越大(2的整数幂)，频谱越精细。把512改成1024后，频率变到1.2左右了。代码：clc;clear;closeall;fs=30;%采样频率N=input('输入FFT点数：');T=1/fs;t=0:T:2*pi;8A=2;P=4;y=A*square(P*t);subplot(2,1,1),plot(t,y)title('方波信号')Fy=abs(fft(y,N));f2=fs*(0:N/2)/N;subplot(2,1,2),plot(f2,Fy(1:(N/2+1)))title('频谱图');set(gcf,'unit','normalized','position',[0011]);set(gca,'xtick',0:0.6:8);axis([0,8,0300]);输入N分别为128，256，512，1024看看，是不是点越多，曲线越平滑。