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

Question

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?

Shot in the dark - have you tried taking the absolute value of the complex value, instead of slicing off just the real component? — antiduh, Apr 15 '15 at 20:17
Most FFTs have an implicit scale factor in at least one direction (usually a scale factor of N in the forward direction) - I'm not familiar with your particular FFT routine but you should probably check which convention it uses as you may well need to scale your output values by 1/N. — Paul R, Apr 15 '15 at 22:11
How big is your `sound` array? It may not necessarily be the same problem that you're having, but I got a sudden jump in the maximum absolute error when comparing input & output for a pure sinusoidal tone as the input size increased from 46341 to 46342 (way more than the typical slowly increasing error with FFT size). The rms error remained under wraps though. — SleuthEye, Apr 16 '15 at 13:14
Start checking whether the entries of your complexInput vector before the Fourier.Forward and after the Forier.Inverse operations. If they're identical, then the Fourier part is ok. Otherwise the Fourier part is bad. — Jens, Apr 16 '15 at 19:33
@Jens I checked the max value before and after, and before is bigger (0,999..) than after (0,800..), after is with the noise. — Dejan Skledar, Apr 21 '15 at 15:49
@SleuthEye the sound array's length is about 800 000 samples. — Dejan Skledar, Apr 21 '15 at 15:51
Turns out that a similar problem had been reported on [gitub](https://github.com/mathnet/mathnet-numerics/issues/286) — SleuthEye, Apr 23 '15 at 01:25

SleuthEye · Accepted Answer · 2015-06-01T01:00:32.450

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 n²) 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.

This looks promising, but now I am stuck with the convolution part, what I am doing is looping through the complex array, and creating a new one using complexNew[i] = complexInput[i] * complexEffect[i]; and then again making the IFFT on the complexNew, but I am getting some "shattering" at some high amplitudes... — Dejan Skledar, Apr 22 '15 at 17:23
Multiplication in the frequency domain corresponds to _circular_ convolution, which would result in unexpected artifacts when comparing against linear convolution. To avoid this, you need to zero pad your time-domain signals to be convoluted to the sum of their lengths (i.e. twice the size if they are both the same length). — SleuthEye, Apr 22 '15 at 17:46
I have fixed my this error, and works fine now, thank you. I may have some more questions in the future, I hope to hear from you again. — Dejan Skledar, Apr 22 '15 at 17:47

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

1 Answers1