2

The following C# NAudio code produces a different result to MATLAB by a factor of 4. Why does this occur and is one of them incorrect? 

Complex[] tmp = new Complex[4];
tmp[0].X = 1.0f;
tmp[1].X = 0.5f;
tmp[2].X = 1.0f;
tmp[3].X = 0.25f;
tmp[0].Y = 0.0f;
tmp[1].Y = 0.0f;
tmp[2].Y = 0.0f;
tmp[3].Y = 0.0f;
FastFourierTransform.FFT(true, 2, tmp);

NAUDIO OUTPUT:

0.6875 + 0.0000i
0.0000 - 0.0625i
0.3125 + 0.0000i
0.0000 + 0.0625i

MATLAB OUTPUT:

2.7500 + 0.0000i
0.0000 - 0.2500i
1.2500 + 0.0000i
0.0000 + 0.2500i
Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
AlexS
  • 510
  • 2
  • 7
  • 23
  • I trust matlab implementation. To find the problem maybe you could use a longer signal with a specific frequency and then observe the spectrum. – Filipe Pinto Jun 11 '18 at 21:46

1 Answers1

5

The Discrete Fourier transform and its inverse require a certain normalization so that ifft(fft(x))==x. How this normalization is done changes from implementation to implementation.

It seems that, in this case, NAudio has chosen a different normalization than MATLAB.

MATLAB uses the most common normalization, where fft(x) at k=0 is equal to sum(x), and the inverse transform does the same thing but divides by n (number of samples). This is also the equation as described in the Wikipedia page for the DFT. In this case, the inverse transform matches the equation for the Fourier series.

NAudio seems to do the division by n in the forward transform, such that at k=0 you have mean(x).

Given the above, you can use the first frequency bin (the DC component) to verify what normalization is used (assuming there is a DC component, if the signal has a zero mean this will not work): If the DC component is equal to the sum of all the sample values, then the "common" normalization is used. It can also be equal to the sum divided by sqrt(n), in the case of a symmetric definition, where the forward and inverse transform carry the same normalization. In the case of NAudio it will be equal to the sum divided by n (i.e. the mean of the sample values). In general, take the DC component and divide it by the sum of the sample values. The result q is the normalization term used. The inverse transform should have a normalization term 1/qn.

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
  • I assume this means neither are incorrect and that it would always be unwise to compare the absolute energy/power of am magnitude plot with that computed by 2 different FFT algorithms unless you are sure some kind of level calibration has taken place. Meaning i can continue to use NAudio without the need to make any manual corrections to the complex data. – AlexS Jun 12 '18 at 17:42
  • @AlexS: The FFT doesn't compute energy or power, it computes the DFT. From there, you can easily compute energy/power, as long as you know the definition of the DFT that your FFT algorithm is using. Note that simply looking at the output for frequency 0 you can figure out what normalization is being used. -- What I mean to say is that, before you add units to any plot, you need to know what you have computed... :) – Cris Luengo Jun 12 '18 at 17:45
  • thanks, in that case can you please elaborate in your answer on how to confirm the normalisation method using the first frequency bin 0 Hz (or DC), I think this will then prove they are the same and the answer will be complete. – AlexS Jun 12 '18 at 19:44
  • @Alex: I have added a paragraph on how to examine the normalization. Is this what you had in mind? – Cris Luengo Jun 12 '18 at 21:21