我需要在调用fft或ifft之前调用fftshift吗?

16
在David Voelz所著的《计算傅里叶光学,Matlab教程》中写道,在调用fftifft之前需要调用fftshift命令。但是在MATLAB fftshift文档中只写到这个命令会移动零频率分量到数组中心重新排列fftfft2fftn的输出结果。文档中没有提到必须在调用fft之前调用这个命令,并且我看到一些例子在没有调用fftshift的情况下调用了fft
我的问题是:是否需要在调用fft或者ifft之前调用fftshift
如果调用fft之前不需要调用fftshift,那么我们应该何时(如果需要)使用ifftshift命令来计算数据集的fft

调用 fft 函数无论你是否先调用 fftshift 函数,都会返回相同的结果?这是如何可能的? - user4861528
7个回答

14

为了节省计算时间,Matlab的fft仅计算频谱的一半(正频率和零频率,如果你的样本数是奇数),然后将第二半频谱(即第一半的复共轭)添加到该向量的末尾。

因此,在执行fft之后,您将得到以下向量:

 0 1 2 3 ... Fs/2   -Fs/2 ... -3 -2 -1
<----------------> <------------------>
  positive freq        negative freq

其中Fs是频率采样。

现在,fftshift的作用就是将负频率的数据(频谱的第二部分)移到向量的开头,这样您可以显示一个从-Fs/2开始,到+Fs/2结束的漂亮频谱。移位后的向量变为:

 -Fs/2 ... -3 -2 -1   0 1 2 3 ... Fs/2
<------------------> <---------------->
    negative freq      positive freq

所以,回答你的问题,不需要在调用fftifft前后使用fftshift。但如果你对向量应用了fftshift,则应该通过应用ifftshiftfftshift来撤消它。(我认为这两种调用是等效的。)


2
MATLAB产生FFT输出所进行的计算是无关紧要的。FFT的输出由DFT的定义给出,其中频率k=0..N-1。在此输出中没有“负频率”。DFT是周期性的,这意味着k=0处的值与k=N和k=-N+1处的值相同。有些人更喜欢查看k=-N/2..N/2的输出,fftshift函数提供了这种能力。ifftshiftfftshift不同,它们的功能对于奇数大小的输入不同。 - Cris Luengo
此页面的附录和评论中有一个关于fftshift和ifftshift的不错讨论。尽管这是针对numpy的,但我认为它也适用于Matlab。 - ahmadh

5

如果你在fftshift的文档中往后阅读:

"它对于将零频率分量放在频谱中央进行傅里叶变换的可视化非常有用。"

移位对于执行fft并不是必要的,但它对于可视化傅里叶变换非常方便。因此,是否使用fftshift取决于您是否想要可视化变换。


我猜当你将数组转换为绘图并将其分配到原始数组而不是单独的变量中时。 - Adriaan
调用 fft 函数后,无论之前是否调用了 fftshift 函数,返回的结果都是相同的。这是如何可能的? - user4861528
我刚刚告诉你,你为了绘图而移动了变换,现在希望将数组移回以进行进一步操作。但是,更好的做法是将变换分配给一个单独的数组,然后不使用反向变换,而是使用原始数组。 - Adriaan

3
请注意,ifftshiftfftshift不同,因为它将负数移回到正数。在进行fftshift之前,假设在频域中存在一个简单的3个bin单位脉冲。
[0, exp(-jw1), exp(-(jw1-pi)=exp(-jw1+pi)];

fftshift函数的作用

[exp(-jw1+pi)], 0, exp(-jw1)];

您可以看到相位上的不连续性。如果执行ifftshift,则负频率将被移回正频率:

[0, exp(-jw1), exp(-jw1+pi)];

然而,fftshift 再次给出:

[exp(-jw1), exp(-jw1+pi), 0];

我们可以看到相位单调性不同,即fftshift-ifftshift情况下相位呈现[减少,增加],而fftshift-fftshift情况下相位呈现[增加,减少]。


2
考虑到您正在学习的书名,我猜测您正在处理图像。那么正确的方法是在调用[i]fft之前和之后都调用fftshift和ifftshift。
您的代码应该如下所示:
spectrum = fftshift(fft2(ifftshift(myimage))

当应用傅里叶逆变换时,情况基本相同。

myimage = fftshift(ifft2(ifftshift(spectrum))

这是一个很好的解释:https://groups.google.com/g/comp.soft-sys.matlab/c/rUcc0bRRZf4 - Amel Alhassan

1
简单的答案是在调用fft之前不需要调用fftshift。fftshift不会影响快速傅里叶变换的计算。fftshift重新排列矩阵内部的值。例如:
cameraman = imread('cameraman.tif');
fftshifted_cameraman = fftshift(cameraman);
subplot(1,2,1); imshow(cameraman); title('Cameraman');
subplot(1,2,2); imshow(fftshifted_cameraman); title('FFTShifted Cameraman');

4
但通过重新排列数值,它确实会影响结果。 - Cris Luengo

0
这是我用来测试MATLAB基本声音分析函数的脚本,包括fft()的输出使用fftshift()进行显示。
if ~exist('inputFile', 'var')
    inputFile = 'vibe.wav';
end

[inputBuffer, Fs] = audioread(inputFile);

fileSize = length(inputBuffer);

numSamples = 2.^(ceil(log2(fileSize))); % round up to nearest power of 2

x = zeros(numSamples, 1);                   % zero pad if necessary

x(1:fileSize) = inputBuffer(:,1);           % if multi-channel, use left channel only

clear inputBuffer;                          % free this memory
clear fileSize;

t = linspace(0, (numSamples-1)/Fs, numSamples)';
f = linspace(-Fs/2, Fs/2 - Fs/numSamples, numSamples)';

X = fft(x);

plot(t, x);
xlabel('time (seconds)');
ylabel('amplitude');
title(['time-domain plot of ' inputFile]);
sound(x, Fs);                                           % play the sound
pause;




% display both positive and negative frequency spectrum

plot(f, real(fftshift(X)));
xlabel('frequency (Hz)');
ylabel('real part');
title(['real part frequency-domain plot of ' inputFile]);
pause;

plot(f, imag(fftshift(X)));
xlabel('frequency (Hz)');
ylabel('imag part');
title(['imag part frequency-domain plot of ' inputFile]);
pause;

plot(f, abs(fftshift(X)));                              % linear amplitude by linear freq plot
xlabel('frequency (Hz)');
ylabel('amplitude');
title(['abs frequency-domain plot of ' inputFile]);
pause;

plot(f, 20*log10(abs(fftshift(X))+1.0e-10));            % dB by linear freq plot
xlabel('frequency (Hz)');
ylabel('amplitude (dB)');
title(['dB frequency-domain plot of ' inputFile]);
pause;





% display only positive frequency spectrum for log frequency scale

semilogx(f(numSamples/2+2:numSamples), 20*log10(abs(X(2:numSamples/2))));       % dB by log freq plot
xlabel('frequency (Hz), log scale');
ylabel('amplitude (dB)');
title(['dB vs. log freq, frequency-domain plot of ' inputFile]);
pause;

semilogx(f(numSamples/2+2:numSamples), (180/pi)*angle(X(2:numSamples/2)));      % phase by log freq plot
xlabel('frequency (Hz), log scale');
ylabel('phase (degrees)');
title(['phase vs. log freq, frequency-domain plot of ' inputFile]);
pause;

%
%   this is an alternate method of unwrapping phase
%
%   phase = cumsum([angle(X(1)); angle( X(2:numSamples/2) ./ X(1:numSamples/2-1) ) ]);  
%   semilogx(f(numSamples/2+2:numSamples), phase(2:numSamples/2));                  % unwrapped phase by log freq plot
%

semilogx(f(numSamples/2+2:numSamples), unwrap(angle(X(2:numSamples/2))));       % unwrapped phase by log freq plot
xlabel('frequency (Hz), log scale');
ylabel('unwrapped phase (radians)');
title(['unwrapped phase vs. log freq, frequency-domain plot of ' inputFile]);

如果您正在将音频分段窗口处理并将其传递给FFT,则应在FFT的输入上使用fftshift(),以将窗口化段的中心定义为t = 0点。

[x_input, Fs] = audioread('vibe.wav');     % load time-domain input
N = 2*floor(length(x_input)/2);            % make sure N is even
x = x_input(1:N);

t = linspace(-N/2, N/2-1, N);              % values of time in units of samples
omega = linspace(-pi, pi*(1-2/N), N);      % values of (normalized) angular frequency

X = fftshift( fft( fftshift( x.*hamming(length(x)) ) ) );

[X_max k_max] = max( abs(X) );

figure(1);
plot(t, x, 'g');

figure(2);
plot(omega, abs(X), 'b');
hold on;
plot(omega(k_max), X_max, 'or');
hold off;

0

这取决于您将如何处理转换后的数据。如果在变换之前不执行fftshift,则fft结果将使每个其他值乘以-1。如果您计划查看结果的幅度或幅度平方,则这并不重要。但是,如果您计划比较相邻的频谱值并且相位很重要,则需要在变换之前应用fftshift以避免交替符号。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接