代码之家 › 专栏 › 技术社区 › Hiddenguy

Python 2.7中的FFT合并两个代码

merging-data fft numpy python-2.7 python

Hiddenguy · 技术社区 · 8 年前

我有两个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秒=30万毫秒。

在使用@zck更改代码后,看起来是这样的来自scipy。信号导入welch

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

我应用了一些更改。我只是错过了汉明窗口,有谁能帮我从这张图表中得出:

那一个:

1 回复 | 直到 8 年前

zck 8 年前

从快速查看中可以看出,在上面的代码段中,您忘记了除以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)

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

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

你可以在以下网站找到更多信息: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.welch.html