Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(3) 高速フーリエ変換（FFT）を実行する。FourierOptionsにMatlabとNumerical Recipesがあるのが意外。

さて、前回、前々回はMath.NET numericsを使って

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(1) 複素行列を定義して一次方程式や逆行列、行列式などを計算する。

Visual C# (C_sharp)の数学ライブラリ Math.NET Numericsを使う(2) 補間を行う（Interpolate) リニア、３次スプライン、有理関数などいろいろ使える。

をやってみた。今回はフーリエ変換。

https://numerics.mathdotnet.com/api/MathNet.Numerics.IntegralTransforms/Fourier.htm

ちょっと面白いのはオプションで、

FourierOptions Matlab

FourierOptions NumericalRecipes

MatlabはいいとしてNumerical Recipesって本当に変なFFTの定義していて、ただそれを私もExcel VBAに移植して使っていたりするので実はおなじみ。まあ普通の人はMatlab定義を使えばいいと思います。Numpyとも同じだし。

例題もNumpyからとってこよう。

https://numpy.org/doc/stable/reference/generated/numpy.fft.fft.html

結果はこんな感じで、

ソースコードはこんなの。

ソースはテキストでも書いておこう。

    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
            // 初期化
            chart1.Series.Clear();
            chart1.Titles.Clear();
            chart1.Legends.Clear();
            chart1.ChartAreas.Clear();
        }

        private void button1_Click(object sender, EventArgs e)
        {

            Title title = new Title("FFT", Docking.Top);
            chart1.Titles.Add(title);
            Legend legend = new Legend();
            chart1.Legends.Add(legend);

            Series series1 = new Series();
            series1.ChartType = SeriesChartType.Line;
            series1.BorderWidth = 1;
            series1.LegendText = "Real part";
            chart1.Series.Add(series1);

            Series series2 = new Series();
            series2.ChartType = SeriesChartType.Line;
            series2.BorderWidth = 1;
            series2.LegendText = "Imaginary part";
            chart1.Series.Add(series2);
            chart1.ChartAreas.Add("");
            Axis axisX = new Axis();
            axisX.Title = "X軸";
            //axisX.Minimum = 0;
            //axisX.Maximum = 7;
            //axisX.Interval = 0.5;
            chart1.ChartAreas[0].AxisX = axisX;

            Axis axisY = new Axis();
            axisY.Title = "Y軸";
            //axisY.Minimum = -1;
            //axisY.Maximum = 1;
            //axisY.Interval = 0.2;
            chart1.ChartAreas[0].AxisY = axisY;

            double[] t = new double[256];
            Complex[] y = new Complex[256];
            Complex imag = Complex.ImaginaryOne;

            for (int i = 0; i < t.Length; i++)
            {
                t[i] = Convert.ToDouble(i);
                y[i] = Math.Sin(t[i]);
            }

            Fourier.Forward(y, FourierOptions.Matlab);

            double[] freq = Fourier.FrequencyScale(256, 1.0);

            for (int i = 0; i < freq.Length; i++)
            {
                series1.Points.AddXY(freq[i], y[i].Real);
                series2.Points.AddXY(freq[i], y[i].Imaginary);
            }

        }
    }