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次の位置に出ています。
理由は何となくわかりましたが、ここまでで疲れたのでまた次回にします。