Source.

using System;

using System.Collections.Generic;

using System.Text;

 

namespace FFT

{

    public class RealFourier

    {

        private int size;

        private double[] real;

        private double[] img;

 

        public RealFourier(double[] data)

        {

            if (data == null)

                return;

 

            size = data.Length;

            real = new double[data.Length / 2 + 1];

            img = new double[data.Length / 2 + 1];

 

            if (data.Length == 1)

            {

                real[0] = data[0];

                img[0] = 0.0;

            }

            else if (data.Length == 2)

            {

                double scale = 1 / Math.Sqrt(2);

                real[0] = (data[0] + data[1]) * scale;

                real[1] = (data[0] - data[1]) * scale;

                img[0] = 0.0;

                img[1] = 0.0;

            }

            else

            {

                RealFourier even = new RealFourier(data.Even());

                RealFourier odd = new RealFourier(data.Odd());

 

                double PI = Math.Acos(-1);

                double alpha = -2.0 * PI / data.Length;

                double scale = 1 / Math.Sqrt(2);

                for (int k = 0; k <= data.Length / 2; k++)

                {

                    double cos = Math.Cos( alpha * k);

                    double sin = Math.Sin( alpha * k);

                    real[k] = scale* (even.R(k) + cos * odd.R(k) - sin * odd.Y(k));

                    img[k] = scale* (even.Y(k) + sin * odd.R(k) + cos * odd.Y(k));

                }

            }

        }

 

        public int Size { get { return size; } }

 

        public int Count { get { return size / 2 + 1; } }

 

        public double[] Real { get { return real; } }

 

        public double[] Img { get { return img; } }

 

        public double Power(int index)

        {

            double r = R(index);

            double y = Y(index);

            return r * r + y * y;

        }

 

        public double Module(int index)

        {

            return Math.Sqrt(Power(index));

        }

 

        public double R(int index)

        {

            index %= size;

            if (index <= size / 2)

                return real[index];

            return real[size - index];

        }

 

        public double Y(int index)

        {

            index %= size;

            if (index <= size / 2)

                return img[index];

            return -img[size - index];

        }

 

        public double PeakFreq(double freq)

        {

            int pos = 0;

            double max = Module(0);

            for (int index = 1; index <= size / 2; index++)

            {

                double m = Module(index) + Module(2 * index) + Module(3 * index);

                if (max < m)

                {

                    max = m;

                    pos = index;

                }

            }

 

            if (pos == 0 || pos == size / 2)

                return F(pos, freq);

 

            double yMinus = Math.Log(Module(pos - 1));

            double yZero = Math.Log(Module(pos));

            double yPlus = Math.Log(Module(pos + 1));

 

            double delta = (yMinus - yPlus) / (yMinus + yPlus - 2 * yZero);

 

            return freq * (pos + delta) / (double)size;

        }

 

        public double F(int index, double freq)

        {

            index %= size;

            if (index > size / 2)

                index = size - index;

 

            return freq * index / (double)size;

        }

 

        private void Process(double[] data)

        {

            double PI = Math.Acos(-1);

            double alpha = -2.0 * PI / data.Length;

            double sinAlpha = Math.Sin(alpha);

            double cosAlpha = Math.Cos(alpha);

            double scale = 1.0 / Math.Sqrt(data.Length);

 

            real.Set(0.0);

            img.Set(0.0);

 

            double sinKalpha = 0.0;

            double cosKalpha = 1.0;

            for (int k = 0; k <= data.Length / 2; k++)

            {

                double sin = 0.0;

                double cos = 1.0;

                for (int index = 0; index < data.Length; index++)

                {

                    real[k] += data[index] * cos;

                    img[k] += data[index] * sin;

 

                    IncrementAlpha(ref sin, ref cos, sinKalpha, cosKalpha);

                }

 

                IncrementAlpha(ref sinKalpha, ref cosKalpha, sinAlpha, cosAlpha);

            }

 

            real.Scale(scale);

            img.Scale(scale);

        }

 

        private void IncrementAlpha(ref double sin, ref double cos, double sinAlpha, double cosAlpha)

        {

            double sinTmp = sin;

            double cosTmp = cos;

            sin = sinTmp * cosAlpha + cosTmp * sinAlpha;

            cos = cosTmp * cosAlpha - sinTmp * sinAlpha;

        }

    }

}