fix to PSP metric to account for unbiased cross-correlation #33
Conversation
| @@ -106,6 +106,7 @@ def psp( | |||
| # Vectorized signal extraction and correlation | |||
| sig = windows.transpose("channel", "time", "wavelength", "window").values | |||
| psp = np.zeros((sig.shape[0], sig.shape[1])) | |||
| lags = np.arange(-nsamples + 1, nsamples) | |||
|
|
|||
| for w in range(sig.shape[1]): # loop over windows | |||
There was a problem hiding this comment.
use sig.sizes['time'] for clarity
src/cedalion/sigproc/quality.py
Outdated
|
|
||
| corr /= np.tile(norm_factor, (corr.shape[1],1)).T | ||
| # corr /= np.tile(norm_factor, (corr.shape[1],1)).T | ||
|
|
||
| for ch in range(sig.shape[0]): |
There was a problem hiding this comment.
same: sig.sizes['channel']
src/cedalion/sigproc/quality.py
Outdated
|
|
||
| corr /= np.tile(norm_factor, (corr.shape[1],1)).T | ||
| # corr /= np.tile(norm_factor, (corr.shape[1],1)).T | ||
|
|
||
| for ch in range(sig.shape[0]): | ||
| window = signal.windows.hamming(len(corr[ch,:])) |
There was a problem hiding this comment.
I think corr[ch,:] has the same length for all channels and windows. So this does not need to be recomputed and could be moved out of the two loops.
|
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. |
|
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 |
|
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. |
|
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
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. |
lauracarlton
force-pushed
the
PSP-fix
branch
from
6e9dd97
to
c8069fa
Compare
|
Changes were merged in e802f72. |
Bug fixes in PSP