MATLAB FFT xaxis限制了混乱和fftshift [英] MATLAB FFT xaxis limits messing up and fftshift
问题描述
这是我第一次使用fft函数，并且试图绘制一个简单余弦函数的频谱:
This is the first time I'm using the fft function and I'm trying to plot the frequency spectrum of a simple cosine function:
f = cos(2 * pi * 300 * t)
f = cos(2*pi*300*t)
采样率为220500.我正在绘制函数f的一秒钟.
The sampling rate is 220500. I'm plotting one second of the function 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)(但是我认为他们只进行了单面频谱分析).
Why do I get odd results when I increase the frequency to 900Hz? These odd results can be fixed by increasing the x-axis limits from, say, 500Hz to 1000Hz. Also, is this the correct approach? I noticed many other people didn't use fftshift(X) (but I think they only did a single sided spectrum analysis).
谢谢.
推荐答案
这是我的承诺.
关于为什么您将频率增加到900 Hz时会得到奇怪的结果"的第一个或您的问题与@Castilho所描述的Matlab的绘图重新缩放功能有关.当您更改x轴的范围时，Matlab将尝试提供帮助并重新缩放y轴.如果峰在指定范围之外，则matlab会放大该过程中产生的微小数字误差.如果困扰您，可以使用"ylim"命令对此进行补救.
The first or your questions related to why you "get odd results when you increase the frequency to 900 Hz" is related to the Matlab's plot rescaling functionality as described by @Castilho. When you change the range of the x-axis, Matlab will try to be helpful and rescale the y-axis. If the peaks lie outside of your specified range, matlab will zoom in on the small numerical errors generated in the process. You can remedy this with the 'ylim' command if it bothers you.
但是，您的第二个更开放的问题是这是正确的方法吗?"需要更深入的讨论.请允许我告诉您如何实现更灵活的解决方案，以实现绘制余弦波的目标.
However, your second, more open question "is this the correct approach?" requires a deeper discussion. Allow me to tell you how I would go about making a more flexible solution to achieve your goal of plotting a cosine wave.
您将从以下内容开始:
time = 1;
freq = 220500;
这立刻引起了我的警觉.查看文章的其余部分，您似乎对低于kHz范围内的频率感兴趣.如果是这种情况，则该采样率过高，因为该速率的奈奎斯特极限(sr/2)高于100 kHz.我猜您是要使用22050 Hz的常见音频采样率(但是我在这里可能是错误的)?
This raises an alarm in my head immediately. Looking at the rest of the post, you appear to be interested in frequencies in the sub-kHz range. If that is the case, then this sampling rate is excessive as the Nyquist limit (sr/2) for this rate is above 100 kHz. I'm guessing you meant to use the common audio sampling rate of 22050 Hz (but I could be wrong here)?
无论哪种方式，您的分析最终都能在数值上得到结果.但是，您并不能帮助您了解如何在实际情况下最有效地使用FFT进行分析.
Either way, your analysis works out numerically OK in the end. However, you are not helping yourself to understand how the FFT can be used most effectively for analysis in real-world situations.
允许我发布如何执行此操作.以下脚本几乎可以完全执行您的脚本，但是打开了我们可以建立基础的潜力. .
Allow me to post how I would do this. The following script does almost exactly what your script does, but opens some potential on which we can build . .
%// 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大小.
You calculate a time axis and calculate your number of FFT points from the length of the time axis. This is very odd. The problem with this approach, is that the frequency resolution of the fft changes as you change the duration of your input signal, because N is dependent on your "time" variable. The matlab fft command will use an FFT size that matches the size of the input signal.
在我的示例中，我直接从NFFT计算频率轴.在上述示例的上下文中，这有点无关紧要，因为我将NFFT设置为等于信号中的采样数.但是，使用这种格式有助于使您的想法神秘化，在我的下一个示例中，它变得非常重要.
In my example, I calculate the frequency axis directly from the NFFT. This is somewhat irrelevant in the context of the above example, as I set the NFFT to equal the number of samples in the signal. However, using this format helps to demystify your thinking and it becomes very important in my next example.
**侧面注意:您在示例中使用real(F).除非有充分的理由仅提取FFT结果的实部，否则使用abs(F)提取FFT的幅度会更加普遍.等同于sqrt(real(F).^ 2 + imag(F).^ 2).**
** SIDE NOTE: You use real(F) in your example. Unless you have a very good reason to only be extracting the real part of the FFT result, then it is much more common to extract the magnitude of the FFT using abs(F). This is the equivalent of sqrt(real(F).^2 + imag(F).^2).**
大多数时候，您将需要使用较短的NFFT.这可能是因为您可能正在实时系统中运行分析，或者是因为您希望对多个FFT的结果求平均，以了解随时间变化的信号的平均频谱，或者是因为您想比较具有不同持续时间而又不浪费信息的信号.仅使用值为NFFT< fft的fft命令.信号中元素的数量将导致根据信号的最后NFFT点计算出fft.这有点浪费.
Most of the time you will want to use a shorter NFFT. This might be because you are perhaps running the analysis in a real time system, or because you want to average the result of many FFTs together to get an idea of the average spectrum for a time varying signal, or because you want to compare spectra of signals that have different duration without wasting information. Just using the fft command with a value of NFFT < the number of elements in your signal will result in an fft calculated from the last NFFT points of the signal. This is a bit wasteful.
以下示例与有用的应用程序更为相关.它显示了如何将信号分成多个块，然后处理每个块并对结果求平均:
The following example is much more relevant to useful application. It shows how you would split a signal into blocks and then process each block and average the result:
%//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
I use a hamming window in the above example. The window that you choose should suit the application http://en.wikipedia.org/wiki/Window_function
您选择的重叠量将在一定程度上取决于您使用的窗口类型.在上面的示例中，汉明窗口将每个缓冲区中的样本权重从零移到远离每个帧中心的位置.为了使用输入信号中的所有信息，重要的是要使用一些重叠.但是，如果仅使用普通的矩形窗口，则由于所有样本的权重均等，因此重叠变得毫无意义.您使用的重叠越多，计算平均光谱所需的处理就越多.
The overlap amount that you choose will depend somewhat on the type of window you use. In the above example, the Hamming window weights the samples in each buffer towards zero away from the centre of each frame. In order to use all of the information in the input signal, it is important to use some overlap. However, if you just use a plain rectangular window, the overlap becomes pointless as all samples are weighted equally. The more overlap you use, the more processing is required to calculate the mean spectrum.
希望这有助于您理解.
这篇关于MATLAB FFT xaxis限制了混乱和fftshift的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
