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)) ])
```