在倍频程上绘制 FFT

Question

我知道FFT将时域中的函数更改为频域中显示的函数。

但是，当我尝试在频域中绘制所述图时，我只能通过使用时间作为 X 轴来使其正常工作，而显然它应该不是那个，而是频率。

此外，我只能通过将 y 轴除以某个整数来获得与原始信号中的幅度相匹配的幅度。这是为什么？

这是我的代码

t=0:0.001:2

x=2*sin(20*pi*t) + sin(100*pi*t)
subplot(2,1,1)
plot(1000*t,x)
grid
xlabel("Time in milliseconds")
ylabel("Signal amplitude")

subplot(2,1,2)
y=fft(x)
plot(1000*t,abs(y))
xlabel("Frequency")
ylabel("Signal amplitude")

和图表。

请帮忙=(

Answer 1

频率关系（x轴缩放）

直到奈奎斯特频率（采样率的一半），FFT 产生的每个值的频率通过以下方式与输出值的索引线性相关：

f(i) = (i-1)*sampling_frequency/N

其中 N 是 FFT 点数（即

N=length(y)

）。就你而言，

N=2001

。在奈奎斯特频率之上，频谱显示出环绕负频率分量（来自频谱的周期性扩展）。

可以从

t

的定义中扣除采样频率 1/T，其中 T 是采样时间间隔（在您的情况下 T=0.001）。 所以采样频率为1000Hz。

注意，由于

t(i)

的值也与索引

i

线性相关，通过

t(i) = (i-1)*0.001

可以定义

f = 1000*t*sampling_frequency/N

（尽管不一定建议，因为这只会使你的代码变得模糊）。 请注意，您错过了

sampling_frequency/N

项，这相应地导致音调以错误的频率显示 （根据

x

的定义，在 10Hz 和 50Hz 处应该有峰值，在 -10Hz 和 -50Hz 处应该有相应的别名，环绕后会出现在 990Hz 和 950Hz 处）。

幅度关系（y轴缩放）

请注意，观察到的关系只是近似的，因此以下不是数学证明，而只是直观地可视化时域音调幅度和频域峰值之间的关系。

将问题简化为单一音调：

x = A*sin(2*pi*f*t)

可以使用 Parseval 定理推导出相应峰值的近似幅度：

在时域（等式左边），表达式近似等于

0.5*N*(A^2)

。

在频域（等式右侧），做出以下假设：

• 光谱泄漏影响可以忽略不计
• 音调的频谱内容仅包含在 2 个箱中（频率

  f

  和相应的混叠频率

  -f

  ）占总和（所有其他箱约为 0）。请注意，这通常仅在音调频率是

  sampling_frequency/N

  的精确（或接近精确）倍数时成立。

对于对应于频率

2*(1/N)*abs(X(k))^2

处的峰值的某个

k

值，右侧的表达式近似等于

f

。

将两者放在一起产生

abs(X(k)) ~ 0.5*A*N

。换句话说，正如您所观察到的，输出幅度显示相对于时域幅度的缩放因子

0.5*N

（在您的情况下约为 1000）。

这个想法仍然适用于多个音调（尽管可忽略的频谱泄漏假设最终不成立）。

Answer 2

其他答案表明，本例中的频率响应为 950Hz 和 990Hz。这是对 FFT 代码如何使用索引的误解。这些“高频”尖峰实际上是 -50Hz 和 -10Hz。

频域从-N/2*sampling_Frequency/N延伸到+N/2*Sampling_Frequency/N。但由于历史原因，惯例是前N/2条信息是正频率，中点是零频率，最后N/2条信息倒序是负频率。对于功率谱，不需要显示超过前 1+N/2 条信息。

这个惯例非常令人困惑，因为我不得不从 Press 等人那里弄清楚它。数值配方和手动编码快速哈特利变换，很多年前，当我第一次使用 FFT 时，早于 Cleve Moler 向一些幸运博士生分发的 Matlab 1.0 Beta 测试版本:-)

Answer 3

为了完成并说明其他答案，此代码受 Matlab 文档启发，适用于 Octave：

通过定义@SleuthEye精确的频率轴，我们得到了预期的 50Hz 和 120Hz 音调（及其负别名）：

Fs = 1000;
Ts = 1/Fs;
L = 1500;
t = (0:L-1)*Ts;

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));

figure(1)
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')

Y = fft(X);

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

figure(2)
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

f2 = Fs*(0:(L-1))/L;

figure(3)
plot(f2,P2) 
title('Two-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P2(f)|')

Answer 4

谢谢！对于这个非常全面的脚本。 kr，sepp2gl