MATLAB FFT x 轴限制混乱和 fftshift

Question

这是我第一次使用 fft 函数，我试图绘制一个简单的余弦函数的频谱：

f = cos(2*pi*300*t)

采样率为 220500。我正在绘制函数 f 的一秒。

这是我的尝试：

time = 1;
freq = 220500;
t = 0 : 1/freq : 1 - 1/freq;
N = length(t);
df = freq/(N*time);

F = fftshift(fft(cos(2*pi*300*t))/N);
faxis = -N/2 / time : df : (N/2-1) / time;

plot(faxis, real(F));
grid on;
xlim([-500, 500]);

为什么当我将频率增加到 900Hz 时会得到奇怪的结果？这些奇怪的结果可以通过将 x 轴限制从 500Hz 增加到 1000Hz 来解决。另外，这是正确的方法吗？我注意到许多其他人没有使用fftshift(X)（但我认为他们只进行了单面频谱分析）。

谢谢。

Answer 1

这是我承诺的回复。

第一个或你的问题与为什么“当你将频率增加到 900 Hz 时得到奇怪的结果”有关，这与 @Castilho 描述的 Matlab 的绘图重新缩放功能有关。当您更改 x 轴的范围时，Matlab 将尝试提供帮助并重新调整 y 轴。如果峰值超出您指定的范围，matlab 将放大该过程中产生的小数值误差。如果它困扰您，您可以使用 'ylim' 命令解决此问题。

但是，您的第二个更开放的问题“这是正确的方法吗？”需要更深入的讨论。请允许我告诉您，我将如何制定更灵活的解决方案来实现您绘制余弦波的目标。

你从以下开始：

time = 1;
freq = 220500;

这立即在我的脑海中引发了警报。查看帖子的其余部分，您似乎对 sub-kHz 范围内的频率感兴趣。如果是这种情况，则该采样率过高，因为该采样率的奈奎斯特限制 (sr/2) 高于 100 kHz。我猜你的意思是使用 22050 Hz 的常见音频采样率（但我在这里可能错了）？

无论哪种方式，您的分析最终在数值上都可以。但是，您并不能帮助自己了解如何最有效地使用 FFT 进行实际情况的分析。

请允许我发布我将如何做到这一点。以下脚本几乎与您的脚本完全一样，但打开了我们可以构建的一些潜力。 .

%// These are the user parameters
durT = 1;
fs = 22050;
NFFT = durT*fs;
sigFreq = 300;

%//Calculate time axis
dt = 1/fs;
tAxis = 0:dt:(durT-dt);

%//Calculate frequency axis
df = fs/NFFT;
fAxis = 0:df:(fs-df);

%//Calculate time domain signal and convert to frequency domain
x = cos(  2*pi*sigFreq*tAxis  );
F = abs(  fft(x, NFFT)  /  NFFT  );

subplot(2,1,1);
plot(  fAxis, 2*F  )
xlim([0 2*sigFreq])
title('single sided spectrum')

subplot(2,1,2);
plot(  fAxis-fs/2, fftshift(F)  )
xlim([-2*sigFreq 2*sigFreq])
title('whole fft-shifted spectrum')

您计算时间轴并根据时间轴的长度计算您的 FFT 点数。这很奇怪。这种方法的问题在于，fft 的频率分辨率会随着您更改输入信号的持续时间而变化，因为 N 取决于您的“时间”变量。 matlab fft 命令将使用与输入信号大小匹配的 FFT 大小。

在我的示例中，我直接从 NFFT 计算频率轴。这在上述示例的上下文中有些无关紧要，因为我将 NFFT 设置为等于信号中的样本数。但是，使用这种格式有助于揭开你的思维的神秘面纱，这在我的下一个示例中变得非常重要。

** 旁注：您在示例中使用 real(F)。除非您有充分的理由只提取 FFT 结果的实部，否则使用 abs(F) 提取 FFT 的幅度更为常见。这相当于 sqrt(real(F).^2 + imag(F).^2).**

大多数情况下，您会希望使用较短的 NFFT。这可能是因为您可能正在实时系统中运行分析，或者因为您希望将许多 FFT 的结果平均在一起以了解随时间变化的信号的平均频谱，或者因为您希望比较具有不同持续时间的信号，而不会浪费信息。只需使用值为 NFFT

以下示例与有用的应用程序更相关。它展示了如何将信号拆分为块，然后处理每个块并平均结果：

%//These are the user parameters
durT = 1;
fs = 22050;
NFFT = 2048;
sigFreq = 300;

%//Calculate time axis
dt = 1/fs;
tAxis = dt:dt:(durT-dt);

%//Calculate frequency axis
df = fs/NFFT;
fAxis = 0:df:(fs-df);

%//Calculate time domain signal 
x = cos(  2*pi*sigFreq*tAxis  );

%//Buffer it and window
win = hamming(NFFT);%//chose window type based on your application
x = buffer(x, NFFT, NFFT/2); %// 50% overlap between frames in this instance
x = x(:, 2:end-1); %//optional step to remove zero padded frames
x = (  x' * diag(win)  )'; %//efficiently window each frame using matrix algebra

%// Calculate mean FFT
F = abs(  fft(x, NFFT)  /  sum(win)  );
F = mean(F,2);

subplot(2,1,1);
plot(  fAxis, 2*F  )
xlim([0 2*sigFreq])
title('single sided spectrum')

subplot(2,1,2);
plot(  fAxis-fs/2, fftshift(F)  )
xlim([-2*sigFreq 2*sigFreq])
title('whole fft-shifted spectrum')

我在上面的例子中使用了一个汉明窗。您选择的窗口应该适合应用程序http://en.wikipedia.org/wiki/Window_function

您选择的重叠量在一定程度上取决于您使用的窗口类型。在上面的示例中，汉明窗将每个缓冲区中的样本从每帧的中心向零加权。为了使用输入信号中的所有信息，使用一些重叠很重要。但是，如果您只使用一个普通的矩形窗口，则重叠变得毫无意义，因为所有样本的权重均等。您使用的重叠越多，计算平均光谱所需的处理就越多。

希望这有助于您的理解。

Answer 2

你的结果是完全正确的。你的频率轴计算也是正确的。问题出在 y 轴刻度上。当您使用函数 xlims 时，matlab 会自动重新计算 y 比例，以便您可以看到“有意义”的数据。当余弦峰值超出您选择的限制时（当 f>500Hz 时），没有峰值可显示，因此比例是根据一些非常小的噪声计算的（在我的计算机上，使用 matlab 2011a，y 比例为 10 -16)。

改变极限确实是正确的做法，因为如果你不改变它，你就看不到频谱上的峰值。

不过，我注意到了一件事情。你有理由绘制变换的实部吗？通常，绘制的是abs(F)，而不是真实的部分。

编辑：实际上，您的频率轴是正确的，因为在这种情况下，df 是 1。faxis 是正确的，但 df 计算不是。

FFT 计算从 -Fs/2 到 Fs/2 的 N 个点。因此，在 Fs 范围内的 N 个点产生 Fs/N 的 df。因为 N/时间 = Fs => 时间 = N/Fs。将其替换为您使用的 df 表达式：your_df = Fs/N*(N/Fs) = (Fs/N)^2。由于Fs/N = 1，最后的结果是对的：P