python中的FFT谱图

Question

我有 python 3.4。

我发送了一个 2MHz（例如）频率并随着时间的推移接收到空化（直到我停止测量）。我想得到一个频谱图（空化与频率），更有趣的是在次谐波 (1MHz) 频率时间内的空化频谱图。

数据保存在sdataA（=空化）和t（=测量时间）

我试图在 FFTA 中保存 fft

FFTA = np.array([])
FFTA = np.fft.fft(dataA)
FFTA = np.append(FFTA, dataA)

我得到了实数和复数然后我只取了一半（从 0 到 1MHz）并保存了实数和复数数据。

nA = int(len(FFTA)/2)
yAre = FFTA[range(nA)].real
yAim = FFTA[range(nA)].imag

我试图通过以下方式获取频率：

FFTAfreqs = np.fft.fftfreq(len(yAre))

但这是完全错误的（我打印了数据print (FFTAfreqs)）

我还绘制了数据，但它又是错误的：

plt.plot(t, FFTA[range(n)].real, 'b-', t, FFTA[range(n)].imag, 'r--')
plt.legend(('real', 'imaginary'))
plt.show()

如何在次谐波 (1MHz) 频率范围内输出空化频谱图？

编辑：

数据示例：

查看“dataA”和“time”的示例：

dataA = [6.08E-04,2.78E-04,3.64E-04,3.64E-04,4.37E-04,4.09E-04,4.49E-04,4.09E-04,3.52E-04,3.24E-04,3.92E-04,3.24E-04,2.67E-04,3.24E-04,2.95E-04,2.95E-04,4.94E-04,4.09E-04,3.64E-04,3.07E-04]
time = [0.00E+00,4.96E-07,9.92E-07,1.49E-06,1.98E-06,2.48E-06,2.98E-06,3.47E-06,3.97E-06,4.46E-06,4.96E-06,5.46E-06,5.95E-06,6.45E-06,6.94E-06,7.44E-06,7.94E-06,8.43E-06,8.93E-06,9.42E-06]

编辑二： 从@Martin 示例中，我尝试了以下代码，如果我做得对，请告诉我。

如果 dataA 和 Time 保存为 h5 文件（或我已经发布的数据）

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

dfdata = pd.read_hdf("C:\\data_python\\DataA.h5")
dft = pd.read_hdf("C:\\data_python\\time.h5")
dft_cor = int((len(dft)-2)*4.96E-6)   # calculating the measured time

fs = 2000000 #sampling frequency 2MHz
CHUNK = 10000
signal_time = dft_cor # seconds

def sine(freq,fs,secs):
    data=dfdata
    wave = np.sin(freq*2*np.pi*data)
    return wave

a1 = sine(fs,fs,120)
a2 = sine(fs/2,fs,120)

signal = a1+a2
afft = np.abs(np.fft.fft(signal[0:CHUNK]))
freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
spectrogram_chunk = freqs/np.amax(freqs*1.0)

# Plot spectral analysis
plt.plot(freqs[0:1000000],afft[0:1000000]) # 0-1MHz
plt.show()

number_of_chunks = 1000

# Empty spectrogram
Spectrogram = np.zeros(shape = [CHUNK,number_of_chunks])

for i in range(number_of_chunks):
    afft = np.abs(np.fft.fft(signal[i*CHUNK:(1+i)*CHUNK]))
    freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
    spectrogram_chunk = afft/np.amax(afft*1.0)
    try:
        Spectrogram[:,i]=spectrogram_chunk
    except:
        break


import cv2
Spectrogram = Spectrogram[0:1000000,:]
cv2.imshow('spectrogram',np.uint8(255*Spectrogram/np.amax(Spectrogram)))
cv2.waitKey()
cv2.destroyAllWindows()

Answer 1

看来您的问题不在于 Python，而在于理解什么是 Spectrogram。

频谱图是信号的频谱分析序列。

1) 您需要在 CHUNKS 中切断您的信号。

2) 对这些 CHUNKS 进行光谱分析并将它们粘在一起。

例子：

您有 1 秒的音频录制时间（44100 HZ 采样）。这意味着录音将有 1s * 44100 -> 44100 个样本。您定义 CHUNK 大小 = 1024（例如）。

对于每个块，您将进行 FFT，并将其粘贴到 2D 矩阵中（X 轴 - CHUNK 的 FFT，Y 轴 - CHUNK 编号，）。44100个样本/CHUNK~44个FFT，每个FFT覆盖1024/44100~0.023秒的信号

CHUNK 越大，频谱图越准确，但“实时性”越低。

CHUNK 越小，频谱图就越不准确，但是当您“更频繁地”测量频率时，您有更多的测量结果。

如果你需要 1MHZ——实际上你不能使用高于 1MHZ 的任何东西，你只需要得到 FFT 阵列的一半——而且哪一半并不重要，因为 1MHZ 只是你采样频率的一半，而 FFT 正在镜像任何高于采样频率的 1/2。

关于 FFT，您不需要复数。你想做

FFT = np.abs(FFT) # 编辑 - 我刚刚注意到你使用了“.real”，但我会保留在这里

因为你想要实数。

频谱图的准备- 频谱图示例

具有 150HZ 波和 300HZ 波的音频信号

import numpy as np
import matplotlib.pyplot as plt


fs = 44100#sampling frequency
CHUNK = 10000
signal_time = 20 # seconds

def sine(freq,fs,secs):
    data=np.arange(fs*secs)/(fs*1.0)
    wave = np.sin(freq*2*np.pi*data)
    return wave

a1 = sine(150,fs,120)
a2 = sine(300,fs,120)

signal = a1+a2
afft = np.abs(np.fft.fft(signal[0:CHUNK]))
freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
spectrogram_chunk = freqs/np.amax(freqs*1.0)

# Plot spectral analysis
plt.plot(freqs[0:250],afft[0:250])
plt.show()

number_of_chunks = 1000

# Empty spectrogram
Spectrogram = np.zeros(shape = [CHUNK,number_of_chunks])

for i in range(number_of_chunks):
    afft = np.abs(np.fft.fft(signal[i*CHUNK:(1+i)*CHUNK]))
    freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
    #plt.plot(spectrogram_chunk[0:250],afft[0:250])
    #plt.show()
    spectrogram_chunk = afft/np.amax(afft*1.0)
    #print(signal[i*CHUNK:(1+i)*CHUNK].shape)
    try:
        Spectrogram[:,i]=spectrogram_chunk
    except:
        break


import cv2
Spectrogram = Spectrogram[0:250,:]
cv2.imshow('spectrogram',np.uint8(255*Spectrogram/np.amax(Spectrogram)))
cv2.waitKey()
cv2.destroyAllWindows()

单个 CHUNK 的光谱分析

频谱图