MATLAB FFT xaxis限制了困惑和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

这是我的 promise 。

关于为什么为什么“将频率增加到900 Hz时会得到奇怪的结果”的第一个或您的问题与@Castilho所描述的Matlab的绘图缩放比例功能有关。当您更改x轴的范围时，Matlab将尝试提供帮助并重新缩放y轴。如果峰在指定范围之外，则matlab会放大该过程中产生的微小数字误差。如果困扰您，可以使用“ylim”命令对此进行补救。

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

您将从以下内容开始:

time = 1;
freq = 220500;

这立刻在我的脑海中发出了警报。查看文章的其余部分，您似乎对低于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点的数量。这很奇怪。这种方法的问题在于，当您更改输入信号的持续时间时，ftf的频率分辨率会发生变化，因为N取决于您的“时间”变量。 matlab fft命令将使用与输入信号大小匹配的FFT大小。

在我的示例中，我直接从NFFT计算频率轴。在上述示例的上下文中，这有点无关紧要，因为我将NFFT设置为等于信号中的采样数。但是，使用这种格式有助于使您的想法神秘化，在我的下一个示例中，它变得非常重要。

**注意:您在示例中使用real(F)。除非有充分的理由仅提取FFT结果的实部，否则使用abs(F)提取FFT的幅度会更加普遍。这等效于sqrt(real(F)。^ 2 + imag(F)。^ 2)。**

大多数时候，您会希望使用较短的NFFT。这可能是因为您可能正在实时系统中运行分析，或者是因为您希望对多个FFT的结果求平均，以了解随时间变化的信号的平均频谱，或者是因为您想比较持续时间不同而又不浪费信息的信号。仅使用值为NFFT <信号中的元素数的fft命令，将产生根据信号的最后NFFT点计算出的fft。这有点浪费。

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

%//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

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

希望这有助于您的理解。