Time-frequency analysis, Morlet wavelet

Hi Mike,

I’m following your book and trying to do Time-frequency analysis, here below is my script, I combined it with your many-trials code (lectures). I’ve tried it on my data, however, I didn’t get any result. I don’t know where I did it incorrectly. It says continue entering statement but I put all for’s end!

min_freq = 2;
max_freq = 40;
num_frex = 40;
frex = linspace(min_freq,max_freq,num_frex);

% which channel to plot
for channel2use = 1:EEG_im_IClabel.nbchan 

    % other wavelet parameters
    range_cycles = [ 4 10 ];

    s = logspace(log10(range_cycles(1)),log10(range_cycles(end)),num_frex) ./ (2*pi*frex);
    wavtime = -2:1/EEG_im_IClabel.srate:2;
    half_wave = (length(wavtime)-1)/2;

    % FFT parameters
    nWave = length(wavtime);
    nData = EEG_im_IClabel.pnts;
    nConv = nWave + nData - 1;

    % initialize output time-frequency data
    tf = zeros(length(frex),EEG_im_IClabel.pnts,EEG_im_IClabel.trials);
    
    % Baseline
    baselinetime = [-150 0]; % in ms 

    % convert baseline window time to indices
    [~,baselineidx(1)]=min(abs(EEG_im_IClabel.times-baselinetime(1)));
    [~,baselineidx(2)]=min(abs(EEG_im_IClabel.times-baselinetime(2)));


    % SNR
    snr_bs = zeros(length(frex),EEG_im_IClabel.pnts, EEG_im_IClabel.trials);
    snr_tf = zeros(length(frex),EEG_im_IClabel.pnts, EEG_im_IClabel.trials);

    for fi=1:length(frex)
    
    % create wavelet and get its FFT
    % the wavelet doesn't change on each trial...
        wavelet  = exp(2*1i*pi*frex(fi).*wavtime) .* exp(-wavtime.^2./(2*s(fi)^2));
        waveletX = fft(wavelet,nConv);
        waveletX = waveletX ./ max(waveletX);
    
    % now loop over trials...
         for triali=1:EEG_im_IClabel.trials
        
            dataX = fft(squeeze(EEG_im_IClabel.data(channel2use,:,triali)), nConv);

            % run convolution
            as = ifft(waveletX .* dataX);
            as = as(half_wave+1:end-half_wave);

            % put power data into time-frequency matrix
            tf(fi,:,triali) = abs(as).^2;

            % extract SNR in two ways
            snr_tf(fi,:,triali) = mean(abs(as).^2,2)./std(abs(as).^2,[],2);
            snr_bs(fi,:,triali) = mean(abs(as).^2,2)./std(mean(abs(as(baselineidx(1):baselineidx(2),:)).^2,1),[],2);

         end
     
    end 
end

tfTrialAve = mean(tf,3);
snr_tfAve = mean(snr_tf,3);
snr_bsAve = mean(snr_bs,3);

baseline_power = mean(tfTrialAve(:,baselineidx(1):baselineidx(2)),2);
baselinediv = tfTrialAve ./ repmat(baseline_power,1,EEG_im_IClabel.pnts);




% plot SNR
figure
subplot(121)
contourf(EEG_im_IClabel.times,frex,snr_bsAve,40,'linecolor','none')
set(gca,'clim',[.5 2],'xlim',[-150 1298],'ylim',[frex(1) 40])

% colorbar
title('SNR_b_a_s_e_l_i_n_e (mean/std)'), ylabel('Frequency (Hz)'), xlabel('Time (ms)')
axis square

subplot(122)
plot(snr_bsAve(1:3:end),baselinediv(1:3:end),'.')
xlabel('SNR_b_a_s_e_l_i_n_e'), ylabel('Power (/baseline)')
axis square


figure
subplot(121)
contourf(EEG_im_IClabel.times,frex,snr_tfAve,40,'linecolor','none')
set(gca,'clim',[.5 1.25],'xlim',[-150 1298],'ylim',[frex(1) 40])
% colorbar
title('SNR_t_f (mean/std)'), ylabel('Frequency (Hz)'), xlabel('Time (ms)')
axis square

subplot(122)
plot(snr_tfAve(1:3:end),baselinediv(1:3:end),'.')
xlabel('SNR_t_f'), ylabel('Power (/baseline)')
axis square

% time-series of SNR

figure
plot(EEG_im_IClabel.times,abs(mean(EEG_im_IClabel.data(8,:,:),3) ./ std(EEG_im_IClabel.data(8,:,:),[],3)))
title('Time-domain SNR time series')
set(gca,'xlim',[-150 1298])
xlabel('Time (ms)'), ylabel('SNR')

% now compute SNR of peak compared to prestim noise
stimeidx = dsearchn(EEG_im_IClabel.times',-150);
etimeidx = dsearchn(EEG_im_IClabel.times',400);
disp([ 'ERP SNR between 150 and 400 ms at FCz: ' num2str(max(mean(EEG_im_IClabel.data(8,stimeidx:etimeidx,:),3)) / std(mean(EEG_im_IClabel.data(8,baselineidx(1):baselineidx(2),:),3),[],2)) ])

Hi Lili. Variable as is a vector, but when you try to get the standard deviation to create snr_bs, you’re indexing a second dimension. If you template that after the previous line (snr_tf), the code will run.

Mike

Thank you so much Mike. I changed a little my code. I used your convolution over all trials code and I created one long vector. Now I want to compare the arms of the baseline period with the stimulus period of all trials. However, I’m confused with the baseline concept. If the time range of each epoch is [0 1148] and I put x.min to -150 to consider [-150 0] as my baseline and [0 998] stimulus period, should I do preprocessing on only [0 998] time period? Because I did all of the ICA, artifact rejection, the morletwavelet convolution on the total time. Since I cleaned the baseline period is it correct to compare the arms of the stimulus period with the baseline period as a signal-to-noise ratio?

This is my code if it is correct:
%baseline parameters

baselinetime = [ -150 0 ]; % in ms
stimulustime = [0 998]

% convert baseline window time to indices

[~,baselineidx(1)]=min(abs(EEG.times-baselinetime(1)));
[~,baselineidx(2)]=min(abs(EEG.times-baselinetime(2)));
[~,stimidx(1)]= min(abs(EEG.times-stimulustime(1)));
[~,stimidx(2)]= min(abs(EEG.times-stimulustime(2)));

after convolution over all trials and before reshaping the long vector of convolution result, I wrote this:
#rms
baseline(fi,:slight_smile: = sqrt(mean(abs(as(baselineidx(1):baselineidx(2)).^2 ,2)
stim(fi,:slight_smile: = sqrt(mean(abs(as(stimidx(1):stimidx(2)).^2 ,2)

   # averaging and squaring all rms 

    baselineAve = mean(abs(baseline(fi,:).^2,2);
    stimAve = mean(abs(stim(fi,:).^,2);

  snr(fi,:) = baselineAve/stimave

To have one value should I get an average from snr(fi,:)? Is this correct in your view?

Since I’m new to these areas, my questions may sound too easy. Thank you so much.

Yes, you want to have a clean baseline, so it’s correct to process and clean all of your data.

I think your RMS code is missing a parenthesis. It probably should be:
sqrt(mean(abs(as(baselineidx(1):baselineidx(2))).^2 ,2))

Otherwise, I think the code looks OK. You are computing SNR as base/stim; maybe that’s what you want but normally I would think it should be stim/base, assuming that the baseline is your estimate of “noise”.

Thanks! I did it and it worked however since I want one value for signal noise ratio for each subject and experimental condition. And with this code SNRs are calculated for each frex and for each channel should I get an average on both fi and channels dimensions?
Also, when I use the channels loop I don’t get any result and I see size errors. When we use channels loop where should we incorporate channel variable? because I want to have all results of loops to be able to get an average and since I didn’t use the channel variable I’ll only have the last channel result.

This is my code:

Power_signal = zeros(EEG.nbchan,length(frex),EEG.pnts);

for channels = 1:length(channel_labels) 
    channel2use = channel_labels{channels}
    
% other wavelet parameters
range_cycles = [ 4 10 ];

s = logspace(log10(range_cycles(1)),log10(range_cycles(end)),num_frex) ./ (2*pi*frex);
wavtime = -2:1/EEG_im_IClabel.srate:2;
half_wave = (length(wavtime)-1)/2;

%baseline parameters

baselinetime = [ -150 0 ]; % in ms
stimulustime = [0 998]

% convert baseline window time to indices

[~,baselineidx(1)]=min(abs(EEG.times-baselinetime(1)));
[~,baselineidx(2)]=min(abs(EEG.times-baselinetime(2)));
[~,stimidx(1)]= min(abs(EEG.times-stimulustime(1)));
[~,stimidx(2)]= min(abs(EEG.times-stimulustime(2)));

% FFT parameters
nWave = length(wavtime);
nData = EEG.pnts * EEG.trials; 
nConv = nWave + nData - 1;

% now compute the FFT of all trials concatenated
alldata = reshape( EEG.data(strcmpi(channel2use,{EEG.chanlocs.labels}),:,:) ,1,[]);
dataX   = fft( alldata ,nConv );


% loop over frequencies
for fi=1:length(frex)

    % create wavelet and get its FFT
    wavelet  = exp(2*1i*pi*frex(fi).*wavtime) .* exp(-wavtime.^2./(2*s(fi)^2));
    waveletX = fft(wavelet,nConv);
    waveletX = waveletX ./ max(waveletX);

    % run convolution 
    as = ifft(waveletX .* dataX);
    as = as(half_wave+1:end-half_wave);
    
    % and reshape back to time X trials
    as = reshape( as, EEG.pnts, EEG.trials );

    
    % compute power and average over trials
    Power_signal(channels,fi,:) = mean( abs(as).^2 ,2);
    
    % SNR
    baseline_power(channels,:,:) = mean(Power_signal(:,baselineidx(1):baselineidx(2)),3);
    stimulus_power(channels,:,:) = mean(Power_signal(:,stimidx(1):stimidx(2)),3);
    SNR(channels,:,:) = stimulus_power./baseline_power
    SNR_db(channels,:,:) = 10*log10(stimulus_power./baseline_power)
    
    
end


% plot results
figure(channels), clf
contourf(EEG.times,frex,Power_signal,40,'linecolor','none')
set(gca,'clim',[0 5],'ydir','normal','xlim',[-150 1000])
title(['Power averaged over trials for channel ', channel_labels{channels}]) 
xlabel('Time (ms)')
ylabel('Frequency')

end

SNRAvechan = mean(SNR(channels,fi,:),[1,2])
SNRAvechan_db = mean(SNR_db(channels,fi,:),[1,2])


Thank you so much!

You’re discovering the dimension-explosion issue in EEG analysis :wink:

You can compute SNR for every channel and frequency. Averaging over all channels and all frequencies is probably not the best approach, because most dynamics in EEG are mainly present in some electrodes and some frequencies.

So the best approach might be to pick one or a few channels, and a frequency band, and then compute SNR only on those data. Averaging together a few neighboring channels and a few neighboring frequencies is usually a good idea.