はじめに

　周波数スペクトルを求めたいと思う時、やるべき処理は当然フーリエ変換である。
　だが現実世界でフーリエ変換したいときフーリエ変換したい波形は数式では表現できず離散的な波形データであることが大概である。つまり

　X(ω)=∫x(t) exp{-jωt}dt
　　t: 時間、ω: ω=2πfで表現される角周波数、f: 周波数、j: 虚数単位
　　x(t): フーリエ変換したい信号
　　exp{}: 自然関数
　　X(ω): 信号x(t)のフーリエ変換結果
　　　x(t)にexp{-jωt)をかけたものを-∞から+∞の範囲で時間積分したもの

で示すフーリエ変換は現実で使えることは早々ないということである。
　そこで活用するのが離散フーリエ変換である。

離散フーリエ変換とは

　離散的な波形データをのためのフーリエ変換が離散フーリエ変換である。
　総データ数Nの離散波形x[n]から離散フーリエ変換によって得られる周波数スペクトルX[i]は次式で求められる。

　X[i] = Σx[n]W^(in)
　　W=exp{-j2π/N}
　　{Σの計算範囲は0~N-1のN個)

　離散波形x[n]を通常のフーリエ変換へ代入。離散波形であるが故に積分は総和の計算に置き換わり、積分範囲は-∞~+∞から0~N-1に変わる。
　なお、X[i]は基本角周波数をW0=(2π)/(NTs)とするiW0に対する周波数スペクトルを示している。Tsは離散波形のサンプリング周波数である。計算上はサンプリング周波数Tsが関係ない数式になっているが、周波数スペクトルの周波数軸の算出にはTsを活用することとなる。
　数式中のW^(in)はin=Nで位相が一周するため、in=0~N-1の間で値がループする。つまりW^(in)の取りうる値はNパターン存在する。
　因みにX[i]は見てわかる通り複素数である。我々が一般的に見かける周波数スペクトルはX[i](またはX(ω))の絶対値である。この周波数スペクトルの絶対値に関していえば、
　|X[i]| = |X[N-i]|
が成り立ち、0~N/2の周波数スペクトルとN/2~N-1の周波数スペクトルは同じ形状が折り返したものとなる。つまり有効なスペクトルはN/2までである。

離散フーリエ変換の算出アルゴリズム

　離散フーリエ変換の数式を用いて手計算で周波数スペクトルというのは中々現実的ではない。離散波形はPCなどで扱うデータとしてあるのだから、プログラミング等で離散フーリエ変換するのが最適だろう。
　Pythonを用いたプログラムの一例を以下に示す。
　本プログラムでは\n区切り(改行区切り)で離散波形データを保存したテキストファイルを読み込んで離散フーリエ変換するようにしている。
　算出したスペクトルのインデックスとなる周波数は、入力変数として設けた離散波形のサンプリング周波数から算出する。数式中は角周波数をインデックスとしていたがプログラム中は周波数にしている(ω=2πf)。
　計算結果であるスペクトルは周波数と共にテキストファイルで出力する。エクセルにコピペしてコンマ区切りすれば即座にグラフ化できる状態である。
　計算において離散フーリエ変換式中のW^(in)は複素数であるため、プログラム中では実部と虚部に分けて扱っている。また、Wはinの値によりN通りの値を持つため事前にN通りのリストを算出している。
　本プログラムで算出しているのは周波数スペクトルの絶対値であり、有効なスペクトルのみを算出している。
　離散フーリエ変換の数式を紐解いてそのまま記述しているため難解ではないだろう。

#inport
import math

#入力変数
Input_Pass = 'xxx\WaveData.txt'  #読み込む離散波形ファイルのパス
Output_Pass = 'xxx\SpectrumData.txt'  #出力先ファイルのパス
Fs = 100*1000  #離散波形のサンプリング周波数[Hz]

#波形データ読み込み
print('Data Import')
Data_file = open(Input_Pass,mode='r')  #読み込みモード
Data_origin = Data_file.readlines()  #一行ごとに全データを読み込み
Data_file.close()
N = len(Data_origin)  #データ総数
Data_wave = [0]*N  #離散波形データの定義
for i in range(N):
   Data_wave[i] = float(Data_origin[i].split('\n')[0])  #データ格納

#変数定義
N_list = list(range(N))  #0~N-1の数列
Frequency = [Fs/N*x for x in N_list]  #周波数リスト
W_Re = [math.cos(2*math.pi/N*x) for x in N_list]  #Wの実部リスト
W_Im = [math.sin(2*math.pi/N*x) for x in N_list]  #Wの虚部リスト
Spectrum_abs = [0]  *N#スペクトルの絶対値の定義

#計算
print('DFT Start')
for i  in range(int(N/2)):
   print('DFT ' + str(i+1) + '/' + str(int(N/2)))
   Spectrum_Re = float(0)  #スペクトル実部
   Spectrum_Im = float(0)  #スペクトル虚部
   for n in range(N):
       m = (i * n) % N  #Wに必要な位相を計算
       Spectrum_Re += Data_wave[n] * W_Re[m]
       Spectrum_Im += Data_wave[n] * W_Im[m]
   Spectrum_abs[i] = abs(complex(Spectrum_Re,Spectrum_Im))

#データ出力
print('Data Export')
Spectrum_output = ''  #データ出力用変数
for i in range(int(N/2)):
   Spectrum_output += str(Frequency[i]) + ',' + str(Spectrum_abs[i]) +'\n'
Output = open(Output_Pass,mode='w')  #書き込みモードで読み込む
Output.write(Spectrum_output)
Output.close()

最後に

　離散フーリエ変換のアルゴリズムを考えるところまで行ったが、離散フーリエ変換の数式をそのまま用いるアルゴリズムはかなり計算時間のかかるものとなっている。離散フーリエ変換を高速化したものに高速フーリエ変換(FFT)が存在する。