This project has moved and is read-only. For the latest updates, please go here.

Iridium: Fourier Transforms

Math.NET Iridium supports fast Fourier transformations (FFT) in both real and complex domains. Multidimensional scenarios are also supported (e.g. to transform an image). However, note that the Fourier transformation naturally lives in the complex domain; the transformation of a real function in time space always returns a complex function in frequency space. Apart from optimization, in the end all real transformations are mapped to complex transformations (although some tricks based on the symmetries below allow us for example to transform two real signals with a single complex transformation).

For more details about the Fourier transformation have a look at Wikipedia:

Symmetries between the time and the frequency space

Time Space h(t) (formal) <=> Frequency Space H(f) (formal)
real Im{h} = 0 <=> Re{H} part even, Im{H} odd H(-f) = (H(f))*
imaginary Re{h} = 0 <=> Re{H} part odd, Im{H} even H(-f) = -(H(f))*
even h(-t) = h(t) <=> even H(-f) = H(f)
odd h(-t) = -h(t) <=> odd H(-f) = -H(f)
real even Im{h} = 0, h(-t) = h(t) <=> real even Im{H} = 0, H(-f) = H(f)
real odd Im{h} = 0, h(-t) = -h(t) <=> imaginary odd Re{H} = 0, H(-f) = -H(f)
imaginary even Re{h} = 0, h(-t) = h(t) <=> imaginary even Re{H} = 0, H(-f) = H(f)
imaginary odd Re{h} = 0, h(-t) = -h(t) <=> real odd Im{H} = 0, H(-f) = -H(f)

Transformation Conventions

There are various definitions around for what a discrete Fourier transformation and its inverse shall look like. They all agree on the general formula, but variate slightly on scaling factors and the sign of exponent. Iridium offers flags to specify the convention you expect, including some templates:
  • Default: Uses a negative exponent sign in forward transformations, and symmetric scaling (that is, sqrt(1/N) for both forward and inverse transformation). This is the convention used in Maple and is widely accepted in the educational sector (due to the symmetry).
  • AsymmetricScaling: Set this flag to suppress scaling on the forward transformation but scale the inverse transform with 1/N.
  • NoScaling: Set this flag to suppress scaling for both forward and inverse transformation. Note that in this case if you apply first the forward and then inverse transformation you won't get back the original signal (by factor N/2).
  • InverseExponent: Uses the positive instead of the negative sign in the forward exponent, and the negative (instead of positive) exponent in the inverse transformation.
  • Matlab: Use this flag if you need Matlab compatibility. Equals to setting the AsymmetricScaling flag.
  • NumericalRecpies: Use this flag if you need Numerical Recipes compatibility. Equal to setting both the InverseExponent and the NoScaling flags.

Data Structures

Iridium does not yet implement any prime factor subdivision FFT algorithm and thus does not support samples of arbitrary size. All samples sets have to have a length of a power of two (2,4,8,16,...). You may try padding sets of other length with zero or with the sample mean, but remember that this influences the frequency space, and I recommend to avoid it where possible. Clearly the medium-time target of Iridium is to provide and implementation that supports sets of arbitrary length, but unfortunately we can't guarantee a point in time. Maybe we'll provide a wrapper to FFTW (probably the best implementation these days, but not managed code/.Net, and unfortunately not LGPL-compliant) in the meantime.

The ordering is equivalent to that of Numerical Recipes:

Complex numbers may be represented either as MathNet.Numerics.Complex instances, or as pairs of doubles. In the second case, they're always ordered in pairs: first the real, and then imaginary value. So the final sample set is alternating: Re{z1}, Im{z1}, Re{z2}, Im{z2}, Re{z3}, Im{z3}, ...

In multiple dimensions, the data is always (expected to be) ordered such that the last index (the last dimension) changes most rapidly. If you have an image with two dimensions (y,x) and is indexed it with [y,x], this means the array has to be stored row-by-row, so the x coordinate changes more rapidly than the y coordinate.

In frequency space Iridium always returns the whole spectrum (both positive and negative part, even if it's redundant) in an array of length 2*N, where the negative parts are "wrapped" around. The DC zero frequency component is at index 0/1 (real/imag). The most positive and most negative frequency component are always equal and at index N/N+1. The smallest positive frequency is at 2/3, the smallest (in magnitude) negative frequency at 2*N-2/2*N-1. Hence, if we have frequencies -40,-30,..,30,40, they're stored in the following order (r=real,i=imaginary): (0r,0i,10r,10i,20r,20i,30r,30i,40r=-40r,40i=-40i,-30r,-30i,-20r,-20i,-10r,-10i)

You may generate the expected scales for both frequency and time space using the two helper methods GenerateTimeScale(f,N) and GenerateFrequencyScale(f,N). They return an array with the expected time/frequency of each sample. This might also be helpful to generate the axis of 2D plots.

Complex Fourier Transformation

In the following we'll give a short introduction how to compute complex Fourier transformations with Iridium.

FFT of a real signal (packed in a complex signal)

Of course we may also use the complex version to transform a real signal.

We first generate a real even signal. In Matlab this could be done like this: samples_t = 1.0 ./ (([-16:1:15] ./ 16) .^ 2 + 1.0)

int numSamples = 32;
int length = 2 * numSamples;
double[] data = new double[length];
for(int i = 0; i < length; i += 2)
{
    double z = (double)(i - numSamples) / numSamples;
    data[i] = 1.0 / (z * z + 1.0);  // real part
    data[i + 1] = 0.0;              // imaginary part
}

Result (only real values): 0.5000 0.5322 0.5664 0.6024 0.6400 ... 0.6400 0.6024 0.5664 0.5322

Then we can simply create complex transformation and use the TransformForward method. Note that if you do multiple transformations it is recommended to reuse the transformation class instance, as it caches some precomputations locally. I chose Matlab convention to be able to compare the result with Matlab (generate it with: samples_f = fft(samples_t).

ComplexFourierTransformation cft = new ComplexFourierTransformation(TransformationConvention.Matlab);
cft.TransformForward(data);

Result (only real values): 25.1275 -3.6230 -0.3105 -0.1888 -0.1064 ... -0.1064 -0.1888 -0.3105 -3.6230

Note that the result is again real and even as expected.

The inverse transformation is just as simple:

cft.TransformBackward(data);

FFT of a complex signal

Works exactly like above, except you also add non-zero values to the imaginary part. For example you may want to try the following sample set with (the same) even real component but also with an odd imaginary component. Of course these components don't have to be even or odd, they can be arbitrary.

Matlab:{BR}
samples_t = 1.0i * (<nowiki>[-16:1:15]</nowiki> ./ 16) ./ (([-16:1:15] ./ 16) .^ 2 + 1.0)
 + 1.0 ./ ((<nowiki>[-16:1:15]</nowiki> ./ 16) .^ 2 + 1.0)
samples_t(1) = real(samples_t(1))
samples_f = fft(samples_t)


int numSamples = 32;
int length = 2 * numSamples;
double[] data = new double[length];
for(int i = 0; i < length; i += 2)
{
    double z = (double)(i - numSamples) / numSamples;
    data[i] = 1.0 / (z * z + 1.0);
    data[i + 1] = z / (z * z + 1.0);
}
data[1] = 0.0; // peridoic continuation; force odd

cft.Convention = TransformationConvention.Matlab;
cft.TransformForward(data);

Real Fourier Transformation

When you only work with real signals it may be easier and more performant to work directly with the real version of the transformation. There's a drawback though: Since the result of a real transformation is still complex we can't operate inplace directly on the array. Instead, the result is provided in two output parameter arrays, one containing the real and one the imaginary part.

FFT of a single Real Signal

To use the real version of the transformation, use the RealFourierTransformation class:

int numSamples = 32;
int half = numSamples >> 1;
double[] data = new double[numSamples];
for(int i = 0; i < numSamples; i++)
{
    double z = (double)(i - half) / half;
    data[i] = 1.0 / (z * z + 1.0);
}

double[] freqReal, freqImag;
RealFourierTransformation rft = new RealFourierTransformation(); // default convention
rft.TransformForward(data, out freqReal, out freqImag);

And the inverse transformation:

double[] data2
rft.TransformBackward(freqReal, freqImag, out data2);

FFT of two Real Signals

Thanks to a trick you can also transform two real signals of the same length with a single complex transformation. We again use an even and an odd signal, but both signals could have an arbitrary form:

int numSamples = 32;
int half = numSamples >> 1;

double[] dataEven = new double[numSamples];
double[] dataOdd = new double[numSamples];

for(int i = 0; i < numSamples; i++)
{
    double z = (double)(i - half) / half;
    dataEven[i] = 1.0 / (z * z + 1.0);
    dataOdd[i] = z / (z * z + 1.0);
}
dataOdd[0] = 0.0; // peridoic continuation; force odd

double[] evenReal, evenImag, oddReal, oddImag;

// Forward Transform
rft.Convention = TransformationConvention.Default;
rft.TransformForward(dataEven, dataOdd,
    out evenReal, out evenImag, out oddReal, out oddImag);

// Backward Transform
double[] dataEven2, dataOdd2;
rft.TransformBackward(evenReal, evenImag, oddReal, oddImag,
    out dataEven2, out dataOdd2);

Last edited Jul 29, 2012 at 4:43 PM by cdrnet, version 2

Comments

No comments yet.