Zenn
✏️

【復習】Unity で DFT(離散フーリエ変換)を実装してみる

2021/01/16に公開
Unity
tech

はじめに

フーリエ変換の復習に、DFT(離散フーリエ変換)をUnityで実装してみました。

フーリエ変換とは

連続な関数 f(x) を 周波数の関数 F(\omega) に変換することをフーリエ変換と呼びます。

フーリエ変換は以下のような数式で表現されます。

F(\omega) =\int_{-\infty}^{\infty}f(x) \cdot exp(-2 \pi i x \omega) dx

離散フーリエ変換は積分を利用したものになっていますが、
Unity(コンピューター)では積分を直接計算することができません。

Unityで実装する際は離散データを利用する 離散フーリエ変換 を使います。

離散フーリエ変換(DFT)

関数 f(x) を N箇所 x = x_0, x_1, x_2, ... , x_{N-1} でサンプリングすると、
N個のデータ f(x_0), f(x_1), f(x_2), ... , f(x_{N-1}) を得ます。

DFT(離散フーリエ変換)を利用すると、これらのデータ点から周波数スペクトルF(\omega)を求めることができます。

DFTは以下のような式になります。

F(\omega)=\sum_{i=0}^{N-1}f(x_i) \cdot exp(-i \frac{2\pi \omega x_i}{N})

オイラーの公式 exp(ix) = cos x + i sin x より、DFTは以下で計算することができます。

F(\omega)=\sum_{i=0}^{N-1}f(x_i) \cdot \{ cos(\frac{2\pi \omega x_i}{N})- i sin(\frac{2\pi \omega x_i}{N}) \}

UnityでDFTを実装してみる

今回は、Sin波から作った信号データにDFTを適用して周波数スペクトル強度を求めるC#スクリプトを書いてみました。

using UnityEngine;

[RequireComponent(typeof(LineRenderer))]
public class FourierTransform : MonoBehaviour
{
    [SerializeField] private int dataNum = 1024; // 信号のデータ数
    [SerializeField] private int spectrumNum = 32; // 周波数スペクトル数
    [SerializeField] private float sineFreq = 4f; // Sinの周波数(Hz)
    private LineRenderer lineRenderer;

    /// <summary>
    /// スペクトル強度の計算
    /// </summary>
    private void ComputeSpectrum()
    {
        dataNum = Mathf.Max(2, dataNum);
        spectrumNum = Mathf.Max(2, spectrumNum);
        
        // 信号データ 初期化
        var x = new float[dataNum]; // サンプリング点 x
        var data = new float[dataNum]; // 信号のサンプリングして得られるデータ f(x)
        for (int i = 0; i < dataNum; i++)
        {
            x[i] = (float) i / dataNum; // 区間[0, 1] を N分割したものを、サンプリング点とする
            data[i] = Mathf.Sin(x[i] * Mathf.PI * 2f * sineFreq); // Sine波をデータ点とする
        }

        // 周波数スペクトルを計算
        var spectrum = new float[spectrumNum];
        for (int i = 0; i < spectrumNum; i++)
        {
            float t = i * 2f * Mathf.PI / spectrumNum; // 区間 [0, 2π] を分割する点でのスペクトル強度を計算
            spectrum[i] = DFT(t, x, data);
        }

        // LineRendererに周波数スペクトルを反映
        lineRenderer.positionCount = spectrumNum;
        for (int i = 0; i < spectrumNum; i++)
        {
            float t = (float) i / (spectrumNum - 1);
            
            // タテ軸(y)を周波数スペクトルとしてプロット
            float y = spectrum[i] / dataNum; // 見やすいようにdataNumで割っておく
            lineRenderer.SetPosition(i, new Vector3(t, y, 0f));　
        }
    }

    /// <summary>
    /// 離散フーリエ変換を利用して周波数スペクトルを求める
    /// </summary>
    /// <param name="t">周波数スペクトルのサンプリング位置</param>
    /// <param name="x">信号データの全サンプリング点</param>
    /// <param name="data">全信号データ</param>
    /// <returns></returns>
    float DFT(float t, float[] x, float[] data)
    {
        float re = 0f;
        float im = 0f;
        for (int i = 0; i < dataNum; i++)
        {
            re += data[i] * Mathf.Cos(-i * 2f * Mathf.PI * t * x[i] / dataNum);
            im += data[i] * Mathf.Sin(-i * 2f * Mathf.PI * t * x[i] / dataNum);
        }
        // 周波数スペクトル強度を計算
        return Mathf.Sqrt(re * re + im * im);
    }
    
    private void Start()
    {
        lineRenderer = GetComponent<LineRenderer>();
        ComputeSpectrum();
    }
    
    /// <summary>
    /// インスペクターの値が変更された時に呼ばれる
    /// </summary>
    private void OnValidate()
    {
        ComputeSpectrum();
    }

}

結果

周波数1Hzの場合

周波数2Hzの場合

周波数4Hzの場合

Discussion