Fast Fourier Transform
This time I want to talk about Fast Fourier Transform functions and how to use these in C# properly.
Intel MKL team
work mainly with C++ and Fortran. They didn’t dedicate a lot of effort toward
C# so sometimes if you work in C# and you want to use MKL you face some
difficulties.
As I said
in the previous article about MKL, the first difficulty is how to declare the
functions in C# code. For this I suggest a document from Intel: 33774-Using_Intel_MKL_IPP_NET_Framework.pdf.
You can find it on line. It is not only a problem of declaration, but also a
question of how to interpret the error messages, sometimes you don’t know
what’s going on.
For example
if you wrong something, you receive messages like these:
- "Inconsistent configuration parameter"
- "Functionality is not implemented".
Where is a
tutorial that explain what I have wronged? You have to contact the Intel
Support and be patient.
With a
little help you can understand that the first error is because you are
assigning a wrong value to a variable (DftiSetValue). The second one is because
you had configured a “Descriptor” with the wrong parameters (DftiCreateDescriptor).
Then after
completed a declarative part you can test something, for example I used this
web site (http://www.sccon.ca/sccon/fft/fft3.htm) with the FFT input and output
certified, as reference.
For some
example of FFT usage this is the page: https://software.intel.com/en-us/node/471390.
The problem
is that, those example are written in Fortran or C++ and contains some little
imprecisions that can confuse you.
Let
consider the paragraph: “One-dimensional Out-of-place FFT”. Where they perform
a real to complex conjugate-even transform.
They use a
real vector for the input and a real vector for the output. But the output is a
complex vector! Also the vector size was wrong. These things and other little
imprecisions make me lost a lot of time. But not desperate because I will
provide you what you need to use FFT in C#.
There are
some possible versions of DftiComputeForward() function that compute the FF
Transform so at a beginning I adopted this one:
status =
DftiComputeForward ( desc_handle, xre_in, xim_in , yre_out, yim_out )
It allows 2
real vectors in input and two real vectors for the output. These vectors represents
real and imaginary part.
By the way,
although I didn’t receive any error from the method execution the result was always
zero! Again problems!
But as I
always say: “it ain't over til it's over” I have tried another version of the
function, this time with a input a complex vector, and output a complex vector.
And like a magic this time worked.
status =
DftiComputeForward( desc_handle, x_in, y_out )
So for your
joy here the code:
using System;
using System.Numerics;
using System.Runtime.InteropServices;
using System.Security;
namespace mkl
{
[SuppressUnmanagedCodeSecurity]
public unsafe class MKLWrapper
{
//==================================================
#region FFT CONSTANTS
enum DFTI_CONFIG_PARAM
{
/* Domain for forward transform. No default value */
DFTI_FORWARD_DOMAIN = 0,
/* Dimensionality, or rank. No default value */
DFTI_DIMENSION = 1,
/* Length(s) of transform. No default value */
DFTI_LENGTHS = 2,
/* Floating point precision. No default value */
DFTI_PRECISION = 3,
/* Scale factor for forward transform [1.0] */
DFTI_FORWARD_SCALE = 4,
/* Scale factor for backward transform [1.0] */
DFTI_BACKWARD_SCALE = 5,
/* Exponent sign for forward transform [DFTI_NEGATIVE] */
/* DFTI_FORWARD_SIGN = 6, ## NOT IMPLEMENTED */
/* Number of data sets to be transformed [1] */
DFTI_NUMBER_OF_TRANSFORMS = 7,
/* Storage of finite complex-valued sequences in complex domain
[DFTI_COMPLEX_COMPLEX] */
DFTI_COMPLEX_STORAGE = 8,
/* Storage of finite real-valued sequences in real domain
[DFTI_REAL_REAL] */
DFTI_REAL_STORAGE = 9,
/* Storage of finite complex-valued sequences in conjugate-even
domain [DFTI_COMPLEX_REAL] */
DFTI_CONJUGATE_EVEN_STORAGE = 10,
/* Placement of result [DFTI_INPLACE] */
DFTI_PLACEMENT = 11,
/* Generalized strides for input data layout [tigth, row-major for
C] */
DFTI_INPUT_STRIDES = 12,
/* Generalized strides for output data layout [tight, row-major
for C] */
DFTI_OUTPUT_STRIDES = 13,
/* Distance between first input elements for multiple transforms
[0] */
DFTI_INPUT_DISTANCE = 14,
/* Distance between first output elements for multiple transforms
[0] */
DFTI_OUTPUT_DISTANCE = 15,
/* Effort spent in initialization [DFTI_MEDIUM] */
/* DFTI_INITIALIZATION_EFFORT = 16, ## NOT IMPLEMENTED */
/* Use of workspace during computation [DFTI_ALLOW] */
DFTI_WORKSPACE = 17,
/* Ordering of the result [DFTI_ORDERED] */
DFTI_ORDERING = 18,
/* Possible transposition of result [DFTI_NONE] */
DFTI_TRANSPOSE = 19,
/* User-settable descriptor name [""] */
DFTI_DESCRIPTOR_NAME = 20, /* DEPRECATED */
/* Packing format for DFTI_COMPLEX_REAL storage of finite
conjugate-even sequences [DFTI_CCS_FORMAT] */
DFTI_PACKED_FORMAT = 21,
/* Commit status of the descriptor - R/O parameter */
DFTI_COMMIT_STATUS = 22,
/* Version string for this DFTI implementation - R/O parameter */
DFTI_VERSION = 23,
/* Ordering of the forward transform - R/O parameter */
/* DFTI_FORWARD_ORDERING = 24, ## NOT IMPLEMENTED */
/* Ordering of the backward transform - R/O parameter */
/* DFTI_BACKWARD_ORDERING = 25, ## NOT IMPLEMENTED */
/* Number of user threads that share the descriptor [1] */
DFTI_NUMBER_OF_USER_THREADS = 26,
/* Limit the number of threads used by this descriptor [0 = don't care] */
DFTI_THREAD_LIMIT = 27,
/* Possible input data destruction [DFTI_AVOID = prevent input data]*/
DFTI_DESTROY_INPUT = 28
};
enum DFTI_CONFIG_VALUE
{
/* DFTI_COMMIT_STATUS */
DFTI_COMMITTED = 30,
DFTI_UNCOMMITTED = 31,
/* DFTI_FORWARD_DOMAIN */
DFTI_COMPLEX = 32,
DFTI_REAL = 33,
/* DFTI_CONJUGATE_EVEN = 34, ## NOT IMPLEMENTED */
/* DFTI_PRECISION */
DFTI_SINGLE = 35,
DFTI_DOUBLE = 36,
/* DFTI_FORWARD_SIGN */
/* DFTI_NEGATIVE = 37, ## NOT IMPLEMENTED */
/* DFTI_POSITIVE = 38, ## NOT IMPLEMENTED */
/* DFTI_COMPLEX_STORAGE and DFTI_CONJUGATE_EVEN_STORAGE */
DFTI_COMPLEX_COMPLEX = 39,
DFTI_COMPLEX_REAL = 40,
/* DFTI_REAL_STORAGE */
DFTI_REAL_COMPLEX = 41,
DFTI_REAL_REAL = 42,
/* DFTI_PLACEMENT */
DFTI_INPLACE = 43, /* Result overwrites input */
DFTI_NOT_INPLACE = 44, /* Have another place for result */
/* DFTI_INITIALIZATION_EFFORT */
/* DFTI_LOW = 45, ## NOT IMPLEMENTED */
/* DFTI_MEDIUM = 46, ## NOT IMPLEMENTED */
/* DFTI_HIGH = 47, ## NOT IMPLEMENTED */
/* DFTI_ORDERING */
DFTI_ORDERED = 48,
DFTI_BACKWARD_SCRAMBLED = 49,
/* DFTI_FORWARD_SCRAMBLED = 50, ## NOT IMPLEMENTED */
/* Allow/avoid certain usages */
DFTI_ALLOW = 51, /* Allow transposition or workspace */
DFTI_AVOID = 52,
DFTI_NONE = 53,
/* DFTI_PACKED_FORMAT (for storing congugate-even finite sequence
in real array) */
DFTI_CCS_FORMAT = 54, /* Complex conjugate-symmetric */
DFTI_PACK_FORMAT = 55, /* Pack format for real DFT */
DFTI_PERM_FORMAT = 56, /* Perm format for real DFT */
DFTI_CCE_FORMAT = 57 /* Complex conjugate-even */
};
#endregion
private MKLWrapper()
{
}
//==================================================
#region FFT (DECLARATIONS)
//DFTI_EXTERN char* DftiErrorMessage(MKL_LONG);
[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern IntPtr DftiErrorMessage(
long error_code
);
//DFTI_EXTERN MKL_LONG DftiCreateDescriptor(DFTI_DESCRIPTOR_HANDLE*,
// enum DFTI_CONFIG_VALUE, /* precision */
// enum DFTI_CONFIG_VALUE, /* domain */
// MKL_LONG, ...);
[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int DftiCreateDescriptor(
ref IntPtr handle_descriptor,
int precision,
int domain,
int dime,
int size
);
//DFTI_EXTERN MKL_LONG DftiCommitDescriptor(DFTI_DESCRIPTOR_HANDLE);
[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int DftiCommitDescriptor(
IntPtr handle_descriptor
);
//DFTI_EXTERN MKL_LONG DftiSetValue(DFTI_DESCRIPTOR_HANDLE, enum DFTI_CONFIG_PARAM, ...);
[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int DftiSetValue(IntPtr handle_descriptor, int config_param, __arglist);
//DFTI_EXTERN MKL_LONG DftiComputeBackward(DFTI_DESCRIPTOR_HANDLE, void*, ...);
[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int DftiComputeBackward(
IntPtr handle_descriptor,
[In] Complex[] x_in,
[Out] Complex[] x_out
);
[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int DftiComputeForward(
IntPtr handle_descriptor,
[In] Complex[] x_in,
[Out] Complex[] x_out
);
//DFTI_EXTERN MKL_LONG DftiFreeDescriptor(DFTI_DESCRIPTOR_HANDLE*);
[DllImport("mkl_rt.dll", ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
internal static extern int DftiFreeDescriptor(
ref IntPtr handle_descriptor
);
#endregion
//==================================================
#region FFT (METHODS)
public static string MKL_FFT_Error(int error, string sText = "")
{
if (error == 0)
return "";
IntPtr s = DftiErrorMessage(error);
string sMsg = Marshal.PtrToStringAnsi(s);
throw new Exception(sText + sMsg);
}
public static Complex[] MKL_FFT(Complex[] in_matrix, double fScalingFactor = 1)
{
int status;
IntPtr des_handle = IntPtr.Zero;
Complex[] out_matrix = new Complex[in_matrix.Length];
try
{
status = DftiCreateDescriptor(ref des_handle, (int)DFTI_CONFIG_VALUE.DFTI_DOUBLE, (int)DFTI_CONFIG_VALUE.DFTI_COMPLEX, 1, in_matrix.Length);
MKL_FFT_Error(status, "DftiCreateDescriptor() ");
status = DftiSetValue(des_handle, (int)DFTI_CONFIG_PARAM.DFTI_FORWARD_SCALE, __arglist(fScalingFactor));
MKL_FFT_Error(status, "DftiSetValue() ");
status = DftiSetValue(des_handle, (int)DFTI_CONFIG_PARAM.DFTI_PLACEMENT, __arglist((int)DFTI_CONFIG_VALUE.DFTI_NOT_INPLACE));
MKL_FFT_Error(status, "DftiSetValue() ");
status = DftiCommitDescriptor(des_handle); //Finalize the descriptor
MKL_FFT_Error(status, "DftiCommitDescriptor() ");
status = DftiComputeForward(des_handle, in_matrix, out_matrix); //Compute the Backward FFT
MKL_FFT_Error(status, "DftiComputeForward() ");
status = DftiFreeDescriptor(ref des_handle); //Free the descriptor
MKL_FFT_Error(status, "DftiFreeDescriptor() ");
return out_matrix;
}
catch (Exception ex)
{
if (des_handle != IntPtr.Zero)
DftiFreeDescriptor(ref des_handle); //Free the descriptor
throw ex;
}
}
public static Complex[] MKL_FFT_Inv(Complex[] in_matrix, double fScalingFactor = 1)
{
int status;
IntPtr des_handle = IntPtr.Zero;
Complex[] out_matrix = new Complex[in_matrix.Length];
try
{
status = DftiCreateDescriptor(ref des_handle, (int)DFTI_CONFIG_VALUE.DFTI_DOUBLE, (int)DFTI_CONFIG_VALUE.DFTI_COMPLEX, 1, in_matrix.Length);
MKL_FFT_Error(status, "DftiCreateDescriptor() ");
status = DftiSetValue(des_handle, (int)DFTI_CONFIG_PARAM.DFTI_BACKWARD_SCALE, __arglist(fScalingFactor));
MKL_FFT_Error(status, "DftiSetValue() ");
status = DftiSetValue(des_handle, (int)DFTI_CONFIG_PARAM.DFTI_PLACEMENT, __arglist((int)DFTI_CONFIG_VALUE.DFTI_NOT_INPLACE));
MKL_FFT_Error(status, "DftiSetValue() ");
status = DftiCommitDescriptor(des_handle); //Finalize the descriptor
MKL_FFT_Error(status, "DftiCommitDescriptor() ");
status = DftiComputeBackward(des_handle, in_matrix, out_matrix); //Compute the Backward FFT
MKL_FFT_Error(status, "DftiComputeBackward() ");
status = DftiFreeDescriptor(ref des_handle); //Free the descriptor
return out_matrix;
}
catch (Exception ex)
{
if (des_handle != IntPtr.Zero)
DftiFreeDescriptor(ref des_handle); //Free the descriptor
throw ex;
}
}
#endregion
}
}
Here the example:
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using mkl;
using System.Numerics;
namespace FFTTest
{
class Program
{
#region TEST INPUT MATRIX
//http://www.sccon.ca/sccon/fft/fft3.htm
//INPUT
//REAL[0] = 1.000, IMAG[0] = 0.000
//REAL[1] = 0.000, IMAG[1] = 0.000
//REAL[2] = 0.000, IMAG[2] = 0.000
//REAL[3] = 0.000, IMAG[3] = 0.000
//REAL[4] = 0.000, IMAG[4] = 0.000
//REAL[5] = 0.000, IMAG[5] = 0.000
//REAL[6] = 0.000, IMAG[6] = 0.000
//REAL[7] = 0.000, IMAG[7] = 0.000
//RESULTS
//REAL[0] = 0.125, IMAG[0] = 0.000
//REAL[1] = 0.125, IMAG[1] = 0.000
//REAL[2] = 0.125, IMAG[2] = 0.000
//REAL[3] = 0.125, IMAG[3] = 0.000
//REAL[4] = 0.125, IMAG[4] = 0.000
//REAL[5] = 0.125, IMAG[5] = 0.000
//REAL[6] = 0.125, IMAG[6] = 0.000
//REAL[7] = 0.125, IMAG[7] = 0.000
#endregion
static void Main(string[] args)
{
Complex[] A = new Complex[8];
Complex[] B = null;
Complex[] C = null;
A[0] = new Complex(0, 0);
A[1] = new Complex(1, 0);
A[2] = new Complex(0, 0);
A[3] = new Complex(0, 0);
A[4] = new Complex(0, 0);
A[5] = new Complex(0, 0);
A[6] = new Complex(0, 0);
A[7] = new Complex(0, 0);
try
{
DebugMatrix(A, "A:");
B = MKLWrapper.MKL_FFT(A, 1.0f/ A.Length);
DebugMatrix(B, "FFT() B:");
C = MKLWrapper.MKL_FFT_Inv(B); //FFT Backward
DebugMatrix(C, "iFFT() C:");
//C = MKLWrapper.MKL_FFT(B);
//DebugMatrix(lstResults, C, "A:");
//MessageBox.Show("Operation completed succesfully!");
}
catch (Exception ex)
{
Console.WriteLine("Error: " + ex.Message);
}
Console.ReadLine();
}
public static void DebugMatrix(Complex[] a, string sCaption, int iRowLimit = 10)
{
if (a == null) return;
string sLine = "";
Console.WriteLine("");
Console.WriteLine(sCaption);
int iDime1 = a.GetUpperBound(0);
int iStart1 = a.GetLowerBound(0);
Console.WriteLine("Size: " + a.GetLength(0).ToString());
for (int i = iStart1; i <= iDime1; i++)
{
sLine = i.ToString() + "> ";
sLine += a[i].ToString("0.###E+000") + "\t";
//sLine += a[i, j].ToString("") + "\t";
Console.WriteLine(sLine);
if (i > iRowLimit)
break;
}
Console.WriteLine("");
Console.WriteLine("");
}
}
}
Part 1 | Part 2 | Part 3 | Part 4
No comments :
Post a Comment