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;