FFTW returning double ghost frequency content

Intro

I am just posting the comment as an answer so in case this is the issue you can accept it and have a solution for future reference.

Confusion

After taking a closer look at your code there seems to be some kind of inconsistency, most probably stemming from bad understanding of the workings of FFTW (either mine or yours), or just missing information. Your fftLength seems to be half the length of the input signal, which I assume you set it like that seeking to discard the negative frequencies. However, it is not clear whether the output variable/array holds the whole spectrum or the positive frequencies only. I assume this depends on the way FFTW calculates things (possibly depending on the options passed to the “plan”).

Below I’ll discuss the bug in your code related to the calculation of the magnitude of the frequency bins, but clarifying this issue (what does the output array actually holds/contains?) is crucial to get correct results (as it is related to the maximum value of your index).

Bug

There seems to be a bug in your code in the last for-loop in your ExamplePerformFFT() function. There is also an error at the comment that gives the length of the array (reads Length = 2 * (2/N - 1) while it should be Length = 2 * (N/2 - 1) but since this does not affect the results won’t be discussed further). I replicate here the for-loop for completeness

for (int bin = 0; bin < output.length - 1; bin++) { finalSpectrum[bin] = (float)Math.Sqrt(complexOutput[bin] * complexOutput[bin] + complexOutput[bin * 2 + 1] * complexOutput[bin * 2 + 1]) * scalingFactor; } 

The bug is in the way you index things. Your bin index goes from 0 to the length of the spectrum array (as mentioned above 2 * (N/2 - 1)) which here I will “bypass” the issue mentiond in “Confusion” above and assume that it goes through all the bins one by one. This means that the length of the complexOutput array must be twice the length of your finalSpectrum variable (since it is supposed to hold real and imaginary values separately). However, your actual bin complex values are intertwined with the real and imaginary parts being interleaved in the complexOutput variable. This means that the real parts will “reside” at indices 2 * bin, which are the even numbered/indexed elements of the complexOutput array and the imaginary parts at indices 2 * bin + 1.

You get the imaginary parts of the bin complex values correctly by indexing with 2 * bin + 1, but the index for the real part is incorrect since you index with just bin. I will demonstrate the first 3 iterations as a demo.

1. Iteration $1$ - bin = 0
   • $\texttt{finalSpectrum[0]} = \sqrt{\texttt{complexOutput[0]}^{2} + \texttt{complexOutput[1]}^{2}}$
2. Iteration $2$ - bin = 1
   • $\texttt{finalSpectrum[1]} = \sqrt{\texttt{complexOutput[1]}^{2} + \texttt{complexOutput[3]}^{2}} \rightarrow$ The indices for the complexOutput should be $2$ (real part) and $3$ (imaginary part) respectively, however you get $1$ (actually the imaginary part of the $1\textrm{st}$ bin complex value) and $3$ (the correct imaginary part of the $2\textrm{nd}$ bin).
3. Iteration $3$ - bin = 2
   • $\texttt{finalSpectrum[2]} = \sqrt{\texttt{complexOutput[2]}^{2} + \texttt{complexOutput[5]}^{2}} \rightarrow$ The indices for the complexOutput should be $4$ (real part) and $5$ (imaginary part) respectively, however you get $3$ (actually the real part of the $2 \textrm{nd}$ bin complex value) and $5$ (the correct imaginary part of the $3\textrm{nd}$ bin).

I hope you get the pattern, which definitely gives you wrong results.

Now, there is a couple of ways you could “solve” the issue. In fact all of them boil down to picking the correct indices, however you could either set the bin indices in such a way as to make the indexing of the complexOutput easier or the finalSpectrum straight forward. I believe they shouldn’t have any impact on performance, but I am not a software engineer (or a skilled DSP engineer either) so if this is not the case, someone hopefully will shed some light on it in the comments.

The two different ways to index are below (well there may be more but those two come to mind that are convenient one way or the other).

1. Index for the finalSpectrum. This entails that you keep the indexing as it is and you correct the indices for the real part, as suggested above. This would make your for-loop look like

for (int bin = 0; bin < output.length - 1; bin++) { finalSpectrum[bin] = (float)Math.Sqrt(complexOutput[2 * bin] * complexOutput[bin] + complexOutput[2 * bin + 1] * complexOutput[2 * bin + 1]) * scalingFactor; } 

Note that the only difference here compared to your original for-loop is the index for the real part of each bin becoming 2 * bin instead of bin. This is an easy tweak that you could test very easily and fast.

2. Alternatively, you could index in a way that will make the indexing of the real and imaginary parts of the bins easy to express. However you have to “compensate” for the index of the finalSpectrum. This would look like

for (int bin = 0; bin < complexOutput.length - 1; bin += 2) { finalSpectrum[bin/2] = (float)Math.Sqrt(complexOutput[bin] * complexOutput[bin] + complexOutput[bin + 1] * complexOutput[bin + 1]) * scalingFactor; } 

Here you change the range of your index to cover the length of the complexOutput array and you increment the index by 2 in each step, effectively moving from one bin to the next (the indices will point to the real parts of the bins). Now the real part you need is exactly at bin and the corresponding imaginary part at bin + 1; and these are the indices you use for the calculation of the magnitude of the bins. However, the index corresponding to the bin in the finalSpectrum is always half the bin value, since the latter increases by $2$ in each step. So the index for the finalSpectrum variable is bin/2. Since bin will always be even (zero inclusive) bin/2 will be an integer so you are good to go.

As you can see, both methods achieve the same results and it is a matter of preference which one you are going to use. However, I would urge you to double-check the actual lengths of the output, finalSpectrum and complexOutput arrays and make sure they are “compatible” and what you want them to be. Furthermore, do make sure that the output array does hold the part of the spectrum you want, or else (if it holds both positive and negative frequencies) you are iterating over the whole spectrum (and not only the right/left part of it). If this is intended behaviour that’s fine, but if not you should act accordingly.

I hope this helps and please don’t hesitate to post comments and ask for clarifications (or suggest corrections if I missed something).