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

Toshihiko Arai

投稿日: 2021-01-29

はじめに

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

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

そこでおすすめの本をご紹介します。

▼ こちらの書籍はPythonでフーリエ解析を学べます。NumPyやmatplotlibの基本を解説し、SciPyを使った応用プログラムまで網羅してます。 とくに、信号処理への応用が学べるので実践的です。

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

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

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

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

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

こちらがそのプログラムです。

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の正弦波データを作るために、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でフーリエ変換してみました。 https://soundcloud.com/user-999425084/piano-tone-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倍音から大きくゲインが下がっているのも気になります。もしかしてこれはピアノの弦とは別の倍音なのでしょうか？

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

ピアノの弦は、中・高音域では1音に対して弦が3本張られてます。しかし、3本の弦をまったく同じ音程にチューニングするのではないのです。 少しずれてたほうが音に厚みが出るからです。ここらへんは大人数で歌うコーラスの厚みだったり、エフェクターのコーラスを知っている方ならすんなり理解できる話ですね。 さて、9倍音以降が理論値とズレてしまったのは、こういった複雑な絡み合いがあるからなのでしょうか？単純な倍音構成でないからこそ、あのピアノらしい音色を私たちに感じさせてくれるのでしょうかね。ピアノの音色の原理を追求するつもりはないので、この辺で考察を終えて先へ進みます。

逆フーリエ変換して元のピアノ音を作れるか？

「フーリエ変換したデータから逆フーリエ変換して元のピアノ音を作れるだろうか？」 非常に興味のあるテーマでした。さっそくやってみました。逆フーリエ変換した結果がこちらです。

残念ながらというか、予想通りというか、元の波形から随分と異なる波形が生まれました。しかし、こちらの音をぜひ聞いてもらいたいです。逆フーリエ変換したデータをWAVファイルに書き出したものになります。 https://soundcloud.com/user-999425084/piano-tone-a-made-from-inverse-fft

音割れは目をつぶってもらって、前半部分だけ聴くとなんとピアノの雰囲気を感じ取れませんでしょうか？少なくとも音程はバッチリなようです。 この結果には正直驚きました。ぜひもう一度、最初のピアノ音源と比較してみてください。

https://soundcloud.com/user-999425084/piano-tone-a https://soundcloud.com/user-999425084/piano-tone-a-made-from-inverse-fft

そもそも「周波数成分が時間的変化する音源を、単純な逆フーリエ変換で再現することは不可能」と疑ってかかった実験でしたが、予想外の良い結果でした。何でもやってみるものですね！ サンプリング時間をもっと微小時間にしてそれぞれをつなぎ合わせれば、つまり積分のような考え方でつなぎ合わせれば、もっと精度良く原音を復元できそうですね。どなたかやってみますか？

これらのデジタル信号処理を勉強後、デジタルフィルタを使って「波動」というシンセサイザーアプリを作ってみました。 {{apps_oscillator}}

https://www.youtube.com/watch?v=QtlBUK8sYtE&feature=emb_title&ab_channel=ToshihikoArai

関連記事

• 【Python】1/fゆらぎ（ピンクノイズ）のデジタルフィルタ
• フーリエ変換を理解するまでの数学的メモ
• M5StickC PLUSの加速度センサで振動測定と周波数特性
• 【自作マイクへの道④】マイクの風切り音対策「ウインドジャマー」の作り方
• ESP32でESP-NOWを使って通信してみよう