Power method does not compute true PSD

[jwilo]

@swharden,

Love your work across both this and other repos of yours such as ScottPlot, have used both extensively.

I've been using FftSharp recently, and stumbled across what may potentially just be a documentation issue, rather than a functional one.

The Power() method's description claims to return a Power Spectrum Density (PSD), which I believe is possibly an interpretation/blending of two concepts; a Power Spectrum, and a Power Spectral Density (PSD). The implementation appears to simply compute the square of an FFT, thus computing a Power Spectrum, not a Power Spectral Density (I appreciate it claims to produce a Power Spectrum Density, but the density portion of this claim is not handled, and the term PSD usually refers to Power Spectral Density).

FftSharp/src/FftSharp/FFT.cs

Lines 242 to 253 in 3f5158f

	/// <summary>
	/// Calculate power spectrum density (PSD) in dB units
	/// </summary>
	public static double[] Power(System.Numerics.Complex[] spectrum, bool positiveOnly = true)
	{
	double[] output = Magnitude(spectrum, positiveOnly);
	for (int i = 0; i < output.Length; i++)
	output[i] = 20 * Math.Log10(output[i]);
	return output;
	}

This Power Spectrum has correct Y axis units dB (as documented), but is not a Power Spectral Density, due to the lack of normalisation. A true PSD represents power per unit frequency, i.e. g^2√Hz (or V, or whatever others units you're using), as opposed to simple power.

Fully appreciate I might just be missing something elsewhere in the source somehow, so wanted to get your take on it?

I need the capability to easily compute PSDs in .NET applications, so would be happy to work on this for merging into FftSharp.

Furthermore, PSDs are normally calculated using Bartlett's method such that the PSD's DoF may be chosen according to the application;

• the time history is broken into smaller segments
• each segment is windowed
• the FFT computed for each segment (with variations on this such as Welch's method, to include overlap between the segments)
• each segment is squared
• the resulting squared FFTs are averaged

Thus, I don't propose 'fixing' Power(), as Power() isn't 'broken', the documentation for what it does is just misleading, IMO. I think this should be rewritten:

Calculate power spectrum density (PSD) in dB units

...replaced with

Calculate Power Spectrum in dB units

I then propose adding a PSD() method, that takes appropriate arguments to set the DoF, overlap % and window type used. Signature would likely take the form of (if computing directly from time history):

public static double[] Psd(double[] samples, int dof, int overlap, FftSharp.WindowT windowType)

Not 100% sure on the return type, I've only thought about this for a few minutes, not entirely sure what we'd need to return. Perhaps a PSD object containing the data + frequency bins... need to think that through some more.

Would love to hear your thoughts?