#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();
}
}
}
}
}
}