usnistgov / microscope

Oscilloscope-style plotting for microcalorimeter data
4 stars 2 forks source link

PSD normalization #53

Open NanoExplorer opened 2 days ago

NanoExplorer commented 2 days ago

I don't think the PSD is being normalized properly. I'm going to start working on this, but making an issue to document.

In FFTComputer.cpp, https://github.com/usnistgov/microscope/blob/47b248a12ad5c5874fbb7ada44cbe8053f6a1729/fftcomputer.cpp#L155 the sampleRate parameter is used to convert the PSD from volts/sqrt(sample) to volts/sqrt(Hz). However here: https://github.com/usnistgov/microscope/blob/47b248a12ad5c5874fbb7ada44cbe8053f6a1729/pulsehistory.cpp#L108-L109 the sample rate is hard-coded to 1, and also here: https://github.com/usnistgov/microscope/blob/47b248a12ad5c5874fbb7ada44cbe8053f6a1729/pulsehistory.cpp#L236

I'll start seeing if I can get the normalization right and then add a pull request when I think I've figured it out.

joefowler commented 2 days ago

@NanoExplorer : I don't know why the sample rate is being fixed to 1.0 in both relevant calls to the computePSD method, but it makes no sense. Maybe it was a placeholder from a time when I had no convenient access to the true sample rate (or period)? Heck, even now, do we have access to the rate or the period at the lines you've highlighted? Is it a matter of transporting that information into a class that does not currently have it?

NanoExplorer commented 2 days ago

It looks like that information is right there. I'm testing a version where I've edited the lines in question to: fftMaster->computePSD(records[i]->data, *psd, 1/records[i]->sampletime, WINDOW, previous_mean); Since the pulseRecord class (in that line records is an array of pulseRecord pointers) has a member variable for the sample time. At least with the demo python scripts, the amplitude of the PSD does change when I change the sample time in the fake pulses, but I'm going to think on it a bit more to make sure it's giving reasonable values.

NanoExplorer commented 2 days ago

OK, I may have dug up another problem? If you look at the following lines: https://github.com/usnistgov/microscope/blob/47b248a12ad5c5874fbb7ada44cbe8053f6a1729/fftcomputer.cpp#L89-L101 dividing the window function by the mean squared doesn't actually do what the comment wants, which is to have average(window[i]^2)==1. In this case the average of the "normalized" window squared is about 2.7. I'm pretty sure that we should be dividing the window by its RMS., i.e. replace line 100 with window[i] /= sqrt(meansq); (though that's not the actual code I'd use)

I think I might still be off by a factor of 2 though, so I'll have to figure out if that's actually the case.

NanoExplorer commented 1 day ago

I did a more thorough test and I am now getting consistent results from Microscope and scipy.signal.periodogram, as I mention in the pull request.