FFT in Python 2.7合并两个代码

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')

我应用了一些更改。我只缺少 Hamming Window，任何人都可以帮助从这张图表中得出：

那个：

Answer 1

乍一看，上面的代码片段似乎忘记除以 N。这是一个数学问题，而不是代码问题。

一般来说，如果您在 spectra/power 光谱之后使用 WOSA（=重叠段平均）方法，如果样本量允许，该方法应用 window 函数和平均值。 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)

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

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

您可以在 https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.welch.html

找到更多相关信息

FFT in Python 2.7合并两个代码

FFT in Python 2.7 merging two codes

python

numpy

fft

merging-data

python-2.7