diff --git a/src/FftSharp.Tests/Inverse.cs b/src/FftSharp.Tests/Inverse.cs index 290ef64..a0befd4 100644 --- a/src/FftSharp.Tests/Inverse.cs +++ b/src/FftSharp.Tests/Inverse.cs @@ -51,5 +51,24 @@ public void Test_IFFT_MatchesOriginal() Assert.AreEqual(original[i].Imaginary, ifft[i].Imaginary, 1e-6); } } + + [Test] + public void Test_InverseReal_MatchesOriginal() + { + Random rand = new Random(0); + double[] original = new double[1024]; + for (int i = 0; i < original.Length; i++) + original[i] = rand.NextDouble() - .5; + + System.Numerics.Complex[] rfft = FftSharp.FFT.ForwardReal(original); + double[] irfft = FftSharp.FFT.InverseReal(rfft); + + Assert.AreEqual(original.Length, irfft.Length); + + for (int i = 0; i < irfft.Length; i++) + { + Assert.AreEqual(original[i], irfft[i], 1e-10); + } + } } } diff --git a/src/FftSharp/FFT.cs b/src/FftSharp/FFT.cs index 42f2755..46e46c5 100644 --- a/src/FftSharp/FFT.cs +++ b/src/FftSharp/FFT.cs @@ -1,4 +1,4 @@ -using System; +using System; namespace FftSharp; @@ -73,6 +73,37 @@ public static System.Numerics.Complex[] ForwardReal(double[] samples) return real; } + /// + /// Calculate IFFT and return just the real component of the spectrum + /// + public static double[] InverseReal(System.Numerics.Complex[] spectrum) + { + int spectrumSize = spectrum.Length; + int windowSize = (spectrumSize - 1) * 2; + + if (!FftOperations.IsPowerOfTwo(windowSize)) + throw new ArgumentException($"{nameof(spectrum)} length must be a power of two plus 1"); + + System.Numerics.Complex[] buffer = new System.Numerics.Complex[windowSize]; + + for (int i = 0; i < spectrumSize; i++) + buffer[i] = spectrum[i]; + + for (int i = spectrumSize; i < windowSize; i++) + { + int iMirrored = spectrumSize - (i - (spectrumSize - 2)); + buffer[i] = System.Numerics.Complex.Conjugate(spectrum[iMirrored]); + } + + Inverse(buffer); + + double[] result = new double[windowSize]; + for (int i = 0; i < windowSize; i++) + result[i] = buffer[i].Real; + + return result; + } + /// /// Apply the inverse Fast Fourier Transform (iFFT) to the Complex array in place. /// Length of the array must be a power of 2.