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

 

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

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

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

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

はじめに

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

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

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

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

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

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

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

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

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

それでは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の正弦波
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)で対数変換を行ないます。

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

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

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

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

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

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

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

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

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

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

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

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

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

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倍音から大きくゲインが下がっているのも気になります。もしかしてこれはピアノの弦とは別の倍音なのでしょうか?

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

倍音理論値(Hz)目分量で実測(Hz)
1440440
2880880
313201308
417601784
522002200
626402695
730803138
835203563
939604122
1044004650
1148405115
1252805956
1357206634

ピアノの弦
ピアノの弦

ピアノの弦は、中・高音域では1音に対して弦が3本張られています。しかも、3本をまったく同じ音程にチューニングするのではないようです。若干の音程差をつけて音に厚みを出すためです。大人数で歌うコーラスの厚みだったり、エフェクターのコーラスを知っている方は理解できると思います。

9倍音以降が理論値とズレているのは、こういった複雑な絡み合いがあるからなのでしょうか?単純な倍音構成でないからこそ、あのピアノらしい音色を私たちに感じさせてくれるのでしょうね。

さて、ピアノの音色の原理を追求するつもりはないので、この辺で考察を終え、先へ進みましょう。

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

「フーリエ変換したデータから逆フーリエ変換して元のピアノ音を作れるだろうか?」

非常に興味のあるテーマでした。さっそくやってみました。逆フーリエ変換した結果がこちらです。

ピアノ音の周波数特性から逆フーリエ変換
ピアノ音の周波数特性から逆フーリエ変換

残念ながらというか、予想通りというか、元の波形から随分と異なる波形が生まれました。しかしながら、こちらの音をぜひ聞いてもらいたいです。逆フーリエ変換したデータをWAVファイルに書き出したものになります。

音割れは目をつぶってもらって、前半部分だけ聴くとなんとピアノの雰囲気を感じ取れませんでしょうか?少なくとも音程はバッチリなようです。

この結果には正直驚きました。ぜひもう一度、最初のピアノ音源と比較してみてください。

そもそも「周波数成分が時間的変化する音源を、単純な逆フーリエ変換で再現することは不可能」と疑ってかかった実験でした。何でもやってみるものですね!

サンプリング時間をもっと微小時間にしてそれぞれをつなぎ合わせれば、つまり積分のような考え方でつなぎ合わせれば原音を復元できそうですね。

どなたかやってみませんか?笑

▼ デジタルの信号処理を詳しく勉強したい方はこちらの書籍がオススメです。高速フーリエ変換(FFT)や離散フーリエ変換(DFT)の仕組みから、IIRやFIRなどのデジタルフィルタが学べます。数学やC言語がメインなので難易度は高めですが、サンプルプログラムが豊富で実践しながら勉強できます。

これらのデジタル信号処理を勉強後、デジタルフィルタを使って「波動」というシンセサイザーアプリを作ってみました。かなりマニアックな内容のアプリですが(^_^;)

興味のある方はこちらのYouTube動画をご覧ください。

記事に関するご質問などがあればTwitterへお返事ください。
この記事で紹介した商品
Python学習にオススメの本
関連記事