One extra array element from Magnitude and Power methods?

I am passing 2048 PCM samples to either FFT.Magnitude or FFT.Power, and I'm getting back 1025 elements instead of 1024.

This isn't a big deal since, of course, most of the "far end" (for music) is zero or near-zero, but I wasn't expecting that. Is this known/expected for some reason? Or am I doing something wrong? Nothing fancy in the code, reading 16-bit short samples via OpenAL:

private void FrequencyProcessing()
{
    // FFT buffer is 2048, "slide" two passes of PCM 1024 data through it.
    // Copy the second half of the FFT buffer "back" to the first half:
    Array.Copy(BufferFFT, BufferLength, BufferFFT, 0, BufferLength);
    // Next, copy new PCM data to the second half of the FFT buffer:
    Array.Copy(BufferPCM, 0, BufferFFT, BufferLength, BufferLength);

    var window = new Hanning();
    double[] windowed = window.Apply(BufferFFT);
    double[] zeroPadded = Pad.ZeroPad(windowed);
    System.Numerics.Complex[] spectrum = FFT.Forward(zeroPadded);

    // this returns 1025 elements?
    BufferFreq = FFT.Magnitude(spectrum);
}

Hi @MV10,

This extra point is only present when the positiveOnly argument is true (the default case).

I think the lowest number represents the DC component, and the highest number is the magnitude of the signal at the Nyquist frequency (sample rate / 2). When the full array is returned the data is mirrored around these points, and the length is equal to the input length.

FftSharp/src/FftSharp/FFT.cs

Lines 157 to 174 in 6f341bb

    
           /// <summary> 
        
           /// Calculate power spectrum density (PSD) in RMS units 
        
           /// </summary> 
        
           public static double[] Magnitude(System.Numerics.Complex[] spectrum, bool positiveOnly = true) 
        
           { 
        
               int length = positiveOnly ? spectrum.Length / 2 + 1 : spectrum.Length; 
        
               double[] output = new double[length]; 
        
               // first point (DC component) is not doubled 
        
               output[0] = spectrum[0].Magnitude / spectrum.Length; 
        
               // subsequent points are doubled to account for combined positive and negative frequencies 
        
               for (int i = 1; i < output.Length; i++) 
        
                   output[i] = 2 * spectrum[i].Magnitude / spectrum.Length; 
        
               return output; 
        
           }

However, if convention in other software packages does something different, I'd be happy to modify the functionality of FftSharp to return data more consistent with what others do

FYI I'm looking over https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf page 2 "Converting from a Two-Sided Power Spectrum to a Single-Sided Power Spectrum" to see if I can get some insights here

I'm adding resources here as I find them...

Here's some matlab that shows the N/2+1 behavior FftSharp currently uses

https://www.mathworks.com/matlabcentral/answers/225440-how-to-extract-phase-and-amplitude-information-from-fft

load('fb2040'); % loading the data
x = fb2040(3,:);
y_filt = filter(b,a,x); % filtering the received signal
nfft = length(y_filt);
res = fft(y_filt,nfft)/ nfft; % normalizing the fft
f = fs/2*linspace(0,1,nfft/2+1); % choosing correct frequency axes
res = res(1:nfft/2+1);
figure, plot(f,abs(res));
xlabel('Frequency in MHz');ylabel('Amplitude');
return

also... omg internet matlab code is so hard to read 💀

Okay, I think I figured this out my original reasoning. My understanding was that the DC component (gray square) and Nyquist frequency (yellow square) are shared by both the positive (blue) and negative (orange) magnitudes.

https://www.gaussianwaves.com/2015/11/interpreting-fft-results-complex-dft-frequency-bins-and-fftshift/

... however, is this the best behavior? I'm still deciding 😅

I confirmed the length N/2+1 behavior is shared by Python (SciPy), so I'm satisfied with that! I think thinks are good as they are 😎

Thanks for raising the original question @MV10! This topic isn't very intuitive, and I'm glad to have some written justification now about why things are the way they are.

import numpy as np
import scipy.fft

xs = np.arange(1024) / 10
data = np.sin(xs)
print(f"data length: {len(data)}")

fft = scipy.fft.fft(data)
print(f"fft length: {len(fft)}")

rfft = scipy.fft.rfft(data)
print(f"rfft length: {len(rfft)}")

data length: 1024
fft length: 1024
rfft length: 513

Gotcha. Interesting. Thank you, this library is incredibly handy!

	/// <summary>
	/// Calculate power spectrum density (PSD) in RMS units
	/// </summary>
	public static double[] Magnitude(System.Numerics.Complex[] spectrum, bool positiveOnly = true)
	{
	int length = positiveOnly ? spectrum.Length / 2 + 1 : spectrum.Length;

	double[] output = new double[length];

	// first point (DC component) is not doubled
	output[0] = spectrum[0].Magnitude / spectrum.Length;

	// subsequent points are doubled to account for combined positive and negative frequencies
	for (int i = 1; i < output.Length; i++)
	output[i] = 2 * spectrum[i].Magnitude / spectrum.Length;

	return output;
	}