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、反变换得到的频谱既有实部也有虚部，我们只取实部。3、 实验过程：具体见如下程序%地球物理信息处理11月15日上机作业clcclfclear all%绘制波形图N=500; %500个采样点ddt=0.002; %采样间隔是0.002dt=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; endendfzhen=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; endendf1zhen=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 f2c=f2(n)*cos(-2*pi/N*(k-1)*(n-1); real_f2(k)=real_f2(k)+f2c; f2s=f2(n)*sin(-2*pi/N*(k-1)*(n-1); imag_f2(k)=imag_f2(k)+f2s; endendf2zhen=sqrt(real_f2.2+imag_f2.2); %f2的振幅f2xiang=atan(real_f2./imag_f2); %f2的相位figure(4) %f2的图像，振幅谱和相位谱subplot(311);plot(t,f2);title(f2的波形图);subplot(312);plot(f2zhen);title(f2的振幅谱);subplot(313);plot(f2xiang);title(f2的相位谱);%对f3变换，并绘制波形图，振幅谱和相位谱，见图五for k=1:N real_f3(k)=0; imag_f3(k)=0; for n=1:N f3c=f3(n)*cos(-2*pi/N*(k-1)*(n-1); real_f3(k)=real_f3(k)+f3c; f3s=f3(n)*sin(-2*pi/N*(k-1)*(n-1); imag_f3(k)=imag_f3(k)+f3s; endendf3zhen=sqrt(real_f3.2+imag_f3.2); %f3的振幅f3xiang=atan(real_f3./imag_f3); %f3的相位figure(5) %f3的图像，振幅谱和相位谱subplot(311);plot(t,f3);title(f3的波形图);subplot(312);plot(f3zhen);title(f3的振幅谱);subplot(313);plot(f3xiang);title(f3的相位谱); %设置滤波器把50Hz的成分滤掉fx=50; %fx表示要滤掉的频率成分df=1/(N*ddt);K=fx/df;lvbo=ones(1,N);lvbo(:,K+1)=0; %把50Hz的K值充零lvbo(:,(N-K+1)=0; %将与50Hz对称的值也充零%设置滤波器：不仅在50Hz的地方充零，在50Hz左右各两个点都充零，对比效果% lvbo(:,K+1-2:K+1+2)=0; %把50Hz的K值左右充零% lvbo(:,(N-K+1)-2:(N-K+1)+2)=0;%将与50Hz对称的值左右也充零fn_lvbohou=fzhen.*lvbo; %滤波fn_lvbohou_z=fn.*lvbo;fn_lvbohou_zhen=sqrt(real(fn_lvbohou_z).2+(imag(fn_lvbohou_z).2);fn_lvbohou_xiang=atan(imag(fn_lvbohou_z)./real(fn_lvbohou_z);%绘制滤波后的振幅谱figure(6)subplot(311);plot(t,fzhen);title(滤波前的振幅谱);subplot(312);plot(t,fn_lvbohou);title(滤波后的振幅谱);subplot(313);plot(t,fn_lvbohou_xiang);title(滤波后的相位谱);axis(0 N -2 2); %将滤波后的频谱反变换回时间域，见图七for n=1:N ifn(n)=0; for k=1:Nifn_lvbohou=fn_lvbohou(k)*exp(1i*2.*pi/N*(k-1)*(n-1); ifn(n)=(ifn(n)+ifn_lvbohou); end ifn(n)=ifn(n)./N;endfigure(7) %对比滤波后反变换的波形图与直接在原始信号中去掉50Hz成分的波形图subplot(311);plot(t,ifn);title(滤波后反变换得到的波形图);subplot(312);plot(t,f1+f3);title(原始信号去掉50Hz的波形图);subplot(313);plot(t,f);title(原始信号的波形图);axis(0 N -3 7); 4、 资料分析：1. 根据奈奎斯特采样定理，采样频率dt1/(2*fmax)，在此我取了ddt=0.002。2. 总信号含有三个部分，分别是直流分量，含有50Hz成分，75Hz成分的余弦信号。做离散时间傅氏变换的时候将实虚部分开来计算，这样可以方便的计算振幅谱和相位谱。在实现dft的过程中，我虽然理解公式的意义和思路，但在matlab编写实现的过程中仍然遇到一些问题，主要就是我所写的程序没有很好的表达出我要表达的意思，可见对matlab的掌握还不够熟练。对这一部分，我思考了许多做法。（1） 、用矩阵实现：Xk=xn*exp(-i*2*pi/N).(n*k); （2） 、用两重for循环计算，做累加。如上述程序所示。（3） 、用了sum函数：for k=1:N ffc(1:N)=cos(-2.*pi./N.*(k-1).*t); %ffc是e的欧拉展开的实部 ffs(1:N)=sin(-2.*pi./N.*(k-1).*t); %ffs是e的欧拉展开的虚部 %总信号的实虚部分别做dft fnc(k)=sum(f.*ffc); %实部做dft fns(k)=sum(f.*ffs); %虚部做dftfzhen=sqrt(fnc.2+fns.2); %总信号的振幅谱fxiang=atan(fns./fnc); fn=fzhen.*exp(i*fxiang);sum函数是我一开始遍的程序，后来才改为比较大众的for循环。实质是一样的，但是实现形式不同。对于实现这些思路，都是我一些想法，可能有不成熟的地方。 3. 坐标横轴是1N，其实代表着0*2pi。其中1N/2代表0pi，而N/2N代表-pi0。 4. 在初始参数的设置中，尽量设采样点数为N，而不要设置时间节点。五、结论：分析如何获得高分辨率谱与减弱截断效应。1. 高分辨率谱：此信号里含有两个频率成分f2=50hz和f3=75hz，ff=|f2-f1|=25hz，频率分辨率的概念是如果频率分辨率大于f