fix to PSP metric to account for unbiased cross-correlation

ibs-lab / cedalion

fNIRS analysis toolbox

MIT License

21 stars 8 forks source link

fix to PSP metric to account for unbiased cross-correlation #33

Closed lauracarlton closed 1 month ago

lauracarlton commented 2 months ago

Bug fixes in PSP

updated to match the qtnirs matlab implementation
unbiased scaling of the cross-correlation timeseries
calculate the power spectrum density instead of just power spectrum

emiddell commented 2 months ago

Hi @lauracarlton, one more idea how we could avoid the two nested loops over windows and channel. Welch's method improves the estimate of the spectral density by using a sliding window approach on the passed time series. This we don't use. nperseg and nfft are both set to the length of the whole time series (here the window length len(corr[ch,:]). So from signal.welch we only use the application of the Hamming window and the scaling of the density estimate. If we apply these two things ourselves, we could use numpy's implementation of the fft for real valued signals numpy.fft.rfft. This function operates on n-dimensional arrays and one can specify via the axis parameter along which array the fft should be calculated. Should speed things up.

lauracarlton commented 2 months ago

Thanks for that idea @emiddell ! I can try implementing it and see how it goes:) I realized yesterday I am still getting different answers than the qt-nirs toolbox is returning. Meryem found this forum on the mne-nirs page https://github.com/mne-tools/mne-nirs/issues/503 Our implementation is a little different than how mne does it but I was using the same method for calculating PSD as they are so I think I am running into a similar issue documented there. Maybe your implementation suggestion above could help address this.

lauracarlton commented 2 months ago

I have attempted an implementation using the method you described above. I am still getting different results compared to the qt-nirs matlab implementation ... Meryem suggested we reach out to Luca Pollonini to see if he has any ideas. Let me know if something jumps out to you though.

emiddell commented 2 months ago

In the linked thread it was pointed out that the difference between qt-nirs and mne-nirs could stem from the former calculating the power spectrum and the latter calculating the power spectral density

qt-nirs calculates the power spectrum by using spectrumtype='power'
mne-nirs calculates the power spectral density by using scipy.signal.periodogram with default parameters.

Scipy's different scaling constants to get the PSD or the power spectrum are defined here.

Your current implementation seems to calculate the PSD, too. Hence, I would expect that you get other result than qt-nirs.

emiddell commented 1 month ago

Changes were merged in e802f72.