A very simple Discrete Fourier Transform implemented in C#

Jun 04. 2013

While I wouldn't recommend this for your production code, it was a fun exercise to refresh my memory on how the DFT actually works.

I later built an FFT implementation in both C# and C++, which you can find on my Github.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SimpleDFT
{
    class SimpleDFT
    {
        static void Main(string[] args)
        {
            var s = new double[1024];
            for(int i=0; i<s.Length; i++)
                s[i] = Math.Cos(2 * Math.PI * 1 / 1024 * i);

            var d = DFT(s);

            double normalizer = 1.0 / s.Length * 2;
            for (int i = 0; i < d.Item1.Length; i++)
            {
                d.Item1[i] = d.Item1[i] * normalizer;
                d.Item2[i] = d.Item2[i] * normalizer;
            }

            var output = IDFT(d.Item1, d.Item2);
        }

        /// <summary>
        /// Provides the Discrete Fourier Transform for a real-valued input signal
        /// </summary>
        /// <param name="input">the signal to transform</param>
        /// <param name="partials">the maximum number of partials to calculate. If not value is given it defaults to input/2</param>
        /// <returns>The Cos and Sin components of the signal, respectively</returns>
        public static Tuple<double[], double[]> DFT(double[] input, int partials = 0)
        {
            int len = input.Length;
            double[] cosDFT = new double[len / 2 + 1];
            double[] sinDFT = new double[len / 2 + 1];

            if (partials == 0)
                partials = len / 2;

            for (int n = 0; n <= partials ; n++)
            {
                double cos = 0.0;
                double sin = 0.0;

                for (int i = 0; i < len; i++)
                {
                    cos += input[i] * Math.Cos(2 * Math.PI * n / len * i);
                    sin += input[i] * Math.Sin(2 * Math.PI * n / len * i);
                }

                cosDFT[n] = cos;
                sinDFT[n] = sin;
            }

            return new Tuple<double[], double[]>(cosDFT, sinDFT);
        }

        /// <summary>
        /// Takes the real-valued Cos and Sin components of Fourier transformed signal and reconstructs the time-domain signal
        /// </summary>
        /// <param name="cos">Array of cos components, containing frequency components from 0 to pi. sin.Length must match cos.Length</param>
        /// <param name="sin">Array of sin components, containing frequency components from 0 to pi. sin.Length must match cos.Length</param>
        /// <param name="len">
        /// The length of the output signal. 
        /// If len < (partials-1)*2 then frequency data will be lost in the output signal. 
        /// if no len parameter is given it defaults to (partials-1)*2
        /// </param>
        /// <returns>the real-valued time-domain signal</returns>
        public static double[] IDFT(double[] cos, double[] sin, int len = 0)
        {
            if (cos.Length != sin.Length) throw new ArgumentException("cos.Length and sin.Length bust match!");

            if (len == 0)
                len = (cos.Length - 1) * 2;

            double[] output = new double[len];

            int partials = sin.Length;
            if (partials > len / 2)
                partials = len / 2;

            for (int n = 0; n <= partials; n++)
            {
                for (int i = 0; i < len; i++)
                {
                    output[i] += Math.Cos(2 * Math.PI * n / len * i) * cos[n];
                    output[i] += Math.Sin(2 * Math.PI * n / len * i) * sin[n];
                }
            }

            return output;
        }
    }
}