DFT变换.doc

地球物理信息处理基础 DFT--实验作业报告 一、 实验目的： 1、 信号：f(t)=2+3cos(2∏*50t-30)+1.5cos(2∏*75t+90)。对该信号波形有直观的认识。 2、 利用DFT将时间信号变换至频率域， 并绘制振幅谱、相位谱。 3、 将50Hz频率成分滤去，绘制滤波后的振幅谱、相位谱。 4、 将频谱反变换回时间域，对比前后波形。 5、 分析如何获得高分辨率谱与减弱截断效应。 二、实验方法： 1、利用matlab编写程序。 2、根据原理，不直接使用函数，实现dft变换，dft反变换以及滤波过程。 3、将这部分利用欧拉公式分成实虚部，然后分别做dft。那么变换到频率域的振幅谱就可以直接利用实部的平方加虚部的平方再开根号来求，相位谱也可以直接利用虚部除以实部再求反正切得到。这里需要注意时间域的离散造成频率域的周期化，出现对称现象。 4、滤波过程的关键在于求取50Hz对应的K值，以及对称的K1值，将滤波器充零。 5、反变换得到的频谱既有实部也有虚部，我们只取实部。 实验过程： 具体见如下程序 %地球物理信息处理11月15日上机作业 clc clf clear all %绘制波形图 N=500; %500个采样点 ddt=0.002; %采样间隔是0.002 dt=0:ddt:N*ddt-ddt; df=1/(N*ddt); t=N.*dt; f1=2*ones(1,N); f2=3*cos(2*pi*50*dt-pi/6); f3=1.5*cos(2*pi*75*dt+pi/2); f=f1+f2+f3; %叠加总信号 figure(1) %见图一 subplot(411);plot(t,f1);title(信号f1=2的波形图); subplot(412);plot(t,f2);title(信号f2的波形图); subplot(413);plot(t,f3);title(信号f3的波形图); subplot(414);plot(t,f);title(总信号f的波形图); %对原始信号做dft变换，并绘制波形图，振幅谱和相位谱 for k=1:N real_f(k)=0; imag_f(k)=0; fn(k)=0; for n=1:N %分成实虚部来计算 fc=f(n)*cos(-2*pi/N*(k-1)*(n-1)); real_f(k)=real_f(k)+fc; fs=f(n)*sin(-2*pi/N*(k-1)*(n-1)); imag_f(k)=imag_f(k)+fs; end end fzhen=sqrt(real_f.^2+imag_f.^2); %总信号的振幅谱 fxiang=atan(real_f./imag_f); %总信号的相位谱 fn=fzhen.*exp(1i.*fxiang); figure(2) %总信号的图像，振幅谱和相位谱，见图二 subplot(311);plot(t,f);title(总信号f的波形图); subplot(312);plot(fzhen);title(总信号的振幅谱); subplot(313);plot(fxiang);title(总信号的相位谱); %对f1变换，并绘制波形图，振幅谱和相位谱 for k=1:N real_f1(k)=0; imag_f1(k)=0; for n=1:N f1c=f1(n)*cos(-2*pi/N*(k-1)*(n-1)); real_f1(k)=real_f1(k)+f1c; f1s=f1(n)*sin(-2*pi/N*(k-1)*(n-1)); imag_f1(k)=imag_f1(k)+f1s; end end f1zhen=sqrt(real_f1.^2+imag_f1.^2); %f1的振幅 f1xiang=atan(real_f1./imag_f1); %f1的相位 figure(3) %f1的图像，振幅谱和相位谱，见图三 subplot(311);plot(t,f1);title(f1的波形图); subplot(312);plot(f1zhen);title(f1的振幅谱); subplot(313);plot(f1xiang);title(f1的相位谱); %对f2变换，并绘制波形图，振幅谱和相位谱，见图四 for k=1:N real_f2(k)=0; imag_f2(k)=0; for n=1:N