HOME > Programming > Python

Pythonではじめてのフーリエ変換・逆フーリエ変換

2021-09-09　 2021-01-29

Pythonではじめてフーリエ変換・逆フーリエ変換を使ってみましたので、その時の忘備録とともに、これからPythonでFFTを使いたい方がこの記事を役立ていただければ幸いです。

Pythonでフーリエ変換を行う場合は、デジタル信号を扱うことになるため離散フーリエ変換（DFT）のアルゴリズムを使うことになります。これを理解するには非常に難しいです。正直私もDFTはよく分かっていません。(^_^;)

ですが、仕組みがわからなくてもnumpyのfft関数を使えば誰でもフーリエ変換できると思いますので、ぜひチャレンジしてみてください。

▼ ArduinoやRaspberry Piを使えば、FFTを使ったガジェットを作れます！

M5StickC PLUSの加速度センサで振動測定と周波数特性

2021-07-03 ESP32

はじめに

フーリエ変換とは

まず、フーリエ変換は音声データなどで時間領域を周波数領域に変換する数学というかアルゴリズムです。フーリエ変換を行えば、サンプルデータの中にどの周波数がどのくらいの大きさで含まれているかを調べることが簡単にできるようになります。このことを周波数特性と呼びます。

フーリエ変換の知識が無くてもPythonならプログラミングできると思いますが、ある程度の知識は身につけておいたほうが何かにつまづいたときの助けになると思います。

そこでオススメの本をご紹介します。

▼ こちらの書籍は、Pythonでフーリエ解析を学ぶ内容となっています。NumPyやmatplotlibの基本を解説し、SciPyを使った応用プログラムまで網羅しています。

とくに、信号処理への応用が学べるので実践的です。

Pythonで学ぶフーリエ解析と信号処理

Amazonで見る楽天市場で見る

また、こちらの本は、難解なフーリエ変換の数学をとてもわかり易く解説されています。フーリエ変換を理解する導入としてとても役立ちました。

マンガでわかるフーリエ解析

Amazonで見る Kindleで見る楽天市場で見る

上記の本の著者は渋谷道雄氏で、同じ方が書かれたこちらの本も大変読みやすかったです。合わせて読むとフーリエ変換の理解がだいぶ進みますので、ご参考になさってみてください。

Excelで学ぶフーリエ変換

Amazonで見る楽天市場で見る

▼ 上記の2冊を読んで、自分なりにフーリエ変換をまとめてみました。

フーリエ変換を理解するまでの数学的メモ

2019-07-02 数学

フーリエ変換・逆フーリエ変換のプログラム

それではPythonでフーリエ変換・逆フーリエ変換のプログラムを行っていきます。プログラム中でやっていることは次のとおりです。

サンプリング周波数44.1kで1000Hzの正弦波の配列を用意
配列をフーリエ変換して、周波数特性をグラフへプロット
フーリエ変換したデータから逆フーリエ変換して❶の波形を作り出せるかを確認

こちらがそのプログラムとなっています。

import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':
    # 正弦波のデータ作成
    f = 1000
    rate = 44100
    T = np.arange(0, 0.01, 1 / rate)
    s = []
    for t in T:
        v = np.sin(2 * np.pi * f * t)
        s.append(v)

    plt.plot(T, s)
    plt.xlabel('Time')
    plt.ylabel('Gain')
    plt.show()

    # フーリエ変換
    fft_data = np.abs(np.fft.rfft(s))
    freqList = np.fft.rfftfreq(len(s), 1.0 / rate)  # 横軸
    plt.loglog(freqList, 10 * np.log(fft_data))
    plt.xlabel('Frequency')
    plt.ylabel('Power')
    plt.show()

    # 逆フーリエ変換
    r = np.fft.irfft(fft_data, len(T))
    plt.plot(T, r)
    plt.xlabel('Time')
    plt.ylabel('Gain')
    plt.show()

プログラムの解説

先ほどのプログラムの内容を詳しくみていきましょう。

1000Hzの正弦波データ

1000Hzの正弦波

まずは1000Hzの正弦波データを作るために、np.arangeを使って0秒〜0.01秒までの間を、1/44100秒間隔で配列Tを生成します。正弦波の周波数をfとすれば、\(sin(2πft)\)で信号データを生成できますね。よって、下記のプログラムのように書けますよね。

f = 442
rate = 44100
T = np.arange(0, 0.01, 1 / rate)
s = []
for t in T:
    v = np.sin(2 * np.pi * f * t)
    s.append(v)

フーリエ変換で周波数特性

フーリエ変換を行うと実数と虚数のデータの2つが得られますが、虚数データは必要ありません。ですのでnp.fft.rfftを使って実数データのみを取得しています。

また、グラフの横軸を周波数にしたいので、np.fft.fftfreq(サンプル数, サンプリング周期)で周波数軸のリストを作成しました。

そして、音の振幅をdBに直すため10 * np.log(fft_data)で対数変換を行ないます。

fft_data = np.abs(np.fft.rfft(s))
freqList = np.fft.rfftfreq(len(s), 1.0 / rate)  # 横軸

こちらのグラフのように、1kHzの周波数成分だけを含むことが確認できました。元の信号の周波数と一致しましたね。

フーリエ変換で周波数特性

逆フーリエ変換で元の波形を再現

ここでは、フーリエ変換されたデータを元に逆フーリエ変換を行っています。

r = np.fft.irfft(fft_data, len(T))

すると元のデータと同じ数の山を持った、つまり同じ周波数で同じ振幅の波形を作り出すことができました。

逆フーリエ変換から元の波形を再現

ただし、位相は考慮されていないため余弦波となってしまいました。フーリエ変換では周波数位相特性も調べることができるので、もしかしたらそれを使えば元の位相と同じ波形を作れるかもしれませんね。

ところで、人間の耳は音の位相に鈍感です。余弦波も正弦波も同じ周波数であれば区別が付きません。同じ音色に聞こえても、まったく違う波形になることもあるのです。

ピアノの音をフーリエ変換してみた

ここからは、興味本位でさらに実験を進めてみました。

次のリンクにありますように、ピアノのA音をwavファイルで用意し、Pythonでフーリエ変換してみました。

基音Aのピアノ音

import numpy as np
import matplotlib.pyplot as plt
import wave
from scipy.io.wavfile import write


def load_wav(filename, N):
    with wave.open(filename, 'r') as wav:
        fr = wav.getframerate()
        data = wav.readframes(N)
        s = np.frombuffer(data, dtype=np.int16)
        return s, fr


def output_wav(filename, rate, data):
    write(filename, rate, data.astype(np.int16))


if __name__ == '__main__':
    N = 44100
    s, rate = load_wav('piano_a.wav', N)
    T = np.arange(0, N / rate, 1 / rate)

    plt.plot(T, s)
    plt.xlabel('Time')
    plt.ylabel('Gain')
    plt.show()

    # フーリエ変換
    fft_data = np.abs(np.fft.rfft(s))
    freqList = np.fft.rfftfreq(len(s), 1.0 / rate)  # 横軸
    plt.loglog(freqList, 10 * np.log(fft_data), '.-')
    plt.xlabel('Frequency')
    plt.ylabel('Power')
    plt.show()

    # 逆フーリエ変換
    r = np.fft.irfft(fft_data, len(T))
    plt.plot(T, r)
    plt.xlabel('Time')
    plt.ylabel('Gain')
    plt.show()
    output_wav('piano_a_irfft.wav', rate, r)

ピアノ音の波形

プログラムを実行してプロットしたピアノ音は、図のような波形になります。

ピアノ音をフーリエ変換

ピアノ音をフーリエ変換して周波数特性を調べてみます。ただし、音源の1秒間だけをサンプルしてフーリエ変換を行っています。

周波数特性をグラフにすると次のとおりです。

ピアノ音の周波数特性

ここで、基音のA音が440Hzにあり、そこから倍音が2倍、3倍と続いているのがグラフから確認できます。また、8倍音あたりまでは計算通りの周波数で倍音があらわれますが、9倍音以降になると大きくズレてきます。9倍音から大きくゲインが下がっているのも気になります。もしかしてこれはピアノの弦とは別の倍音なのでしょうか？

先ほどの周波数特性をグラフにして整理してみます。