感動！数学って素敵やん！　〜離散フーリエ変換したら芸術的な絵が出てきた！〜

どうもこんにちは、

前回の記事で離散フーリエ変換をしていた時に、フーリエ変換がどうなっているのかを途中の結果を抜き出してプロットしてみたら芸術的な絵が出てきたので、みんなと共有したく、ブログを書き始めました。

また、芸術的な絵が出てきたのは離散フーリエ変換という美しい数式を作ってくれた人がいるからであり、なぜか一瞬作った人に感謝を伝えたくなりました（笑）。

では、まず前回の復習ということで一応離散フーリエ変換の説明をのしてから私が感動したグラフの紹介に移りたいと思います。

離散フーリエ変換

アナログ信号のフーリエ変換では、x(t)とX(f)のどちらも無限長の連続信号となりますが、コンピュータでは無限長の信号を扱うことが出来ません。そのため、コンピュータを使ってフーリエ変換を計算するには、tとfをそれぞれ０からN-1までの整数nとkに置き換えた「離散フーリエ変換（DFT）」を用います。

ここで、t_sは標本化周期、f_saは標本化周波数とします。

$t = nt_s \ (0 \le n \le N-1)} \tag{3$

$f = \frac{kf_s}{N} \ (0 \le k \le N-1)} \tag{4$

DFT（離散フーリエ変換）とIDFT（逆離散フーリエ変換）：

$X(k) = \sum_{n=0}^{N-1} x(n)exp(\frac{-j2πkn}{N} ) \ (0 \le k \le N-1)} \tag{5$

$x(n) = \sum_{n=0}^{N-1} X(k)exp(\frac{j2πkn}{N} ) \ (0 \le n \le N-1)} \tag{5$

ここで、x(n)は時間nを変数とするディジタル信号、X(k) は周波数kを変数とするx(n)の周波数特性を表します。

音声データの場合、x(n)は実数となり、X(k)は一般的に複素数となります。

感動した実験

以下のサイン波をDFTしてみました。

振幅：１

サンプリング周波数：１００

周波数：１０

コード


def dft(sig):
    X = []
    N = len(sig)
    tmpReal = []
    tmpImag = []
    for k in range(N):
        w = np.exp(-1j * 2 * np.pi * k / N)
        X_k = 0
        print(k)
        tmpReal = []
        tmpImag = []
        for n in range(N):
            X_k += sig[n] * (w**n)
            tmpReal.append(X_k.real)
            tmpImag.append(X_k.imag)
            
        plt.plot(tmpReal, tmpImag)
        plt.show()
        X.append(abs(X_k))
    return X

変数：

　X:　dftの結果を格納する

　N:　入力データの長さ

for文（dft処理部）:

１つ目のforで各周波数の計算を行っていく。

そして２つ目のforで、ある周波数kにおけるdftを音データ全てにかける。

そして、前回と違うのが、変数のtmpReal, tmpImagです。

この２つの変数を使って、各周波数kのdftを実部と虚部に分けてプロットしました。

つまり、以下の式を実部と虚部に分けたということです。

${\sum_{n=0}^{N-1} x(n)exp(\frac{-j2πkn}{N} ) \ (0 \le k \le N-1)}$