FFT/Fourier Transforms QuickStart Sample (C#)

Illustrates how to compute the forward and inverse Fourier transform of a real or complex signal using classes in the Extreme.Mathematics.SignalProcessing namespace in C#.

Visual Basic code F# code IronPython code Back to QuickStart Samples

using System;

namespace Extreme.Mathematics.QuickStart.CSharp
{
    // We'll need real vectors...
    using Extreme.Mathematics.LinearAlgebra;
    // and complex vectors...
    using Extreme.Mathematics.LinearAlgebra.Complex;
    // The FFT classes reside in the Extreme.Mathematics.SignalProcessing
    // namespace.
    using Extreme.Mathematics.SignalProcessing;

    /// <summary>
    /// Illustrates the use of the FftProvider and Fft classes for computing 
    /// the Fourier transform of real and complex signals.
    /// </summary>
    class FourierTransforms
    {
        /// <summary>
        /// The main entry point for the application.
        /// </summary>
        [STAThread]
        static void Main(string[] args)
        {
            // This QuickStart sample shows how to compute the Fouier
            // transform of real and complex signals.

            // Some vectors to play with:
            DenseVector r1 = Vector.Create(1000);
            for (int i = 0; i < r1.Length; i++)
                r1[i] = 1.0 / (1 + i);
            ComplexDenseVector c1 = ComplexVector.Create(1000);
            for (int i = 0; i < c1.Length; i++)
                c1[i] = new DoubleComplex(Math.Sin(0.03 * i), Math.Cos(0.07 * i));

            DenseVector r2 = Vector.Create(1.0, 2.0, 3.0, 4.0);
            ComplexDenseVector c2 = ComplexVector.Create(new DoubleComplex(1, 2),
                new DoubleComplex(3, 4), new DoubleComplex(5, 6), new DoubleComplex(7, 8));

            //
            // One-time FFT's
            //

            // The Vector and ComplexVector classes have static methods to compute FFT's:
            ComplexConjugateSignalVector c3 = Vector.FourierTransform(r2);
            Vector r3 = Vector.InverseFourierTransform(c3);
            Console.WriteLine("fft(r2) = {0:F3}", c3);
            Console.WriteLine("ifft(fft(r2)) = {0:F3}", r3);
            // The ComplexConjugateSignalVector type represents a complex vector
            // that is the Fourier transform of a real signal. 
            // It enforces certain symmetry properties:
            Console.WriteLine("c3[i] == conj(c3[N-i]): {0} == conj({1})", c3[1], c3[3]);

            //
            // FFT Providers
            //

            // FFT's require a fair bit of pre-computation. Using the FftProvider class,
            // you can get an Fft object that caches these computations.

            // Here, we create an FFT implementation for a real signal:
            Fft realFft = FftProvider.ManagedProvider.Create1DRealFft(r1.Length);
            // For a complex to complex transform:
            Fft complexFft = FftProvider.ManagedProvider.Create1DComplexFft(c1.Length);

            // You can set the scale factor for the forward transform.
            // The default is 1/N.
            realFft.ForwardScaleFactor = 1.0 / Math.Sqrt(c1.Length);
            // and the backward transform, with default 1:
            realFft.BackwardScaleFactor = realFft.ForwardScaleFactor;

            // The ForwardTransform method performs a forward transform:
            ComplexVector c4 = realFft.ForwardTransform(r1);
            Console.WriteLine("First 5 terms of fft(r1):");
            for (int i = 0; i < 5; i++)
                Console.WriteLine("   {0}: {1}", i, c4[i]);
            c4 = complexFft.ForwardTransform(c1);
            Console.WriteLine("First 5 terms of fft(c1):");
            for (int i = 0; i < 5; i++)
                Console.WriteLine("   {0}: {1}", i, c4[i]);

            // ForwardTransform has many overloads for real to complex and
            // complex to complex transforms.

            // A one-sided transform returns only the first half of the FFT of
            // a real signal. The rest can be deduced from the symmetry properties.
            // Here's how to compute a one-sided FFT:
            ComplexVector c5 = ComplexVector.Create(r1.Length / 2 + 1);
            realFft.ForwardTransform(r1, c5, RealFftFormat.OneSided);

            // The BackwardTransform method has a similar set of overloads:
            Vector r4 = Vector.Create(r1.Length);
            realFft.BackwardTransform(c5, r4, RealFftFormat.OneSided);

            // As the last step, we need to dispose the two FFT implementations.
            realFft.Dispose();
            complexFft.Dispose();

            //
            // 2D transforms
            //

            // 2D transforms are handled in a completely analogous way.
            Matrix m = Matrix.Create(36, 56, delegate (int i, int j)
              { return Math.Exp(-0.1 * i) * Math.Sin(0.01 * (i * i + j * j - i * j)); });
            ComplexMatrix mFft = new ComplexDenseMatrix(m.RowCount, m.ColumnCount);

            using (Fft2D fft2 = FftProvider.CurrentProvider.Create2DRealFft(m.RowCount, m.ColumnCount))
            {
                fft2.ForwardTransform(m, mFft);

                Console.WriteLine("First few terms of fft(m):");
                for (int i = 0; i < 4; i++)
                {
                    string comma = string.Empty;
                    for (int j = 0; j < 4; j++)
                    {
                        Console.Write(comma);
                        Console.Write("{0}", mFft[i, j].ToString("F4"));
                        comma = ", ";
                    }
                    Console.WriteLine();
                }

                // and the backward transform:
                fft2.BackwardTransform(mFft, m);

                // Dispose is called automatically.
            }

            Console.Write("Press Enter key to exit...");
            Console.ReadLine();
        }
    }
}