Search code examples
c#fftnaudiomathnet-numerics

What is wrong with my Fourier Transformation (Convolution) in MathNet.Numerics? - C#

I am trying to do a simple Convolution between 2 audio files using the MathNet.Numerics's FFT (Fast Fourier Transformation), but I get some weird background sounds, after the IFFT.

I tested if it's the Convolution or the Transformations, thats causing the problem, and I found out that the problem shows already in the FFT -> IFFT (Inverze FFT) conversion.

My code for a simple FFT and IFFT:

float[] sound; //here are stored my samples

Complex[] complexInput = new Complex[sound.Length];
for (int i = 0; i < complexInput.Length; i++)
{
      Complex tmp = new Complex(sound[i],0);
      complexInput[i] = tmp;
 }

MathNet.Numerics.IntegralTransforms.Fourier.Forward(complexInput);

//do some stuff

MathNet.Numerics.IntegralTransforms.Fourier.Inverse(complexInput);

float[] outSamples = new float[complexInput.Length];

for (int i = 0; i < outSamples.Length; i++)
     outSamples[i] = (float)complexInput[i].Real;

After this, the outSamples are corrupted with some wierd background sound/noise, even though I'm not doing anything between the FFT and IFFT.

What am I missing?

Solution

  • The current implementation of MathNet.Numerics.IntegralTransform.Fourier (see Fourier.cs and Fourier.Bluestein.cs) uses Bluestein's algorithm for any FFT lengths that are not powers of 2.

    This algorithm involves the creation of a Bluestein sequence (which includes terms proportional to n2) which up to version 3.6.0 was using the following code:

    static Complex[] BluesteinSequence(int n)
{
  double s = Constants.Pi/n;
  var sequence = new Complex[n];

  for (int k = 0; k < sequence.Length; k++)
  {
    double t = s*(k*k); // <--------------------- (k*k) 32-bit int expression
    sequence[k] = new Complex(Math.Cos(t), Math.Sin(t));
  }

  return sequence;
}

    For any size n greater than 46341, the intermediate expression (k*k) in this implementation is computed using int arithmetic (a 32-bit type as per MSDN integral type reference table) which results in numeric overflows for the largest values of k. As such the current implementation of MathNet.Numerics.IntegralTransfom.Fourier only supports input array sizes which are either powers of 2 or non-powers of 2 up to 46341 (included).

    Thus for large input arrays, a workaround could be to pad your input to the next power of 2.

    Note: this observation is based on version 3.6.0 of MathNet.Numerics, although the limitation appears to have been present in earlier releases (the Bluestein sequence code has not changed significantly going as far back as version 2.1.1).

    Update 2015/04/26:

    After I posted this and a comment on an similar issue on github bugtracking, the issue was quickly fixed in MathNet.Numerics. The fix should now be available in version 3.7.0. Note however that you may still want to pad to a power of two for performance reasons, especially since you already need to zero pad for the linear convolution.