numpy.fftのテンプレート

matplotlibのテンプレートにつづいて、numpyでFFTするときのテンプレートです。1.5 MHzのsin関数を5.8MHz 4ビット分解能のADCで5000点サンプルしたときの結果です。np.fft.fftfreq()関数で周波数軸も求められるので、Hzで表示しています。

GistにはJupyter Notebookをそのままアップロードしておいたので、プロット結果、sin.textの中身とあわせて見ることができます。

import matplotlib.pyplot as plt
import numpy as np

# Constants
dataFile = "sin.text"
samplingFreq = 5.8e6
samplingPeriod = 1/samplingFreq

# Read data
# 数字が文字として並んだテキストファイルを読み込む。
# sin.textは、sinを4ビットでサンプルした値が5000個並んだデータ。
dataString = open(dataFile).read()
n = len(dataString)
d = np.zeros(n)
for i in np.arange(n-1):
    d[i]=int(dataString[i])

# Subtract the DC component
d = d - np.mean(d)

# FFT - X axis
fftfreq = np.fft.fftfreq(n, samplingPeriod)

# FFT - Y axis
fft = np.fft.fft(d)[0:n/2]
fftmag = np.abs(fft)

# FFT - Plot
fig, panel = plt.subplots(1,1)
panel.set_xlim(0, fftfreq[:len(fftmag)][-1])
# panel.set_ylim(1e-1,1e4)
# panel.set_yscale("log")
panel.plot(fftfreq[:len(fftmag)], fftmag)

# Save figure
fig.savefig("fft_sample.png", dpi=200)

# Maximum power bin and frequency
binIndexMaxPower = np.argmax(fftmag)
freqResolution = samplingFreq / n
freqMaxPower = binIndexMaxPower*freqResolution
print "Maximum-power bin       : {:d}".format(binIndexMaxPower)
print "Maximum-power frequency : {:.3f}".format(freqMaxPower)