LoginSignup
0
1

More than 3 years have passed since last update.

@ryujimiya2361(ryujimiya)

Ooura FFTをネイティブDLL(x64)にしてC#から使ってみた

Last updated at Posted at 2019-09-17

1. はじめに

FFTを使う予定がありそうなので以前使ってたILNumericsでもNuGetしてみようかなと思ったら、いつのまにか有料になっていた。
FFTのライブラリを探していると、大浦FFTなるものが元になっているライブラリがいくつかあることがわかった。
　　大浦FFTのパッケージページ：
　　http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html
これらはライブラリ化はしておらず、プロジェクトにファイルを追加することを想定しているようです。

2. スタティックライブラリの作成

大浦FFTのソースコードをVisual Studio 2019でスタティックライブラリにしました。
大浦FFTには次のの3ファイルあります。
　fft4g.c fft8g.c fftsg.c

これらは同じ関数名が使われているのでそのままで共用できませんでした。外部から参照する関数名は語尾に4gなどをつけることで解消しました。
また外部参照しない関数にはstaticを付けてファイル内スコープになるようにしました。
これでスタティックライブラリの作成に成功しました。

OouraFFTLib.h
#pragma once

void cdft4g(int n, int isgn, double* a, int* ip, double* w);
void rdft4g(int n, int isgn, double* a, int* ip, double* w);
void ddct4g(int n, int isgn, double* a, int* ip, double* w);
void ddst4g(int n, int isgn, double* a, int* ip, double* w);
void dfct4g(int n, double* a, double* t, int* ip, double* w);
void dfst4g(int n, double* a, double* t, int* ip, double* w);

void cdft8g(int n, int isgn, double* a, int* ip, double* w);
void rdft8g(int n, int isgn, double* a, int* ip, double* w);
void ddct8g(int n, int isgn, double* a, int* ip, double* w);
void ddst8g(int n, int isgn, double* a, int* ip, double* w);
void dfct8g(int n, double* a, double* t, int* ip, double* w);
void dfst8g(int n, double* a, double* t, int* ip, double* w);

void cdftsg(int n, int isgn, double* a, int* ip, double* w);
void rdftsg(int n, int isgn, double* a, int* ip, double* w);
void ddctsg(int n, int isgn, double* a, int* ip, double* w);
void ddstsg(int n, int isgn, double* a, int* ip, double* w);
void dfctsg(int n, double* a, double* t, int* ip, double* w);
void dfstsg(int n, double* a, double* t, int* ip, double* w);

3. DLLの作成

2.で作成したlibをリンクするようにし、必要な外部参照関数を定義しました。

OouraFFTDLL.h
#pragma once

#ifdef OOURAFFTDll_EXPORTS
#define DLL_API __declspec(dllexport)
#else
#define DLL_API __declspec(dllimport)
#endif


extern "C"
{
    DLL_API void OouraFFT_cdft4g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_rdft4g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_ddct4g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_ddst4g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_dfct4g(int n, double* a, double* t, int* ip, double* w);
    DLL_API void OouraFFT_dfst4g(int n, double* a, double* t, int* ip, double* w);

    DLL_API void OouraFFT_cdft8g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_rdft8g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_ddct8g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_ddst8g(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_dfct8g(int n, double* a, double* t, int* ip, double* w);
    DLL_API void OouraFFT_dfst8g(int n, double* a, double* t, int* ip, double* w);

    DLL_API void OouraFFT_cdftsg(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_rdftsg(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_ddctsg(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_ddstsg(int n, int isgn, double* a, int* ip, double* w);
    DLL_API void OouraFFT_dfctsg(int n, double* a, double* t, int* ip, double* w);
    DLL_API void OouraFFT_dfstsg(int n, double* a, double* t, int* ip, double* w);
}
OouraFFTDll.cpp
#include "pch.h"

#include "OouraFFTDll.h"
#include "OouraFFTLib.h"

void OouraFFT_cdft4g(int n, int isgn, double* a, int* ip, double* w)
{
    cdft4g(n, isgn, a, ip, w);

}

void OouraFFT_rdft4g(int n, int isgn, double* a, int* ip, double* w)
{
    rdft4g(n, isgn, a, ip, w);
}
(省略)

4. C#でOouraFFTDll.dllを使うときのメソッドの定義

FFT.ImportedFunctions.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Runtime.InteropServices; // DllImport

namespace IvyFEM.FFT
{
    public class ImportedFunctions
    {
        ///////////////////////////////////////////////////////////////////
        // OouraFFT
        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_cdft4g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_rdft4g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_ddct4g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_ddst4g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_dfct4g(int n, double* a, double* t, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_dfst4g(int n, double* a, double* t, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_cdft8g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_rdft8g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_ddct8g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_ddst8g(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_dfct8g(int n, double* a, double* t, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_dfst8g(int n, double* a, double* t, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_cdftsg(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_rdftsg(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_ddctsg(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_ddstsg(int n, int isgn, double* a, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_dfctsg(int n, double* a, double* t, int* ip, double* w);

        [DllImport("OouraFFTDll.dll")]
        public static extern unsafe void OouraFFT_dfstsg(int n, double* a, double* t, int* ip, double* w);
    }
}

5. FFT関数を使う

        public static void DoFFT(
            double[] times,
            double[] timeDomainDatas,
            out double[] freqs,
            out System.Numerics.Complex[] freqDomaindatas)
        {
            freqs = null;
            freqDomaindatas = null;

            System.Diagnostics.Debug.Assert(times.Length == timeDomainDatas.Length);
            int dataCnt = timeDomainDatas.Length;
            double dt = times[1] - times[0];

            // FFTを実行する
            {
                int n = dataCnt;
                int isgn = 1; // 1 or -1;
                double[] a = new double[n * 2];
                /*
                a[0...2 * n - 1]   :input / output data(double *)
                        input data
                            a[2 * j] = Re(x[j]), 
                            a[2 * j + 1] = Im(x[j]), 0 <= j < n
                        output data
                            a[2 * k] = Re(X[k]), 
                            a[2 * k + 1] = Im(X[k]), 0 <= k < n

                X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
                */
                for (int i = 0; i < n; i++)
                {
                    double data = timeDomainDatas[i];
                    a[2 * i] = data;
                    a[2 * i + 1] = 0;
                }

                int ipSize = (int)(2 + Math.Sqrt(n)) + 1;
                int[] ip = new int[ipSize];
                ip[0] = 0;
                // w[],ip[] are initialized if ip[0] == 0

                int wSize = n / 2;
                double[] w = new double[wSize];

                unsafe
                {
                    fixed (double* aP = &a[0])
                    fixed (int* ipP = &ip[0])
                    fixed (double* wP = &w[0])
                    {
                        IvyFEM.FFT.ImportedFunctions.OouraFFT_cdftsg(2 * n, isgn, aP, ipP, wP);
                    }
                }

                freqDomaindatas = new System.Numerics.Complex[n];
                for (int i = 0; i < n; i++)
                {
                    double real = a[2 * i];
                    double imag = a[2 * i + 1];
                    // フーリエ複素振幅
                    System.Numerics.Complex data = new System.Numerics.Complex(real, imag);
                    freqDomaindatas[i] = data;
                }
            }

            // 時間長さ
            double tl = dt * dataCnt;
            // 周波数刻み
            double df = 1.0 / tl;
            // 周波数アレイ
            freqs = new double[dataCnt];
            for (int i = 0; i < dataCnt; i++)
            {
                freqs[i] = i * df;
            }
        }

6. 計算例

【注意】サンプル数nは2^mにしないと駄目なようです。

        public void FFTExample(MainWindow mainWindow)
        {
            int n = 8;
            double[] times = new double[n];
            double[] datas = new double[n];
            for (int i = 0; i < n; i++)
            {
                times[i] = i;
                datas[i] = i < 4 ? 1.0 : 0.0;
            }

            string CRLF = System.Environment.NewLine;
            string ret =
                "FFT Example" + CRLF +
                "---------------------------------" + CRLF +
                "datas" + CRLF;
            for (int i = 0; i < n; i++)
            {
                ret += string.Format("{0}", datas[i]) + CRLF;
            }
            ret += "---------------------------------" + CRLF;

            double[] freqs;
            System.Numerics.Complex[] freqDomainDatas;
            IvyFEM.FFT.Functions.DoFFT(times, datas, out freqs, out freqDomainDatas);
            ret +=
                "Result" + CRLF +
                "---------------------------------" + CRLF +
                "frequency domain data" + CRLF;
            for (int i = 0; i < n; i++)
            {
                ret += string.Format("{0}", freqDomainDatas[i]) + CRLF;
            }
            ret += "---------------------------------" + CRLF;

            System.Diagnostics.Debug.WriteLine(ret);
        }

実行結果
20190917_FFT_example.jpg

7. まとめ

Ooura FFTをネイティブDLL(x64)にしてC#から使えるようにしました。
このOouraFFTDll.dllは下記からダウンロードできます。
https://github.com/ryujimiya/OouraFFTDll

[追記]
このdllを作成した目的は、下記の計算を行いたかったからです。
IvyFEM開発日誌（10）.NET 有限要素法 IvyFEM 0.0.0.21 - 時間領域FEMによるH面導波管直角コーナーベンドの散乱係数の計算(Givoli-Neta-Patlashenko's High order ABC)

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
1