MATLAB中的fft后为何要用fftshift?
来源:学生作业帮 编辑:神马作文网作业帮 分类:综合作业 时间:2024/11/20 15:46:36
MATLAB中的fft后为何要用fftshift?
fft是一维傅里叶变换,即将时域信号转换为频域信号fftshift是针对频域的,将FFT的DC分量移到频谱中心\x0d即对频域的图像,(假设用一条水平线和一条垂直线将频谱图分成四块)对这四块进行对角线的交换与反对角线的交换FFTSHIFT Shift zero-frequency component to center of spectrum.\x0dFor vectors,FFTSHIFT(X) swaps(交换) the left and right halves of\x0dX.For matrices,FFTSHIFT(X) swaps the first and third\x0dquadrants and the second and fourth quadrants.For N-D\x0darrays,FFTSHIFT(X) swaps "half-spaces" of X along eachdimension.\x0dFFTSHIFT(X,DIM) applies the FFTSHIFT operation along the\x0ddimension DIM.\x0dFFTSHIFT is useful for visualizing the Fourier transform with\x0dthe zero-frequency component in the middle of the spectrum.fftshift就是对换数据的左右两边比如\x0dx=[1 2 3 4]\x0dfftshift(x) ->[3 4 1 2]\x0dIFFTSHIFT Inverse FFT shift.(就是fftshift的逆)x=[1 2 3 4 5];y=fftshift(x)y = 4 5 1 2 3ifftshift(y)ans = 1 2 3 4 5\x0dIFFTSHIFT undoes the effects of FFTSHIFT.注意:在使用matlab的fft及fftshift时,应注意.假定采样频率fs,采样间隔dt,采样点数N.fft后,频率为(0:N-1)/N/dt进行fftshift后,频率为if mod(N,2)==0n1=(0:N-1)-N/2;elsen1=(0:N-1)-(N-1)/2;end实际上,频率为N点为周期的,所以(0:N-1)所以,对于频率0,1,2,3,4,实际上为0,1,2,-2(3-5),-1(4-5).fftshift后的频率为-2,-1,0,1,2对于二维fftshift,其与直接用下面的结果一样if mod(tempN,2)==0 kx=(0:tempM-1)/tempM/dx-tempM/2/tempM/dx;% kx=kx*2*pielse kx=(0:tempM-1)/tempM/dx-(tempM-1)/2/tempM/dx;% kx=kx*2*piendkx=kx*2*pi;if mod(tempM,2)==0 ky=(0:tempN-1)/tempN/dy-tempN/2/tempN/dy;% kx=kx*2*pielse ky=(0:tempN-1)/tempN/dy-(tempN-1)/2/tempN/dy;% kx=kx*2*piendky=ky*2*pi;temp1=sqrt(kx.^2+ky.^2);k1=temp1;[kx,ky]=meshgrid(kx,ky);如下面程序表明上面两个相同:dx=50e3; dy=50e3;% % % % % % % % % % % tempN=41;tempM=41;% % % % % % % % % % % % % % % %determining the wavenumber kx and kyif mod(tempM,2)==0 kx=(0:tempM-1)-tempM/2;% kx=kx*2*pielse kx=(0:tempM-1)-(tempM-1)/2;% kx=kx*2*piendkx=kx*2*pi/tempM/dx;if mod(tempN,2)==0 ky=(0:tempN-1)-tempN/2;% kx=kx*2*pielse ky=(0:tempN-1)-(tempN-1)/2;% kx=kx*2*piendky=ky*2*pi/tempN/dy;[kxx,kyy]=meshgrid(kx,ky);k00=sqrt(kx.^2+ky.^2);% % % % % % % % % % % % % % % % if mod(tempM,2)==0 temp1=tempM/2-1; temp2=(temp1+1):(tempM-1); temp2=temp2-tempM; temp3=[0:temp1,temp2]; kx=temp3/tempM/dx;% kx=kx*2*pielse temp1=(tempM-1)/2; temp2=(temp1+1):(tempM-1); temp2=temp2-tempM; temp3=[0:temp1,temp2]; kx=temp3/tempM/dx;% kx=kx*2*piendkx=kx*2*pi;if mod(tempN,2)==0 temp1=tempN/2-1; temp2=(temp1+1):(tempN-1); temp2=temp2-tempN; temp3=[0:temp1,temp2]; ky=temp3/tempN/dy;% kx=kx*2*pielse temp1=(tempN-1)/2; temp2=(temp1+1):(tempN-1); temp2=temp2-tempN; temp3=[0:temp1,temp2]; ky=temp3/tempN/dy;% kx=kx*2*piendky=ky*2*pi;[kx,ky]=meshgrid(kx,ky);kx=fftshift(kx);ky=fftshift(ky);k=sqrt(kx.^2+ky.^2);figuresubplot(3,1,1),contourf(kxx-kx)subplot(3,1,2),contourf(kyy-ky)subplot(3,1,3),contourf(k00-k)%%%%%%%%%%%fft及fftshift示例:\x0dclf;fs=100;N=256; %采样频率和数据点数n=0:N-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y1=fft(x,N); %对信号进行快速Fourier变换y2=fftshift(y1);mag1=abs(y1); %求得Fourier变换后的振幅mag2=abs(y2); f1=n*fs/N; %频率序列f2=n*fs/N-fs/2;%这个未必正确subplot(3,1,1),plot(f1,mag1,'r'); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('图1:usual FFT','color','r');grid on;subplot(3,1,2),plot(f2,mag1,'b'); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('图2:FFT without fftshift','color','b');grid on;subplot(3,1,3),plot(f2,mag2,'c'); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('图3:FFT after fftshift','color','c');grid on;