Pythonで高速フーリエ変換を実装

前書き

音声処理でよくでてくる高速フーリエ変換ですが、大体の人がnumpyを使ってると思います。どんなアルゴリズムなのか気になったのでPythonで雑に実装してみました。

高校時代に書いたので間違ってたらすまん。

コード

今回は振幅100周波数300、振幅200周波数100、振幅500周波数240の正弦波を合成したものをdataとして実験してみます。

#-*- coding:utf-8 -*-
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import math
import numpy as np
import cmath
import time

class FFT():

    def __init__(self):
        self.N = 1024
        self.n = int(math.log(self.N)/math.log(2)+0.1)
        self.data = np.zeros(self.N,dtype=np.complex)
        self.datar = np.zeros(self.N,dtype=np.complex)
        self.datai = np.zeros(self.N,dtype=np.complex)
        self.tempr = np.zeros(self.N,dtype=np.complex)
        self.tempi = np.zeros(self.N,dtype=np.complex)
        self.dtemp = np.zeros(self.N,dtype=np.complex)
        self.temp = 0
        self.d = 0
        self.a = 1
        self.AMP = np.zeros(self.N,dtype=np.float32)
        for i in range(self.N):
            self.data[i] = 100*math.sin(300*2* np.pi/self.N*i)+200*math.sin(100*2 * np.pi/self.N*i)+500*math.sin(240*2 * np.pi/self.N*i)
        self.bitReverse()
        self.butterfly()


    #ビット反転
    def bitReverse(self):
        for j in range(1,self.N):
            self.d = 0
            self.a = 1
            for i in range(self.n):
                if (self.a & j) != 0:
                    self.d += int(self.N /(2**(i+1)))
                self.a*=2
            if(j < self.d):
                self.temp = self.data[j]
                self.data[j] = self.data[self.d]
                self.data[self.d] = self.temp

    #バタフライ演算
    def butterfly(self):
        #ステージ
        for i in range(self.n):
            #子ブロック
            for m in range(int(self.N/(2**(i+1)))):
                #ブロック内順番
                for k in range(2**(i+1)):
                    W = np.exp(-1j*2*np.pi*k/(2**(i+1)))
                    if k < 2**(i+1)/2:
                        self.dtemp[int(k+2**(i+1)*m)] = self.data[int(k+2**(i+1)*m)]+self.data[int(k+2**(i+1)*m+2**i)]*W
                    else:
                        self.dtemp[int(k+2**(i+1)*m)] = self.data[int(k+2**(i+1)*m)]*W+self.data[int(k+2**(i+1)*m-2**i)]
            self.data = self.dtemp
            self.dtemp = np.zeros(self.N,dtype=np.complex)
        for i in range(self.N):
            self.AMP[i] = np.sqrt(self.data[i].real**2+self.data[i].imag**2)
            print(self.AMP[i])

if __name__ == '__main__':
    fft = FFT()

ちょっとした解説

高速フーリエ変換の肝はビット反転とバタフライ演算です。

• ビット反転

ビット列を反対に入れ替える処理です。ex)1011→1101
各ビットに1とand演算をしてビットの存在を確認して反転させます。その後2倍して左にシフトしています。

• バタフライ演算

有名なあのフローを素直にプログラムに起こしただけです。

参照<高速フーリエ変換 - Wikipedia>

蝶形のフローの各列をステージとし、行をブロック、その中の計算順序をブロック内順番としています。

実行結果

出力結果を適当なファイルにリダイレクトしてexelでみてみましょう。

結果

フーリエ変換の対称性がちゃんと現れています。

感想

高校レベルの知識でも高速フーリエ変換を実装(理解したとは言っていない)できることがわかりました。

Cythonとかを使って高速化してもnumpyには勝てないので
みんなもnumpyを使おう!

参考文献

www.amazon.co.jp