FFTプログラムの作成（４）

FFTのピーク位置について調べてみました。
FFT実施時にサンプリング個数が2のべき乗に足りない場合は元データを0詰めして2のべき乗にする、というテクニックがありますのでデータ数65536の50001以降のデータを強制的に0にして再度FFT実施しました。
自作、MathNet.Nemricsともにピーク位置は16次に変化ありません。

 続いて強制0ではなく、すべてデータで埋めた50000のサンプリング数でFFTします。2のべき乗ではないので自作プログラムではFFTできませんが、MathNet.NemricsではFFT可能です。
以下がFFT結果です。ピークの位置が12次に移動しました。

どうやらMathNet.Nemrics内部でFFTの手法を切り替えてますね。

調べたところChirp z-変換アルゴリズム（Bluestein FFT）という手法を用いれば2のべき乗でなくともFFTできるようなので、それを使っているのかもしれません。

今回の結果から以下のように理解しました（個人の感想です！）。

情報が不明な信号を解析する場合は0埋めFFTで良いと思いますが、基準の信号波長が既知でその強度を調査したい、等の用途には0埋めは不向きな場合があるようです。
といいますか、65536のデータのうち1/4近くの15536個も0詰めしているのが駄目なんでしょうね。既知の信号を解析する場合は、できるだけ2のべき乗に近しいサンプリング数になるようサンプリングレートを決める必要がある、という当たり前の結論なのか。
いちおう納得したので2のべき乗でないFFTの実装チャレンジしてみます。

FFTプログラムの作成（３）

少し暇ができたので見直し。Complex構造体を使ってみました。ビット反転をもう少し効率化したい。

        public static void FFT2Radix2(ref Complex[] Wave)
        {
            int Number = Wave.Length;
            const double PI = Math.PI;
            int StageCount;
            int BlockCount;
            int NodeCount;
            int StageNumber;
            int BlockNumber;
            int NodeNumber;
            int[] IndexBitChange = new int[Number];

            int TempN1;
            int TempN2;

            Complex Temp1;
            Complex Temp2;
            Complex RotationFactor;
            Complex ButterflyResolution;

            //計算準備
            //累乗数を計算
            StageCount = PrimeFactors(Number).Length;

            //バタフライ演算の繰り返しループ
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                //回転因子の分割数
                BlockCount = (int)(Math.Pow(2, StageNumber));
                //分割した回転因子内の信号数
                NodeCount = Number / BlockCount;
                //回転因子の刻み度数
                ButterflyResolution = new Complex(0, -2 * PI / Number * BlockCount);

                //回転因子の分割数に達するまでループを回す
                for (BlockNumber = 0; BlockNumber < BlockCount; BlockNumber++)
                {
                    //バタフライ演算
                    for (NodeNumber = 0; NodeNumber < NodeCount / 2; NodeNumber++)
                    {
                        //回転因子
                        RotationFactor = Complex.Exp(ButterflyResolution * NodeNumber);
                        //バタフライ演算
                        TempN1 = NodeNumber + NodeCount * BlockNumber;
                        TempN2 = TempN1 + NodeCount / 2;
                        Temp1 = Wave[TempN1];
                        Temp2 = Wave[TempN2];
                        Wave[TempN1] = Temp1 + Temp2;
                        Wave[TempN2] = (Temp1 - Temp2) * RotationFactor;
                    }
                }
            }
            //ビット反転
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                for (int i = 0; i < Math.Pow(2, StageNumber); i++)
                {
                    IndexBitChange[(int)Math.Pow(2, StageNumber) + i] = IndexBitChange[i] +
 　　　　　　　　(int)Math.Pow(2, StageCount - StageNumber - 1);
                }
            }

            for (int i = 0; i < Number; i++)
            {
                if (IndexBitChange[i] > i)
                {
                    Temp1 = Wave[IndexBitChange[i]];
                    Wave[IndexBitChange[i]] = Wave[i];
                    Wave[i] = Temp1;
                }
            }
        }

FFTプログラムの作成（２）

FFTの理解を深めるため、ネットで先達が提示されていたVBのプログラムをC#に移植しました。もっと効率的な書き方できるんだろうけど、理解深めるために愚直に書きます。
実際にN=8=2³の場合で計算順を辿って考えてみるとバタフライ演算の流れが少しずつ理解できました。
最初に大きな回転因子で計算し、それを半分ずつ小さな回転因子に分割していくようなイメージでしょうか？
実際にサンプルデータを作ってみて検算してみました。
        public static void FFT2Radix(double[] InputWave, out double[] OutputWaveReal, out double[] OutputWaveImaginary)
        {
            int Number = InputWave.Length;
            OutputWaveReal = new double[Number];
            OutputWaveImaginary = new double[Number];

            int StageCount;
            int BlockCount;
            int NodeCount;
            int StageNumber;
            int BlockNumber;
            int NodeNumber;
            int[] index = new int[Number];

            int TempN1;
            int TempN2;

            double ButterflyResolution;
            int[] IndexBitChange = new int[Number];
            double Temp1_Real;
            double Temp1_Imaginary;
            double Temp2_Real;
            double Temp2_Imaginary;
            double RotationFactorCos;
            double RotationFactorSin;

            for (int i = 0; i < Number; i++)
            {
                OutputWaveReal[i] = InputWave[i];
                OutputWaveImaginary[i] = 0;
            }

            //計算準備
            //累乗数を計算
            StageCount = PrimeFactors(Number).Length;

            //バタフライ演算の繰り返しループ
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                //回転因子の分割数
                BlockCount = (int)(Math.Pow(2, StageNumber));
                //分割した回転因子内の信号数
                NodeCount = Number / BlockCount;
                //回転因子の刻み度数
                ButterflyResolution = 2 * PI / Number * BlockCount;

                //回転因子の分割数に達するまでループを回す
                for (BlockNumber = 0; BlockNumber < BlockCount; BlockNumber++)
                {
                    //バタフライ演算
                    for (NodeNumber = 0; NodeNumber < NodeCount / 2; NodeNumber++)
                    {
                        RotationFactorCos = Math.Cos(ButterflyResolution * NodeNumber);
                        RotationFactorSin = Math.Sin(ButterflyResolution * NodeNumber);
                        //バタフライ演算
                        TempN1 = NodeNumber + NodeCount * BlockNumber;
                        TempN2 = TempN1 + NodeCount / 2;
                        Temp1_Real = OutputWaveReal[TempN1];
                        Temp1_Imaginary = OutputWaveImaginary[TempN1];
                        Temp2_Real = OutputWaveReal[TempN2];
                        Temp2_Imaginary = OutputWaveImaginary[TempN2];
                        OutputWaveReal[TempN1] = Temp1_Real + Temp2_Real;
                        OutputWaveImaginary[TempN1] = Temp1_Imaginary + Temp2_Imaginary;
                        OutputWaveReal[TempN2] = (Temp1_Real - Temp2_Real) * RotationFactorCos
                            + (Temp1_Imaginary - Temp2_Imaginary) * RotationFactorSin;
                        OutputWaveImaginary[TempN2] = -(Temp1_Real - Temp2_Real) * RotationFactorSin
                            + (Temp1_Imaginary - Temp2_Imaginary) * RotationFactorCos;
                    }
                }
            }
            //インデックス入れ替え
            for (StageNumber = 0; StageNumber < StageCount; StageNumber++)
            {
                for (int i = 0; i < Math.Pow(2, StageNumber); i++)
                {
                    index[(int)Math.Pow(2, StageNumber) + i] = index[i] + (int)Math.Pow(2, StageCount - StageNumber - 1);
                }
            }

            for (int i = 0; i < Number; i++)
            {
                if (index[i] > i)
                {
                    Temp1_Real = OutputWaveReal[index[i]];
                    Temp1_Imaginary = OutputWaveImaginary[index[i]];
                    OutputWaveReal[index[i]] = OutputWaveReal[i];
                    OutputWaveImaginary[index[i]] = OutputWaveImaginary[i];
                    OutputWaveReal[i] = Temp1_Real;
                    OutputWaveImaginary[i] = Temp1_Imaginary;
                }
            }

MathNet.Numericsの計算結果と照らし合わせて検算しましたが全く同じ結果となり、まずはうまくいったようです。元の波形は以下。
Vector[i] = 3 + 0.5 * Math.Cos(2 * Math.PI * 12.0 * 2.0* (i / 100000))
振幅0.5Vの正弦波にDCが3V乗ったイメージです。
直流成分をサンプル数で割るとほぼ3になるので問題なし。ただし本来なら12次の位置にピークが出ると思っていたら16次の位置に出ています。
理由は何となくわかりましたが、ここまでで疲れたのでまた次回にします。