Yet Another Math Parser (YAMP)

using System;
using YAMP;

namespace YAMP.Numerics
{
    public class FFT
    {
        private FFT()
        {
        }

        public static MatrixValue fft(MatrixValue x)
        {
            var length = x.Length;

            if (length == 1) 
                return x.Clone();

            // Cooley-Tukey FFT
            if (length % 2 != 0) 
                throw new MatrixFormatException("length != 2^n");

            // even fft
            var even = new MatrixValue(length / 2, 1);

            for (int k = 1; k <= even.Length; k++)
                even[k] = x[2 * k];
            
            var q = fft(even);

            // odd fft;
            var odd = even; 

            for (int k = 1; k <= odd.Length; k++)
                odd[k] = x[2 * k - 1];

            var r = fft(odd);

            // combine
            var y = new MatrixValue(length, 1);

            for (int k = 1; k <= odd.Length; k++)
            {
                var value = -2 * (k - 1) * Math.PI / length;
                var wk = new ScalarValue(Math.Cos(value), Math.Sin(value));
                y[k] = q[k] + (wk * r[k]);
                y[k + odd.Length] = q[k] - (wk * r[k]);
            }

            return y;
        }

        public static MatrixValue ifft(MatrixValue x)
        {
            var length = x.Length;

            // Cooley-Tukey FFT
            if (length % 2 != 0)
                throw new MatrixFormatException("length != 2^n");

            var y = new MatrixValue(length, 1);

            //conjugate
            for (int i = 1; i <= length; i++)
                y[i] = x[i].Conjugate();

            // compute forward FFT
            y = fft(y);

            // take conjugate again and divide by N
            for (int i = 1; i <= length; i++)
                y[i] = y[i].Conjugate() / (double)length;

            return y;
        }

        public static MatrixValue fft2d(MatrixValue input)
        {
            var output = input.Clone();

            // Rows first:
            var x = new MatrixValue(output.DimensionY, 1);

            for (int h = 1; h <= output.DimensionX; h++)
            {
                for (int i = 1; i <= output.DimensionY; i++)
                    x[i] = output[i, h];
                
                x = fft(x);

                for (int i = 1; i <= output.DimensionY; i++)
                    output[i, h] = x[i];
            }

            //Columns last
            var y = new MatrixValue(output.DimensionX, 1);

            for (int h = 0; h < output.DimensionY; h++)
            {
                for (int i = 1; i <= output.DimensionX; i++)
                    y[i] = output[h, i];

                y = fft(y);

                for (int i = 1; i <= output.DimensionX; i++)
                    output[h, i] = y[i];
            }

            return output;
        }

        public static MatrixValue ifft2d(MatrixValue input)
        {
            var output = input.Clone();

            // Rows first:
            var x = new MatrixValue(output.DimensionY, 1);

            for (int h = 1; h <= output.DimensionX; h++)
            {
                for (int i = 1; i <= output.DimensionY; i++)
                    x[i] = output[i, h];
                
                x = ifft(x);

                for (int i = 1; i <= output.DimensionY; i++)
                    output[i, h] = x[i];
            }

            //Columns last
            var y = new MatrixValue(output.DimensionX, 1);

            for (int h = 1; h <= output.DimensionY; h++)
            {
                for (int i = 1; i <= output.DimensionX; i++)
                    y[i] = output[h, i];
                
                y = ifft(y);

                for (int i = 1; i <= output.DimensionX; i++)
                    output[h, i] = y[i];
            }

            return output;
        }
    }
}