diff --git a/Sources/Accord.Math/AForge.Math/FourierTransform.cs b/Sources/Accord.Math/AForge.Math/FourierTransform.cs index 8b226dde2..9e6236fe7 100644 --- a/Sources/Accord.Math/AForge.Math/FourierTransform.cs +++ b/Sources/Accord.Math/AForge.Math/FourierTransform.cs @@ -13,14 +13,27 @@ namespace Accord.Math { using System; using Accord.Compat; - using System.Numerics; - + using System.Numerics; + using Accord.Math.Transforms; + /// - /// Fourier transformation. + /// Original Fourier transform from AForge.NET. If possible, + /// please use instead. /// /// - /// The class implements one dimensional and two dimensional - /// Discrete and Fast Fourier Transformation. + /// + /// + /// The class implements one dimensional and two dimensional Discrete and Fast Fourier + /// Transformation. However, this class works only with square matrices with sizes that + /// are power of 2, and implements a different form of the transform that differs from + /// the implementations in other packages such as Octave and Matlab. For a more general + /// transform that should produce the same results as Octave, see . + /// + /// + /// This class may be deprecated (marked as obsolete) in the future. + /// + /// + /// /// public static class FourierTransform { diff --git a/Sources/Accord.Math/Matrix/Matrix.Complex.cs b/Sources/Accord.Math/Matrix/Matrix.Complex.cs index c68dd235e..178abd33f 100644 --- a/Sources/Accord.Math/Matrix/Matrix.Complex.cs +++ b/Sources/Accord.Math/Matrix/Matrix.Complex.cs @@ -57,7 +57,8 @@ public static Complex[] Abs(this Complex[] x) /// public static Complex Sum(this Complex[] x) { - if (x == null) throw new ArgumentNullException("x"); + if (x == null) + throw new ArgumentNullException("x"); Complex r = Complex.Zero; for (int i = 0; i < x.Length; i++) @@ -71,8 +72,10 @@ public static Complex Sum(this Complex[] x) /// public static Complex[] Multiply(this Complex[] a, Complex[] b) { - if (a == null) throw new ArgumentNullException("a"); - if (b == null) throw new ArgumentNullException("b"); + if (a == null) + throw new ArgumentNullException("a"); + if (b == null) + throw new ArgumentNullException("b"); Complex[] r = new Complex[a.Length]; for (int i = 0; i < a.Length; i++) @@ -88,14 +91,7 @@ public static Complex[] Multiply(this Complex[] a, Complex[] b) /// public static double[] Magnitude(this Complex[] c) { - if (c == null) - throw new ArgumentNullException("c"); - - double[] magnitudes = new double[c.Length]; - for (int i = 0; i < c.Length; i++) - magnitudes[i] = c[i].Magnitude; - - return magnitudes; + return c.Apply((x, i) => x.Magnitude); } /// @@ -104,17 +100,16 @@ public static double[] Magnitude(this Complex[] c) /// public static double[,] Magnitude(this Complex[,] c) { - if (c == null) throw new ArgumentNullException("c"); - - int rows = c.GetLength(0); - int cols = c.GetLength(1); - - double[,] magnitudes = new double[rows, cols]; - for (int i = 0; i < rows; i++) - for (int j = 0; j < cols; j++) - magnitudes[i, j] = c[i, j].Magnitude; + return c.Apply((x, i, j) => x.Magnitude); + } - return magnitudes; + /// + /// Gets the magnitude of every complex number in a matrix. + /// + /// + public static double[][] Magnitude(this Complex[][] c) + { + return c.Apply((x, i, j) => x.Magnitude); } /// @@ -123,13 +118,25 @@ public static double[] Magnitude(this Complex[] c) /// public static double[] Phase(this Complex[] c) { - if (c == null) throw new ArgumentNullException("c"); + return c.Apply((x, i) => x.Phase); + } - double[] phases = new double[c.Length]; - for (int i = 0; i < c.Length; i++) - phases[i] = c[i].Phase; + /// + /// Gets the phase of every complex number in a matrix. + /// + /// + public static double[,] Phase(this Complex[,] c) + { + return c.Apply((x, i, j) => x.Phase); + } - return phases; + /// + /// Gets the phase of every complex number in a matrix. + /// + /// + public static double[][] Phase(this Complex[][] c) + { + return c.Apply((x, i, j) => x.Phase); } /// @@ -142,13 +149,7 @@ public static double[] Phase(this Complex[] c) /// public static double[] Re(this Complex[] c) { - if (c == null) throw new ArgumentNullException("c"); - - double[] re = new double[c.Length]; - for (int i = 0; i < c.Length; i++) - re[i] = c[i].Real; - - return re; + return c.Apply((x, i) => x.Real); } /// @@ -161,18 +162,20 @@ public static double[] Re(this Complex[] c) /// public static double[,] Re(this Complex[,] c) { - if (c == null) - throw new ArgumentNullException("c"); - - int rows = c.GetLength(0); - int cols = c.GetLength(1); - - var re = new double[rows, cols]; - for (int i = 0; i < rows; i++) - for (int j = 0; j < cols; j++) - re[i, j] = c[i, j].Real; + return c.Apply((x, i, j) => x.Real); + } - return re; + /// + /// Returns the real matrix part of the complex matrix c. + /// + /// + /// A matrix of complex numbers. + /// + /// A matrix of scalars with the real part of the complex numbers. + /// + public static double[][] Re(this Complex[][] c) + { + return c.Apply((x, i, j) => x.Real); } /// @@ -186,14 +189,7 @@ public static double[] Re(this Complex[] c) // TODO: Rename to Imaginary public static double[] Im(this Complex[] c) { - if (c == null) - throw new ArgumentNullException("c"); - - double[] im = new double[c.Length]; - for (int i = 0; i < c.Length; i++) - im[i] = c[i].Imaginary; - - return im; + return c.Apply((x, i) => x.Imaginary); } /// @@ -203,18 +199,17 @@ public static double[] Im(this Complex[] c) /// A matrix of scalars with the imaginary part of the complex numbers. public static double[,] Im(this Complex[,] c) { - if (c == null) - throw new ArgumentNullException("c"); - - int rows = c.GetLength(0); - int cols = c.GetLength(1); - - var im = new double[rows, cols]; - for (int i = 0; i < rows; i++) - for (int j = 0; j < cols; j++) - im[i, j] = c[i, j].Imaginary; + return c.Apply((x, i, j) => x.Imaginary); + } - return im; + /// + /// Returns the imaginary matrix part of the complex matrix c. + /// + /// A matrix of complex numbers. + /// A matrix of scalars with the imaginary part of the complex numbers. + public static double[][] Im(this Complex[][] c) + { + return c.Apply((x, i, j) => x.Imaginary); } /// diff --git a/Sources/Accord.Math/Transforms/FourierTransform2.cs b/Sources/Accord.Math/Transforms/FourierTransform2.cs index 98454a45c..2df738466 100644 --- a/Sources/Accord.Math/Transforms/FourierTransform2.cs +++ b/Sources/Accord.Math/Transforms/FourierTransform2.cs @@ -59,12 +59,17 @@ namespace Accord.Math.Transforms /// /// /// - /// This fourier transform accepts arbitrary-length matrices and is not - /// restricted only to matrices that have dimensions which are powers of - /// two. It also provides results which are more equivalent with other - /// mathematical packages, such as MATLAB and Octave. + /// + /// The transforms in this class accept arbitrary-length matrices and are not restricted to + /// only matrices that have dimensions which are powers of two. It also provides results which + /// are more equivalent with other mathematical packages, such as MATLAB and Octave. + /// + /// This class had been created as an alternative to AForge.NET's + /// original FourierTransform class that would provide more expected results. /// /// + /// + /// public static class FourierTransform2 { @@ -72,9 +77,11 @@ public static class FourierTransform2 /// 1-D Discrete Fourier Transform. /// /// - /// The data to transform.. + /// The data to transform. /// The transformation direction. /// + /// + /// public static void DFT(Complex[] data, FourierTransform.Direction direction) { int n = data.Length; @@ -124,41 +131,61 @@ public static void DFT(Complex[] data, FourierTransform.Direction direction) /// The data to transform. /// The transformation direction. /// + /// + /// public static void DFT2(Complex[][] data, FourierTransform.Direction direction) { - int n = data.Rows(); int m = data.Columns(); - // process rows - var row = new Complex[m]; - for (int i = 0; i < n; i++) + if (direction == FourierTransform.Direction.Forward) { - // copy row - for (int j = 0; j < row.Length; j++) - row[j] = data[i][j]; + // process rows + for (int i = 0; i < data.Length; i++) + { + // transform it + DFT(data[i], FourierTransform.Direction.Forward); + } - // transform it - DFT(row, direction); + // process columns + var col = new Complex[data.Length]; + for (int j = 0; j < m; j++) + { + // copy column + for (int i = 0; i < col.Length; i++) + col[i] = data[i][j]; - // copy back - for (int j = 0; j < row.Length; j++) - data[i][j] = row[j]; - } + // transform it + DFT(col, FourierTransform.Direction.Forward); - // process columns - var col = new Complex[n]; - for (int j = 0; j < m; j++) + // copy back + for (int i = 0; i < col.Length; i++) + data[i][j] = col[i]; + } + } + else { - // copy column - for (int i = 0; i < col.Length; i++) - col[i] = data[i][j]; + // process columns + var col = new Complex[data.Length]; + for (int j = 0; j < m; j++) + { + // copy column + for (int i = 0; i < col.Length; i++) + col[i] = data[i][j]; - // transform it - DFT(col, direction); + // transform it + DFT(col, FourierTransform.Direction.Backward); - // copy back - for (int i = 0; i < col.Length; i++) - data[i][j] = col[i]; + // copy back + for (int i = 0; i < col.Length; i++) + data[i][j] = col[i]; + } + + // process rows + for (int i = 0; i < data.Length; i++) + { + // transform it + DFT(data[i], FourierTransform.Direction.Backward); + } } } @@ -169,6 +196,8 @@ public static void DFT2(Complex[][] data, FourierTransform.Direction direction) /// The data to transform.. /// The transformation direction. /// + /// + /// public static void FFT(Complex[] data, FourierTransform.Direction direction) { int n = data.Length; @@ -212,6 +241,8 @@ public static void FFT(Complex[] data, FourierTransform.Direction direction) /// The imaginary part of the complex numbers to transform. /// The transformation direction. /// + /// + /// public static void FFT(double[] real, double[] imag, FourierTransform.Direction direction) { if (direction == FourierTransform.Direction.Forward) @@ -237,9 +268,11 @@ public static void FFT(double[] real, double[] imag, FourierTransform.Direction /// 2-D Fast Fourier Transform. /// /// - /// The data to transform.. + /// The data to transform. /// The Transformation direction. /// + /// + /// public static void FFT2(Complex[][] data, FourierTransform.Direction direction) { int n = data.Length; diff --git a/Unit Tests/Accord.Tests.Math/FourierTransformTest.cs b/Unit Tests/Accord.Tests.Math/FourierTransformTest.cs index 2e7746228..1dd4e4635 100644 --- a/Unit Tests/Accord.Tests.Math/FourierTransformTest.cs +++ b/Unit Tests/Accord.Tests.Math/FourierTransformTest.cs @@ -34,6 +34,163 @@ namespace Accord.Tests.Math public class FourierTransformTest { + [Test] + public void dft1d_test() + { + #region doc_dft + // Example from http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf + Complex[] input = new Complex[] { 8, 4, 8, 0 }; + FourierTransform2.DFT(input, FourierTransform.Direction.Forward); + Complex[] output = input; // the output will be { 20, -4j, 12, 4j } + + // We can also get the real and imaginary parts of the vector using + double[] re = output.Re(); // should be { 20, 0, 12, 0 } + double[] im = output.Im(); // should be { 0, -4, 0, 4 } + + // We can recover the initial vector using the backward transform: + FourierTransform2.DFT(output, FourierTransform.Direction.Backward); + #endregion + + double[] re2 = output.Re(); + double[] im2 = output.Im(); + + Assert.IsTrue(re.IsEqual(new[] { 20, 0, 12, 0 }, 1e-8)); + Assert.IsTrue(im.IsEqual(new[] { 0, -4, 0, 4 }, 1e-8)); + Assert.IsTrue(re2.IsEqual(new[] { 8, 4, 8, 0 }, 1e-8)); + Assert.IsTrue(im2.IsEqual(new[] { 0, 0, 0, 0 }, 1e-8)); + } + + [Test] + public void dft2d_test() + { + #region doc_dft2 + // Suppose we would like to transform the matrix + Complex[][] input = new Complex[][] + { + new Complex[] { 1, 4, 8, 0 }, + new Complex[] { 0, 4, 5, 0 }, + new Complex[] { 0, 0, 8, 0 }, + }; + + // We can forward it to the time domain using the Fourier transform + FourierTransform2.DFT2(input, FourierTransform.Direction.Forward); + + // We can also get the real and imaginary parts of the matrix using + double[][] re = input.Re(); // The real part of the output will be: + double[][] expectedRe = new double[][] + { + new double[] { 30.0, -20.0, 14.0, -20.0 }, + new double[] { 4.5, -3.9641, 0.4999, 2.9641 }, + new double[] { 4.5, 2.9641, 0.5000, -3.9641 } + }; + + double[][] im = input.Im(); // The imaginary part of the output will be: + double[][] expectedIm = new double[][] + { + new double[] { 0, -8, 0, 7.9999 }, + new double[] { -0.8660, -4.5980, 6.0621, -0.5980 }, + new double[] { 0.8660, 0.5980, -6.0621, 4.5980 } + }; + + // We can recover the initial matrix using the backward transform: + FourierTransform2.DFT2(input, FourierTransform.Direction.Backward); + #endregion + + string a = re.ToCSharp(); + string b = im.ToCSharp(); + + Assert.IsTrue(re.IsEqual(expectedRe, 1e-3)); + Assert.IsTrue(im.IsEqual(expectedIm, 1e-3)); + var expected = new Complex[][] + { + new Complex[] { 1, 4, 8, 0 }, + new Complex[] { 0, 4, 5, 0 }, + new Complex[] { 0, 0, 8, 0 }, + }; + + Assert.IsTrue(input.Re().IsEqual(expected.Re(), 1e-3)); + Assert.IsTrue(input.Im().IsEqual(expected.Im(), 1e-3)); + } + + [Test] + public void fft1d_test() + { + #region doc_fft + // Example from http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf + Complex[] input = new Complex[] { 8, 4, 8, 0 }; + FourierTransform2.FFT(input, FourierTransform.Direction.Forward); + Complex[] output = input; // the output will be { 20, -4j, 12, 4j } + + // We can also get the real and imaginary parts of the vector using + double[] re = output.Re(); // should be { 20, 0, 12, 0 } + double[] im = output.Im(); // should be { 0, -4, 0, 4 } + + // We can recover the initial vector using the backward transform: + FourierTransform2.FFT(output, FourierTransform.Direction.Backward); + #endregion + + double[] re2 = output.Re(); + double[] im2 = output.Im(); + + Assert.IsTrue(re.IsEqual(new[] { 20, 0, 12, 0 }, 1e-8)); + Assert.IsTrue(im.IsEqual(new[] { 0, -4, 0, 4 }, 1e-8)); + Assert.IsTrue(re2.IsEqual(new[] { 8, 4, 8, 0 }, 1e-8)); + Assert.IsTrue(im2.IsEqual(new[] { 0, 0, 0, 0 }, 1e-8)); + } + + [Test] + public void fft2d_test() + { + #region doc_fft2 + // Suppose we would like to transform the matrix + Complex[][] input = new Complex[][] + { + new Complex[] { 1, 4, 8, 0 }, + new Complex[] { 0, 4, 5, 0 }, + new Complex[] { 0, 0, 8, 0 }, + }; + + // We can forward it to the time domain using the Fourier transform + FourierTransform2.FFT2(input, FourierTransform.Direction.Forward); + + // We can also get the real and imaginary parts of the matrix using + double[][] re = input.Re(); // The real part of the output will be: + double[][] expectedRe = new double[][] + { + new double[] { 30, -20, 14, -20 }, + new double[] { 4.5, -3.9641, 0.4999, 2.9641 }, + new double[] { 4.5, 2.9641, 0.5000, -3.9641 }, + }; + + double[][] im = input.Im(); // The imaginary part of the output will be: + double[][] expectedIm = new double[][] + { + new double[] { 0, -8, 0, 8 }, + new double[] { -0.8660, -4.5980, 6.0621, -0.5980 }, + new double[] { 0.8660, 0.5980, -6.0621, 4.5980 } + }; + + // We can recover the initial matrix using the backward transform: + FourierTransform2.FFT2(input, FourierTransform.Direction.Backward); + #endregion + + string a = re.ToCSharp(); + string b = im.ToCSharp(); + + Assert.IsTrue(re.IsEqual(expectedRe, 1e-3)); + Assert.IsTrue(im.IsEqual(expectedIm, 1e-3)); + var expected = new Complex[][] + { + new Complex[] { 1, 4, 8, 0 }, + new Complex[] { 0, 4, 5, 0 }, + new Complex[] { 0, 0, 8, 0 }, + }; + + Assert.IsTrue(input.Re().IsEqual(expected.Re(), 1e-3)); + Assert.IsTrue(input.Im().IsEqual(expected.Im(), 1e-3)); + } + + [Test] public void gh878() {