ガレスタさんのDIY日記

電子回路、Web、組み込み、自作エフェクターを語るblog

Python:ScipyのFFT(scipy.fftpack)をやってみる。

PythonFFTをする記事です。


FFTは下に示すように信号を周波数スペクトルで表すことができどの周波数をどの程度含んでいるか可視化することができます。

440Hzの場合

f:id:gsmcustomeffects:20180810010502p:plain

2000Hzの場合

f:id:gsmcustomeffects:20180810010524p:plain

コード

numpyとScipy両方に同じようなメソッドがあるけどScipyおじさんなのでscipy.fftpackを使います。

from pylab import *
import numpy as np
import scipy.fftpack as spfft



f0 = 440
fs = 96000
N = 1000
addnum = 5.0

def create_sin_wave(amplitude,f0,fs,sample):
    wave_table = []
    for n in np.arange(sample):
        sine = amplitude * np.sin(2.0 * np.pi * f0 * n / fs)
        wave_table.append(sine)
    return wave_table

wave1 = create_sin_wave(1.0,f0,fs,N)

X = spfft.fft(wave1[0:N])
freqList = spfft.fftfreq(N, d=1.0/ fs)

amplitude = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X]  # 振幅スペクトル

# 波形を描画
subplot(211) 
subplots_adjust(hspace=0.5)
plot(range(0,N), wave1[0:N],label = "wave1")
axis([0, 500, -1.0, 1.0])
xlabel("sample")
ylabel("amplitude")

# 振幅スペクトルを描画
subplot(212)
plot(freqList, amplitude, marker='.', linestyle='-',label = "fft plot")
axis([0, fs / 6, 0, 1000])
xlabel("frequency [Hz]")
ylabel("amplitude")

show()

注意点など

プログラムとしては適当なサンプル分正弦波を作成してそれをfftメソッドに投げるだけ。
ここで注意してほしいのがFFTして帰ってくる値は複素数を含むということ(complex128型)
f:id:gsmcustomeffects:20180810010719p:plain

値の確認はPycharmのデバッガからできる(はじめてつかうライブラリの挙動確認にはほんとに便利です