#region Copyright and License /* This file is part of FFTW.NET, a wrapper around the FFTW library for the .NET framework. Copyright (C) 2017 Tobias Meyer This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. */ #endregion using System; using System.Text; using System.Numerics; using System.Runtime.InteropServices; using System.Linq; namespace FFTW.NET { /// /// This class provides methods for convenience. /// However, for optimal performance you should consider using /// or directly. /// public static class DFT { /// /// Large object heap threshold in bytes /// const int LohThreshold = 85000; /// /// Memory alignment in bytes /// const int MemoryAlignment = 16; static readonly BufferPool _bufferPool = new BufferPool(LohThreshold); /// /// Performs a complex-to-complex fast fourier transformation. The dimension is inferred from the input (). /// /// /// public static void FFT(IPinnedArray input, IPinnedArray output, PlannerFlags plannerFlags = PlannerFlags.Default, int nThreads = 1) => Transform(input, output, DftDirection.Forwards, plannerFlags, nThreads); /// /// Performs a complex-to-complex inverse fast fourier transformation. The dimension is inferred from the input (). /// /// /// public static void IFFT(IPinnedArray input, IPinnedArray output, PlannerFlags plannerFlags = PlannerFlags.Default, int nThreads = 1) => Transform(input, output, DftDirection.Backwards, plannerFlags, nThreads); static void Transform(IPinnedArray input, IPinnedArray output, DftDirection direction, PlannerFlags plannerFlags, int nThreads) { if ((plannerFlags & PlannerFlags.Estimate) == PlannerFlags.Estimate) { using (var plan = FftwPlanC2C.Create(input, output, direction, plannerFlags, nThreads)) { plan.Execute(); return; } } using (var plan = FftwPlanC2C.Create(input, output, direction, plannerFlags | PlannerFlags.WisdomOnly, nThreads)) { if (plan != null) { plan.Execute(); return; } } /// If with no plan can be created /// and is not specified, we use /// a different buffer to avoid overwriting the input if (input != output) { using (var plan = FftwPlanC2C.Create(output, output, input.Rank, input.GetSize(), direction, plannerFlags, nThreads)) { input.CopyTo(output); plan.Execute(); } } else { using (var bufferContainer = _bufferPool.RequestBuffer(input.Length * Marshal.SizeOf() + MemoryAlignment)) using (var buffer = new AlignedArrayComplex(bufferContainer.Buffer!, MemoryAlignment, input.GetSize())) using (var plan = FftwPlanC2C.Create(buffer, buffer, input.Rank, input.GetSize(), direction, plannerFlags, nThreads)) { input.CopyTo(plan.Input); plan.Execute(); plan.Output.CopyTo(output, 0, 0, input.Length); } } } /// /// Performs a real-to-complex fast fourier transformation. /// /// /// public static void FFT(IPinnedArray input, IPinnedArray output, PlannerFlags plannerFlags = PlannerFlags.Default, int nThreads = 1) { if ((plannerFlags & PlannerFlags.Estimate) == PlannerFlags.Estimate) { using (var plan = FftwPlanRC.Create(input, output, DftDirection.Forwards, plannerFlags, nThreads)) { plan.Execute(); return; } } using (var plan = FftwPlanRC.Create(input, output, DftDirection.Forwards, plannerFlags | PlannerFlags.WisdomOnly, nThreads)) { if (plan != null) { plan.Execute(); return; } } /// If with no plan can be created /// and is not specified, we use /// a different buffer to avoid overwriting the input using (var bufferContainer = _bufferPool.RequestBuffer(input.Length * sizeof(double) + MemoryAlignment)) using (var buffer = new AlignedArrayDouble(bufferContainer.Buffer!, MemoryAlignment, input.GetSize())) using (var plan = FftwPlanRC.Create(buffer, output, DftDirection.Forwards, plannerFlags, nThreads)) { input.CopyTo(plan.BufferReal); plan.Execute(); } } public static Complex[] IFFT(double[] data, double[] img) { if (data == null || data.Length == 0 || img == null || img.Length != data.Length) return new Complex[0]; var input = new PinnedArray(Enumerable.Range(0, data.Length).Select(x => new Complex(data[x], img[x])).ToArray()); var output = new PinnedArray(new Complex[data.Length]); IFFT(input, output); var result = output.Buffer.Cast().ToArray(); input.Dispose(); output.Dispose(); return result; } public static Complex[] FFT(double[] data, double[] img) { if (data == null || data.Length == 0 || img == null || img.Length != data.Length) return new Complex[0]; var input = new PinnedArray(Enumerable.Range(0, data.Length).Select(x => new Complex(data[x], img[x])).ToArray()); var output = new PinnedArray(new Complex[data.Length]); FFT(input, output); var result = output.Buffer.Cast().ToArray(); input.Dispose(); output.Dispose(); return result; } /// /// Performs a complex-to-real inverse fast fourier transformation. /// /// /// public static void IFFT(IPinnedArray input, IPinnedArray output, PlannerFlags plannerFlags = PlannerFlags.Default, int nThreads = 1) { if ((plannerFlags & PlannerFlags.Estimate) == PlannerFlags.Estimate) { using (var plan = FftwPlanRC.Create(output, input, DftDirection.Backwards, plannerFlags, nThreads)) { plan.Execute(); return; } } using (var plan = FftwPlanRC.Create(output, input, DftDirection.Backwards, plannerFlags | PlannerFlags.WisdomOnly, nThreads)) { if (plan != null) { plan.Execute(); return; } } /// If with no plan can be created /// and is not specified, we use /// a different buffer to avoid overwriting the input using (var bufferContainer = _bufferPool.RequestBuffer(input.Length * Marshal.SizeOf() + MemoryAlignment)) using (var buffer = new AlignedArrayComplex(bufferContainer.Buffer!, MemoryAlignment, input.GetSize())) using (var plan = FftwPlanRC.Create(output, buffer, DftDirection.Backwards, plannerFlags, nThreads)) { input.CopyTo(plan.BufferComplex); plan.Execute(); } } /// /// Gets the required size of the complex buffer in a complex-to-real /// or rea-to-complex transormation from the size of the real buffer. /// /// public static int[] GetComplexBufferSize(int[] realBufferSize) { int[] n = new int[realBufferSize.Length]; Buffer.BlockCopy(realBufferSize, 0, n, 0, n.Length * sizeof(int)); n[n.Length - 1] = realBufferSize[n.Length - 1] / 2 + 1; return n; } /// /// Provides access to FFTW's wisdom mechanism /// /// public static class Wisdom { /// /// Exports the accumulated wisdom to a file. /// /// /// public static bool Export(string filename) { lock (FftwInterop.Lock) { return FftwInterop.GetDelegate()(filename); } } /// /// Imports wisdom from a file. The Current accumulated wisdom is replaced. /// Wisdom is hardware specific, thus importing wisdom created with different hardware /// can result in sub-optimal plans and should not be done. /// /// public static bool Import(string filename) { lock (FftwInterop.Lock) { return FftwInterop.GetDelegate()(filename); } } /// /// Clears the current wisdom. /// /// public static void Clear() { lock (FftwInterop.Lock) { FftwInterop.GetDelegate()(); } } /// /// Gets or sets the current wisdom. /// public static string Current { get { lock (FftwInterop.Lock) { return FftwInterop.fftw_export_wisdom_to_string(); } } set { lock (FftwInterop.Lock) { if (string.IsNullOrEmpty(value)) FftwInterop.GetDelegate()(); else if (!FftwInterop.GetDelegate()(value)) throw new FormatException(); } } } } } }