我该如何在Python中生成一个一维信号的频谱图？

Question

我不确定如何做这件事，但是我收到了一个例子，例如频谱图，但是这是一个二维的图像。

我有一段代码生成了混合频率，并且在fft中可以选择这些频率，但是我该如何在频谱图中看到这些频率？我知道我的例子中频率并没有随时间变化，那么这是否意味着我会在频谱图上看到一条直线？

以下是我的代码和输出图片：

# create a wave with 1Mhz and 0.5Mhz frequencies
dt = 2e-9
t = np.arange(0, 10e-6, dt)
y = np.cos(2 * pi * 1e6 * t) + (np.cos(2 * pi * 2e6 *t) * np.cos(2 * pi * 2e6 * t))
y *= np.hanning(len(y))
yy = np.concatenate((y, ([0] * 10 * len(y))))

# FFT of this
Fs = 1 / dt  # sampling rate, Fs = 500MHz = 1/2ns
n = len(yy)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n / 2)]  # one side frequency range
Y = fft(yy) / n  # fft computing and normalization
Y = Y[range(n / 2)] / max(Y[range(n / 2)])

# plotting the data
subplot(3, 1, 1)
plot(t * 1e3, y, 'r')
xlabel('Time (micro seconds)')
ylabel('Amplitude')
grid()

# plotting the spectrum
subplot(3, 1, 2)
plot(frq[0:600], abs(Y[0:600]), 'k')
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
grid()

# plotting the specgram
subplot(3, 1, 3)
Pxx, freqs, bins, im = specgram(y, NFFT=512, Fs=Fs, noverlap=10)
show()

- Harry Lime

---

我期望在频谱图中看到条纹。然而，由于第一个图显示振幅随时间变化，这些条纹应该在中间占主导地位，并在两端逐渐消失。从代码中并不清楚为什么会出现明显的振幅调制。我怀疑波浪乘法中有错别字。 - phkahler

当您发布内容时，请注明需要导入的库以使代码正常运行。谢谢！ - Ravaal

一个 [mcve] 需要在顶部包含 import 语句。 - wizzwizz4

Answer 1

你的代码在技术上是正确的，但你需要查看具有有趣频谱图的信号。为此，你需要频率随时间变化。（为了发生这种情况，你需要许多振荡，因为需要几个振荡来确定一种频率，然后需要许多这样的振荡来使频率以有趣的方式随时间改变。）

下面我尽可能少地修改了你的代码，以获得一个产生有趣效果的频率（fscale只是随着时间逐渐增加频率）。我会发布所有代码以使其正常运行，但我只更改了前四行中的三行。

# create a wave with 1Mhz and 0.5Mhz frequencies
dt = 40e-9
t = np.arange(0, 1000e-6, dt)
fscale = t/max(t)
y = np.cos(2 * pi * 1e6 * t*fscale) + (np.cos(2 * pi * 2e6 *t*fscale) * np.cos(2 * pi * 2e6 * t*fscale))
y *= np.hanning(len(y))
yy = np.concatenate((y, ([0] * 10 * len(y))))

# FFT of this
Fs = 1 / dt  # sampling rate, Fs = 500MHz = 1/2ns
n = len(yy)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n / 2)]  # one side frequency range
Y = fft(yy) / n  # fft computing and normalization
Y = Y[range(n / 2)] / max(Y[range(n / 2)])

# plotting the data
subplot(3, 1, 1)
plot(t * 1e3, y, 'r')
xlabel('Time (micro seconds)')
ylabel('Amplitude')
grid()

# plotting the spectrum
subplot(3, 1, 2)
plot(frq[0:600], abs(Y[0:600]), 'k')
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
grid()

# plotting the specgram
subplot(3, 1, 3)
Pxx, freqs, bins, im = specgram(y, NFFT=512, Fs=Fs, noverlap=10)
show()

此外，需要注意的是只有频谱图是有用的。如果您可以看到一个好的波形或光谱图，那么频谱图可能不会有趣：1）如果波形清晰，您可能没有足够的数据和时间来定义频率并且变化足够大以使其有趣；2）如果整个光谱清晰，您可能没有足够的频率变化使频谱图有趣，因为光谱基本上只是在频谱图中随时间变化所见的东西的平均值。

如果您真的想要看到原始信号的频谱图，您只需要在y轴上缩放以查看您期望的峰值（请注意，频谱图y轴为2.5e8，比您的光谱大得多）。

- tom10

Answer 2

为了实现你的目标：

1）以高频率（至少是其最高频率成分的5倍）对1d波形进行采样。

2）使用样本块（如1024、16384等2的幂次方）来计算FFT。

3）对于每个频谱，绘制一个垂直像素线，其颜色表示每个频率的振幅。

4）对每个样本块重复步骤2和3。

在你的情况下，绘图具有整个彩虹色，这不应该只有几个非常明显的频率。你的频谱图在峰值周围有相当宽的带，但这可能是由于低采样率和平滑绘图所致。

- phkahler

---

谢谢phkahler，不知道你能否提供一个这个的代码框架？最好的，Jake - Harry Lime

@phkahler：我给了你一个+1，你说的如何制作频谱图是正确的，但所有这些都将由matplotlib命令specgram完成，而OP已经在使用它。 - tom10

Answer 3

我刚开始学习Python 3.6

感谢您提供的频谱图示例代码！

然而，使用Python 3.6时，我在使这个示例频谱图代码工作方面遇到了一些困难（函数调用和浮点除法）。 我已经编辑了代码，现在它可以在Python 3.6上运行，适用于我的Python新手朋友们。

享受吧！

'''
Original Script for Python 2.7
https://dev59.com/RXfZa4cB1Zd3GeqPQViV
Modified in August 2017 for Python 3.6
Python 2.7 two integers / Division generate Integer
Python 3.6 two integers / Division generate Float
Python 3.6 two integers // Division generate integer
'''


import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt


dt = 40e-9
t = np.arange(0, 1000e-6, dt)
fscale = t/max(t)
y = np.cos(2 * np.pi * 1e6 * t*fscale) + (np.cos(2 * np.pi * 2e6 *t*fscale) * np.cos(2 * np.pi * 2e6 * t*fscale))
y *= np.hanning(len(y))
yy = np.concatenate((y, ([0] * 10 * len(y))))

# FFT of this
Fs = 1 / dt  # sampling rate, Fs = 500MHz = 1/2ns
n = len(yy)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n // 2)]  # one side frequency range
Y = fftpack.fft(yy) / n  # fft computing and normalization
Y = Y[range(n // 2)] / max(Y[range(n // 2)])

# plotting the data
plt.subplot(3, 1, 1)
plt.plot(t * 1e3, y, 'r')
plt.xlabel('Time (micro seconds)')
plt.ylabel('Amplitude')
plt.grid()

# plotting the spectrum
plt.subplot(3, 1, 2)
plt.plot(frq[0:600], abs(Y[0:600]), 'k')
plt.xlabel('Freq (Hz)')
plt.ylabel('|Y(freq)|')
plt.grid()

# plotting the specgram
plt.subplot(3, 1, 3)
Pxx, freqs, bins, im = plt.specgram(y, NFFT=512, Fs=Fs, noverlap=10)
plt.show()

- Pascal G