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=[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 3ifftshift(y)
ans =
1 2 3 4 5IFFTSHIFT 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*pielse
kx=(0:tempM-1)/tempM/dx-(tempM-1)/2/tempM/dx;% kx=kx*2*piend
kx=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*piend
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*pielse
kx=(0:tempM-1)-(tempM-1)/2;% kx=kx*2*piend
kx=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*piend
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*pielse
temp1=(tempM-1)/2;temp2=(temp1+1):(tempM-1);temp2=temp2-tempM;temp3=[0:temp1,temp2];kx=temp3/tempM/dx;% kx=kx*2*piend
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*pielse
temp1=(tempN-1)/2;temp2=(temp1+1):(tempN-1);temp2=temp2-tempN;temp3=[0:temp1,temp2];ky=temp3/tempN/dy;% kx=kx*2*piend
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;
總結
以上是生活随笔為你收集整理的MATLAB中的fft后为何要用fftshift?的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: labview 软件编程规范
- 下一篇: 信号互相关及其应用