Pythonで検証しながらSwiftで離散フーリエ変換を実装する

  • このエントリーをはてなブックマークに追加

離散フーリエ変換は単純な式で実装できる信号処理アルゴリズムです。Pythonは豊富なライブラリについ頼ってしまいますが、他の言語でも自由に実装できるように基本を学びましょう。ここではPythonであえてforループを使った離散フーリエ変換の方法を学び、次に同様の結果が得られるSwiftコードを紹介します。

こんにちは。wat(@watlablog)です。今回は基本の離散フーリエ変換式をPythonで確認しながら、iOSでも使えるようSwiftでも実装してみます

本記事の目標とモチベーション

Swiftで離散フーリエ変換を実装する

WATLABブログでは最近iOS版の信号処理アプリ制作を目標にしています。「SwiftUIとAVFoundationで音声を再生する方法」で紹介した通り、Swiftで録音したデータを再生したりグラフに波形を表示させたりできるようになってきました。
そのため次は周波数分析を実装するためにフーリエ変換を実装することを目標とします。

ただし、離散フーリエ変換DFT:Discrete Fourier Transform)は録音波形等の大規模サンプリングデータを扱うにはあまりおすすめできないことが注意点です。今回はあくまで感覚を掴むために、離散フーリエ変換のアルゴリズムを学ぶ記事という位置付けです。

より実用的な周波数分析には高速フーリエ変換FFT:Fast Fourier Transform)があります。Swiftにもライブラリがあるとのことなので、この記事はその前段とします(あとで書く!)。

Pythonで検証する理由

最終目標はSwiftによる離散フーリエ変換の実装ですが、今回は最初にPythonでアルゴリズムを実装します。

あえて一度Pythonで実装するのはPythonが数式をほぼそのままの形で直感的に書けたり、プログラミング特有の知識がほぼいらないからアルゴリズムのデバッグに向いているというのが理由です。上記の特徴は本職のプログラマからすると色々賛否両論ありそうですが、簡単にグラフが書けたり数式を実装できるのはプロトタイプが作りやすく非常にありがたいです。

Pythonで離散フーリエ変換をフルスクラッチするコード

この記事は基本的なアルゴリズムを学ぶための位置付けであり、計算自体は遅いコードを紹介しています。Pythonによる実用的なDFTやFFTコードはSciPyを使った方法として以下の記事にまとめていますので、是非こちらもご覧ください。
PythonでFFT実装!SciPyのフーリエ変換まとめ

離散フーリエ変換の式

離散フーリエ変換を式(1)に示します。ここで\(F\)は時間信号\(x\)をフーリエ変換した後のフーリエスペクトルです。\(k\)は周波数インデックス、\(n\)は標本点(時間信号のインデックス)を示し共に整数です。\(N\)は任意の自然数で通常時間信号\(x\)の総数を指します。

そして虚数\(j\)を使って指数部が構成されており、フーリエスペクトル\(F\)は複素数となります。

\[ F[k] = \sum_{n=0}^{N-1} x[n] e^{ (-j \frac{2 \pi k n}{N})} \tag{1} \]

周波数軸は式(2)で計算されます。

\[ f = \frac{k}{N \cdot dt} \tag{2} \]

式(1)による離散フーリエ変換の計算量は\(O(N^{2})\)とデータ数の二乗になります。これは録音データ処理にはなかなか適用が難しいレベルの計算量です。

PythonによるDFTコード

DFT関数部

まずは以下にdft関数を示します。この関数は時間軸t時間信号値xを引数として、離散フーリエ変換後の周波数軸freq振幅成分ampを返すものです。

DFTをすると信号が正の周波数軸と負の周波数軸に分散されます。振幅成分を物理的に正確になるよう求めるためにはNによるスケーリングが必要です。振幅成分のスケーリングは正と負の周波数に分散された領域では\(\frac{N}{2}\)で除算しますが、直流成分である0[Hz]のデータは\(N\)で除算します(この部分はなかなかネットに説明がなく、今回ようやく理解したところです…)。

そのため関数内ではif文を用いてスケールファクタsfを分岐させています。

NumPyにはもっと効率的な書き方がありますが、今回はアルゴリズムの理解のためにforループを使って実装してみました。

全コード

以下に全コードを示します。matplotlibで簡単にグラフ化してみましょう。

実行結果の例

こちらが実行結果の例です。上記プログラムを実行すると、サンプル信号として振動数1[Hz]の正弦波を生成します。このサンプル信号は直流成分の検証を行うためにdcを加算しています。
下図の周波数波形を確認すると、0Hzと1Hzの振幅成分はプログラムで指定した値になっていることがわかりました。

実行結果の例(サインカーブとその離散フーリエ変換結果)

if文でスケーリング係数を変更しないと直流成分が合わない…といったことになり得るためご注意ください。

Swiftで離散フーリエ変換をフルスクラッチするコード

基本のコード:正弦波を表示させる

これからSwiftで離散フーリエ変換のコードを実装しますが、まずは基本とするコードを紹介します。
以下のコードはButtonを押すとwavePoints関数で正弦波を生成し、Chart{}でグラフ表示するだけのプログラムです。以下の記事を参考にしています。
Swift入門:最低限覚えておく基礎文法の備忘録

コードを実行してButtonをクリックすると以下の画面になります。正弦波が表示されました。

Swiftによる正弦波表示

複素数型の構造体をつくる

離散フーリエ変換は複素数Complex)を扱いますが、Swiftには標準で複素数型がありません。そのため自分で定義する必要があります。

Swiftにおいてユーザー定義による構造体はstructで作成します。以下のコードはSwiftで複素数を扱えるように定義した構造体です。

static func +」と「static func *」はそれぞれ加算と乗算の演算子が使われた場合の処理を意味しています。ここでそれぞれのstatic funcには「(lhs: Complex, rhs: Complex)」と書かれていますが、これは複素数が演算子の左側(Left hand side)と右側(Right hand side)に配置された場合の処理を指定しています。

離散フーリエ変換の式(1)では複素数同士の加算と乗算が使われているので、今回は最低限この部分の記述をしているというわけです。

複素数同士の加算(a+jb)+(c+jd)は実部がa+c, 虚部が(b+d)jです。Swiftでは標準で複素数型を持っていないため、この構造体の中で実部と虚部に分けて整理する…という処理をしているのですね。

仮にFloat型とDouble型の除算を定義したい時は、「static func /(lhs: Float, rhs: Double)」と書けば良いというわけです。ややこしいですが、わかれば自分で自由に構造体をつくれそうです。

Swiftによる離散フーリエ変換の関数

続いて本題の離散フーリエ変換の関数を以下に示します。何度も書きますが、この関数はアルゴリズムを学ぶための自作コードなので非常に遅いことに注意してください。別記事でライブラリを使った場合も学びます。

先ほどのPythonコードとほぼ対応させているのでコメントアウトの説明を見ていただければ内容がわかると思います。

Pythonは複素数を標準で扱えるため複素指数関数としてネイピア数の冪計算をしていましたが、Swiftでは先の理由から複素指数関数の構造体を作るか工夫しなければなりません。最も簡単なのはオイラーの公式複素指数関数を実数の余弦関数と虚数の正弦関数に変換する方法です。
先ほどのコードは式(3)を使っています。

\[ e^{jx} = \cos x + j \sin x \tag{3} \]

全コード例

コピペ可能な全コードを以下に示します。

こちらがDFT計算結果です。Pythonの結果と同様に直流成分と変動成分の振幅が1になっており、さらに鏡像ピークも確認できました。

SwiftによるDFT結果

指定した1[Hz]のピークが確かに立っていることも拡大することで確認可能です。

SwiftによるDFT結果2

拡大は通常Chartでやるのが良いかもしれませんが、今回はデータそのものをスライスしてみました。以下のコードでデータ範囲を絞ることができます。

まとめ

この記事ではPythonで離散フーリエ変換のアルゴリズムを確認した後、Swiftで同様の機能を有するコードを書いてみました。

今回紹介したコードはアルゴリズムの理解を優先したためまだ実用からは少し遠いものとなっていますが、関数や複素数構造体の定義、Chartへ周波数データを渡すといった一連の作業を体験できました。

SwiftにはAccelerateといったライブラリを使うことで高速フーリエ変換ができるらしいので、次回はそのあたりを書いてみようと思います。

DFT計算をiOSでもできるようになりました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

  • このエントリーをはてなブックマークに追加

SNSでもご購読できます。

コメントを残す

*