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

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

Pythonではじめてフーリエ変換・逆フーリエ変換を使ってみた。

フーリエ変換とは、音声データなどの時間領域を周波数領域に変換する数学だ。フーリエ変換を使えば、データの中にどの周波数がどのくらいの大きさで含まれているかを調べることができる。つまり、周波数特性を調べることができるのだ。

Pythonでフーリエ変換を行う場合は、デジタル信号を扱うことになるため離散フーリエ変換(DFT)のアルゴリズムを使うことになるが、仕組みがわからなくてもnumpyfft関数を使えば誰でも簡単にフーリエ変換できるので安心を。

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

今回実行したプログラムの手順としては次の通り。

  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の正弦波

np.arangeを使って0秒〜0.01秒までを1/44.1k秒の間隔で配列を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)

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

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

フーリエ変換では実数・虚数のデータが得られるが、虚数データは必要ないため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)

ピアノ音の波形

ピアノ音の波形
ピアノ音の波形

ピアノ音のwavデータはこんな感じの波形になっている。プログラムでは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本張られている。すべて同じ音程にチューニングしているのではなく、若干の音程差をつけて音に厚みを増しているのだ。ここら辺は、大人数で歌うコーラスの厚みだったり、エフェクターのコーラスを想像すると納得できる。

そのような理由から、9倍音以降が理論値とズレているのかもしれない。そして単純な倍音構成でないからこそ、あのピアノらしい音色を私たちに感じさせてくれるのだろう。

ところで、周波数特性の図で低域が大分持ち上がってしまっているが、これはピアノの打音(アタック音)がFFTの対象に含まれているからだと思う。

ピアノの音色の原理を追求するつもりはないので、今回はこの程度の考察で先へ進むことにしよう。

逆フーリエ変換して元の波形が作れるか?

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

さて、フーリエ変換したデータから逆フーリエ変換して元のピアノ音が作れるだろうか。そう思って実際に行ってみた。グラフの通りあまりにも元の波形からかけ離れた形となってしまった。

しかしながら、音割れを無視して前半部分だけ聴くとピアノの雰囲気を感じ取れるではないか。これには正直驚いた。しかも音程はバッチリだ。

実際にwavデータに書き出したので、先ほどのピアノ音と比較して聴いてもらいたい。

周波数成分が時間的変化する音源を、単純な逆フーリエ変換で再現することは不可能だと思っていたが実験してみるものである。

おすすめの書籍

こちらの二冊は、フーリエ変換の仕組みを理解した方におすすめの書籍である。同じ著者が書いた書籍で、私の場合は2つ読むことでフーリエ変換の理解がだいぶ進んだ。

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

身近な具体例でフーリエ解析のイメージがつかめる!女子高生3人組のバンド結成の物語を絡めてマンガでフーリエ解析をわかりやすく解説した。フーリエ解析に必要な数学の基礎(微分や積分、三角関数)から解説しており、フーリエ解析が苦手な人も読むことができる。フーリエ解析は必要だけど専門書を読むの苦手な人、電気・情報系の高専・学生、電気・情報系の初級技術者が導入書として読む最初一冊としてお勧め。

KindleAmazon
Excelで学ぶフーリエ変換
Excelで学ぶフーリエ変換

信号を解析する際,なくてはならないのがフーリエ変換だ。しかし,いざ勉強しようと思っても,一から勉強するには敷居が高い。微積分や三角関数など多くの数学的な知識が必要だからだ。また,数式で結果を導けたとしても,イメージが湧かず,なかなか理解できない。このギャップを埋めるのが本書だ。Excelを使い,微積分/三角関数の基礎からフーリエ変換によるスペクトル解析まで解説している。学生時代,専門的に勉強する機会を得られなかった読者にも分かりやすい教材だろう。

Amazon

また、高速フーリエ変換(FFT)や離散フーリエ変換(DFT)のアルゴリズムを学びたい方はこちらの本をおすすめしておく。

C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理
C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理

Amazon

こちらの記事も参考に。

最後まで読んでいただきありがとうございました。

「この記事が参考になったよ」という方は、ぜひ記事をシェアをしていただけるととても嬉しいです。

今後も有益な記事を書くモチベーションにつながりますので、どうかよろしくお願いいたします。↓↓↓↓↓↓↓

あなたにおすすめ