lpollonini / phoebe

PHOEBE (Placing Headgear Optodes Efficiently Before Experiment) is a MATLAB GUI application that measures and displays the optical coupling between fNIRS optodes and the scalp of a subject in real time.
8 stars 1 forks source link

peak power metric as a function of time window length #2

Open BAK262 opened 2 months ago

BAK262 commented 2 months ago

Hi,

Thank you very much for sharing this useful work!

I'm recently exploring your peak power metric on my own fNIRS data. There are two weird facts that I don't understand, both about the setting of the time window length.

  1. It looks like the peak power metric likewise oscillates at 2x the cardiac frequency as the length of the time window changes. Here is an example:

    fs = 10; % sampling frequency
    d = [3*sin(2*pi*((1/fs):(1/fs):10))', 0.1*sin(2*pi*((1/fs):(1/fs):10))'+1]; % simulated 10s filtered signals from same channel
    t = 1:(1/fs):10; % different time window length
    peakPower = nan(length(t),1);
    for tt = 1:length(t)
    x = d(1:(1+t(tt)*fs-1), 1);
    y = d(1:(1+t(tt)*fs-1), 2);
    xy = xcorr(zscore(x),zscore(y),'unbiased');
    pxy = periodogram(xy,hamming(length(xy)),length(xy),fs,'power');
    peakPower(tt) = max(pxy);
    end
    figure;plot((1./fs):(1./fs):t(tt),x);hold on;plot((1./fs):(1./fs):t(tt),y);hold off; % plot of the simulated filtered signals
    figure;plot(t,peakPower) % plot of peak power ~ time window length
  2. Both in the above example and in your publications, the filtered signals are assumed to have unchanged amplitude over time. But real cases are not usually like this assumption -- their amplitudes keep changing based on my observation. And then, the peak power metric may dramatically decrease if the signal amplitudes change within the time window. Here is an example:

    fs = 10; % sampling frequency
    d = [3*sin(2*pi*((1/fs):(1/fs):10))', 0.1*sin(2*pi*((1/fs):(1/fs):10))'+1]; % simulated 10s filtered signals from same channel
    multAmp = (1.005:0.005:1.5)'; % multiplier to the filtered signal
    d = d.*multAmp; % so that the signal amplitudes increase over time
    t = 1:(1/fs):10; % different time window length
    peakPower = nan(length(t),1);
    for tt = 1:length(t)
    x = d(1:(1+t(tt)*fs-1), 1);
    y = d(1:(1+t(tt)*fs-1), 2);
    xy = xcorr(zscore(x),zscore(y),'unbiased');
    pxy = periodogram(xy,hamming(length(xy)),length(xy),fs,'power');
    peakPower(tt) = max(pxy);
    end
    figure;plot((1./fs):(1./fs):t(tt),x);hold on;plot((1./fs):(1./fs):t(tt),y);hold off; % plot of the simulated filtered signals
    figure;plot(t,peakPower) % plot of peak power ~ time window length

As the 2 examples show, the peak power metric is neither 'independent of signal length' nor 'independent of amplitude' in realistic situations. It depends on the appropriate choice of time window length (signal length to compute it) -- which in turn depends on the cardiac frequency (example#1) and the time scale on which the signal amplitude is constant (example#2).

I would like to know if you and your colleagues have any further explanations for these two phenomena. From a practical view, what are the recommendations for setting the length of the time window?

Thanks in advance and best, Ming

Ph.D. candidate Department of Psychological and Cognitive Sciences Tsinghua University

Email: liming16@tsinghua.org.cn

BAK262 commented 2 months ago

I'm very sorry for the error in my #2 example! After a more careful comparison, I realized that it wasn't the change in amplitude, but the drift of the filtered signal that caused the drop in the peak power metric.

Now I have a new question: in the data preprocessing process, we usually filter the continuously recorded signal as a whole. Filtering ensures that the complete signal is de-trended, but the signal is still likely to be drifting within the finite time window used to calculate the peak power metric. Should I de-trend this time window segment again before calculating the peak power metric?

lpollonini commented 1 month ago

Hello Ming,

I am tackling your final question first: you can surely remove a low-frequency drift again after chopping the whole recording into small windows, although I doubt that any residual slow drift on such a short window would have any influence on the metrics.

Regarding the length of the window, I'd rather keep it as short as possible, yet including at least 2-3 cardiac periods. Practically, a 3-second window is what I use the most, since it is a good tradeoff between time resolution (and detection of motion artifacts) and estimation of measures.

I hope it helps, but please let me know if you have additional questions.