Matlab中用fft作频谱后为什么要用fftshift

fft是一维傅里叶变换，即将时域信号转换为频域信号
fftshift
是针对频域的，将FFT的DC分量移到频谱中心
即对频域的图像，（假设用一条水平线和一条垂直线将频谱图分成四块）对这四块进行对角线的交换与反对角线的交换

FFTSHIFT Shift zero-frequency component to center of spectrum.
    For vectors, FFTSHIFT(X) swaps（交换） the left and right halves of
    X. For matrices, FFTSHIFT(X) swaps the first and third
    quadrants and the second and fourth quadrants. For N-D
    arrays, FFTSHIFT(X) swaps "half-spaces" of X along each
    dimension.

    FFTSHIFT(X,DIM) applies the FFTSHIFT operation along the 
    dimension DIM.

    FFTSHIFT is useful for visualizing the Fourier transform with
    the zero-frequency component in the middle of the spectrum.

fftshift就是对换数据的左右两边比如
x=[1 2 3 4]
fftshift(x) ->[3 4 1 2]

IFFTSHIFT Inverse FFT shift.（就是fftshift的逆）

x=[1     2     3     4     5];

y=fftshift(x)

y =

     4     5     1     2     3

ifftshift(y)

ans =

     1     2     3     4     5

   
    IFFTSHIFT undoes the effects of FFTSHIFT.

注意：在使用matlab的fft及fftshift时，应注意。

假定采样频率fs，采样间隔dt，采样点数N。

fft后，频率为（0：N-1）/N/dt

进行fftshift后，频率为

if mod(N,2)==0

n1=(0:N-1)-N/2;

else

n1=(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*pi

else

    kx=(0:tempM-1)/tempM/dx-(tempM-1)/2/tempM/dx;% kx=kx*2*pi

end

kx=kx*2*pi;

if mod(tempM,2)==0

    ky=(0:tempN-1)/tempN/dy-tempN/2/tempN/dy;% kx=kx*2*pi

else

    ky=(0:tempN-1)/tempN/dy-(tempN-1)/2/tempN/dy;% kx=kx*2*pi

end

ky=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 ky

if mod(tempM,2)==0

    kx=(0:tempM-1)-tempM/2;% kx=kx*2*pi

else

    kx=(0:tempM-1)-(tempM-1)/2;% kx=kx*2*pi

end

kx=kx*2*pi/tempM/dx;

if mod(tempN,2)==0

    ky=(0:tempN-1)-tempN/2;% kx=kx*2*pi

else

    ky=(0:tempN-1)-(tempN-1)/2;% kx=kx*2*pi

end

ky=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*pi

else

    temp1=(tempM-1)/2;

    temp2=(temp1+1):(tempM-1);

    temp2=temp2-tempM;

    temp3=[0:temp1,temp2];

    kx=temp3/tempM/dx;% kx=kx*2*pi

end

kx=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*pi

else

    temp1=(tempN-1)/2;

    temp2=(temp1+1):(tempN-1);

    temp2=temp2-tempN;

    temp3=[0:temp1,temp2];

    ky=temp3/tempN/dy;% kx=kx*2*pi

end

ky=ky*2*pi;

[kx,ky]=meshgrid(kx,ky);

kx=fftshift(kx);

ky=fftshift(ky);

k=sqrt(kx.^2+ky.^2);

figure

subplot(3,1,1),contourf(kxx-kx)

subplot(3,1,2),contourf(kyy-ky)

subplot(3,1,3),contourf(k00-k)

%%%%%%%%%%%

fft及fftshift示例：

clf;

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;