MATLAB FFT xaxis限制了混乱和移位

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;

这立刻引起了我的警觉.查看帖子的其余部分,您似乎对亚kHz范围内的频率感兴趣.如果是这种情况,则该采样率过高,因为该速率的奈奎斯特极限(sr/2)高于100 kHz.我猜你的意思是使用22050赫兹的普通音频采样率(但我可能在这里错了)？

无论哪种方式,您的分析最终都会在数值上正常运行.但是,您并没有帮助自己了解如何在实际情况下最有效地使用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设置为等于信号中的样本数.但是,使用这种格式有助于揭开思维的神秘面纱,在下一个例子中它变得非常重要.

**侧面注意:您在示例中使用了真实(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

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

希望这有助于您的理解.