Python 2.7 中的 FFT 合并两个代码

Question

我有两个用于 FFT 的代码，但我需要将它们合并，因为它们似乎在某些部分运行良好。让我解释。 第一个代码：

fft1 = (Bx[51:-14])
fft2 = (By[1:-14])

# NS antena FFT - red
FFTdata = np.sqrt(fft1*fft1)**1
samples = FFTdata.size

# WE antena FFT - blue
FFTdata2 = np.sqrt(fft2*fft2)**1
samples2 = FFTdata2.size

# Adjusting FFT variables
duration = 300 # in seconds
Fs = float(samples)/duration # sampling frequency (sample/sec)
delta_t = 1.0/Fs
t = np.arange(0, samples, 1)*delta_t
FFTdata_freq = np.abs(np.fft.rfft(FFTdata))**1
FFTdata2_freq2 = np.abs(np.fft.rfft(FFTdata2))**1
freq = np.fft.rfftfreq(samples, d=delta_t)
freq2 = np.fft.rfftfreq(samples2, d=delta_t)

# Printing data
plt.semilogy(freq, FFTdata_freq, color='r')
plt.semilogy(freq2, FFTdata2_freq2, color='b')
plt.xticks([0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300, 
320,340,360,380,400,420,440])

第二个代码：

# Number of samplepoints
N = 600
# sample spacing
T = 300.0 / 266336.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))

现在问题来了。我有两个 FFT 代码，但我不知道如何使它们工作。好吧，第一个代码正确地显示了我的数据，在图表上你有上图，但是比例错误，我需要上图位于底部图的位置。我不知道该怎么做。有任何想法吗？ fft1 和 fft2 是数据数组。一切都发生在 300 秒 = 300000 毫秒内。

感谢@zck 更改代码后，它看起来像这样 从 scipy.signal 导入韦尔奇

plt.subplot(212)
plt.title('Fast Fourier Transform')
plt.ylabel('Power [a.u.]')
plt.xlabel('Frequency Hz')
fft1 = (Bx[51:-14])
fft2 = (By[1:-14])

for dataset in [fft1]:
    dataset = np.asarray(dataset)
    psd = np.abs(np.fft.fft(dataset))**2.5
    freq = np.fft.fftfreq(dataset.size, float(300)/dataset.size)
    plt.semilogy(freq[freq>0], psd[freq>0]/dataset.size**2, color='r')

for dataset2 in [fft2]:
    dataset2 = np.asarray(dataset2)
    psd2 = np.abs(np.fft.fft(dataset2))**2.5
    freq2 = np.fft.fftfreq(dataset2.size, float(300)/dataset2.size)
    plt.semilogy(freq2[freq2>0], psd2[freq2>0]/dataset2.size**2, color='b')

我应用了一些更改。我只缺少汉明窗口，任何人都可以帮助，从这张图表中制作：

那个：

Answer 1

快速浏览一下，在上面的代码 sn-p 中，您似乎忘记了除以 N。这是一个数学问题，而不是代码问题。

一般来说，如果您在寻找光谱/功率谱，请使用 WOSA（=重叠段平均）方法，该方法在样本量允许的情况下应用窗口函数和平均值。 welch 方法包含在scipy.signals 中，您应该探索scipy.signal 库，因为它在信号分析中非常方便。

对于您的数据集，以下代码应该可以工作：

from scipy.signal import welch

plt.figure()
for dataset in [Bx, By]:
    dataset = np.asarray(dataset)
    freq, psd = welch(dataset, fs=dataset.size/300, return_onesided=True)
    plt.semilogy(freq, psd/2)

请注意，return_onesided 的 psd 除以 2，如果 False 不需要除以 2。希望这有助于产生好看的图表。上图绘制了功率谱密度。如果您需要功率谱而不是功率谱密度，请传递参数 scaling='spectrum'。

您还可以为窗口函数传递参数，默认为hanning，但它包括最常见的窗口，如 blackman、hamming、boxcart 等。

您可以在https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.welch.html找到更多相关信息