Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

One extra array element from Magnitude and Power methods? #68

Closed
MV10 opened this issue Jun 29, 2023 · 5 comments
Closed

One extra array element from Magnitude and Power methods? #68

MV10 opened this issue Jun 29, 2023 · 5 comments

Comments

@MV10
Copy link

MV10 commented Jun 29, 2023

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);
}
@swharden
Copy link
Owner

swharden commented Jul 8, 2023

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.

/// <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

@swharden

This comment was marked as outdated.

@swharden
Copy link
Owner

swharden commented Jul 8, 2023

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.

image

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 😅

@swharden
Copy link
Owner

swharden commented Jul 8, 2023

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

@swharden swharden closed this as completed Jul 8, 2023
@swharden swharden mentioned this issue Jul 8, 2023
@MV10
Copy link
Author

MV10 commented Jul 8, 2023

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

swharden added a commit that referenced this issue Aug 21, 2023
resolves #69 and extends  #68
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants