There's a discrepancy between Bst and CTF power spectra units:
Looks like it could be a psd normalization factor of 2*pi or something like that. I know there are different normalization conventions when doing Fourier transforms, but probably only one of them gives the units T/sqrt(Hz) with "real frequencies" in the denominator and not angular frequencies. Something we should verify.
Great catch, Marc.
One way to test the consistency of the normalization of the PSD is described under Eq 12 of the following white paper:
"The area under the power spectrum is the variance of the process. So, a straightforward way to test normalization is to compute the PSD for a realization of X(t) with known variance, and zero mean [e.g. N (0, Ļ2 )]; and then calculate the integrated spectrum. For example, the single-sided PSD for a realization of a N (0, 1) process, sampled at 1 Hz, will be flat at 2 units2 Hzā1 across the frequency band [0, 1/2 ], and will have an area equal to one."
Would mind looking into this, please?
Thanks!
Thanks, I'm actually working on this right now. As I suspected the current method gives "Unit^2 / bin" (some call this U^2/(Hz.s) ) and not "U^2 / Hz" as is displayed in the figures. These spectrum amplitudes actually change depending on the window length.
There's also an issue with using the Hamming window but normalizing by N (which would be correct for a rectangular window). I'll do a PR, but I'm not sure how we'll want to implement this change as it wouldn't be comparable to previously computed PSDs.
I also haven't looked if any of this applies to TF.
Thanks Marc!
For legacy purposes, we probably need to have option checkboxes re: PSD methods wrt pre/post normalization fixes -- just like we do with Compute Source processes.
Done!
Here are figures for Gaussian noise, N(0,1) @ 1kHz, for the 3 new unit options:
physical
As expected, the amplitude (0.002 U^2/Hz) times bandwidth (500 Hz) = 1 U^2.
normalized frequency
Again, 2 * 0.5 = 1.
This "unit" of Hz.s is not commonly used, but I found this in a document that explained Welch's method in detail. I don't think many people will use this anyway - we might even want to hide this option...
"old"
The "integral" here, using bin numbers and not Hz, gives 0.79. Indeed there were a few issues that affected the global scaling, so the unit of U^2/"bin" is not exactly correct, but it's pretty close and it's understandable I think.
Great - thanks Marc!
So in the end, the recommendation is now to use the "physical" units I suppose?
Yes!
