Time-frequency analysis

Hi Mike,

I’m following your book and course in Udemy and trying to do Time-frequency analysis, here below is my script, I’ve tried it on my data already, the result looks normal, but I’m still not sure if this is correct? could you help me take a look? thanks in advance!

% wavelet parameters
num_frex = 50; % define the number of frequency
min_freq = 3;
max_freq = 30;

channel2use = ‘FCz’;
baseline_window = [ -300 -100];
timewindow = [200 500 ];
Timeidx = dsearchn(EEG.times’,timewindow’);

% set a few different wavelet widths (number of wavelet cycles)
range_cycles = [ 4 8 ];

% other wavelet parameters
nCycles = logspace(log10(range_cycles(1)),log10(range_cycles(end)),num_frex); % from log10(4) to log10(10), there are 40 points (based on the number of frequency)
frex = linspace(min_freq,max_freq,num_frex); % from 2 to 30 Hz, including 40 points
time = -2:1/EEG1.srate:2; % the length of wavelet
half_wave = (length(time)-1)/2;

num_sub = [1:10]; % the number of subjects
% initialize output time-frequency data
tf = zeros(length(frex),EEG1.pnts);
% frex * timepoints
% according to personal analysis goal: Theta band power, 200-400ms
% Theta band power: 4-8Hz: from 4 to 10 rows
% Time window (200-400ms): [205, 308]
tf_mean_ForEachSub = cell(length(num_sub),1);

% load data for each subject
for ssnr = [1:10,12,14,17:18,20:26,28:34];
if ssnr<10;
m=[‘0’,num2str(ssnr)];
else
m=num2str(ssnr);
end

EEG = pop_loadset('filename',strcat('Induced_subject', m,'_posFB.set'),'filepath', ACPath);   %% load data for each

% FFT parameters
nKern = length(time); % length of wavelet
nData = EEG.pnts*EEG.trials; % length of data (timepoints * trials)
nConv = nKern+nData-1; % the length of results

% FFT of data (doesn’t change on frequency iteration)
dataX = fft(reshape(EEG.data(strcmpi(channel2use,{EEG.chanlocs.labels}),:,:),1,[]),nConv);

% convert baseline time into indices
[~,baseidx(1)] = min(abs(EEG.times-baseline_window(1)));
[~,baseidx(2)] = min(abs(EEG.times-baseline_window(2)));

for fi=1:length(frex)

% create wavelet and get its FFT
s = nCycles(fi)/(2*pi*frex(fi));  % number of cycles
cmw = exp(2*1i*pi*frex(fi).*time) .* exp(-time.^2./(2*s^2));  % create complex wavelet
cmwX = fft(cmw,nConv);  % get its FFT

% run convolution (as is the result of convolution, which contain 3
% indexes: power, phase, filtered signal)
as = ifft(cmwX.*dataX,nConv);
as = as(half_wave+1:end-half_wave);
as = reshape(as,EEG.pnts,EEG.trials);  % this is the data where contain the result of power

% put power data into big matrix
tf(fi,:) = mean(abs(as).^2,2);  % the formula for extracting the power 
% this is the value refers to the distance from the point to the
% origin

end
tf = 10*log10( bsxfun(@rdivide, tf, mean(tf(:,baseidx(1):baseidx(2)),2)) ); % for each subject, updated
tf_index_data = tf([4:8], [Timeidx(1):Timeidx(2)]); % to get the data for 4-8Hz, 280-320ms
tf_mean = mean(tf_index_data(:)); % extract the mean value
tf_mean_ForEachSub {ssnr,1} = tf_mean; % the TF power is sotred in the file ‘tf_mean’, for each subject

end

Hi echo. From a quick glance it looks mostly fine. Two ways to be more confident about it would be to (1) simulate data and run it through your code to make sure the resulting plot accurately reflects what you simulated; (2) use your code to analyze the book sample data to see if you can reproduce one of the figures in the book.

One potential mistake I see is at the end of the code where you extract power in a time-frequency window. I’m not sure that 4:8 in frequency indices corresponds to 4-8 Hz. It would be good to soft-code that.

Hi Mike,

Thanks a lot for your reply and suggestions, I’ll follow your suggestion to re-check.

Best,
Echo

Hi Mike,
As you supposed, the “4:8” in frequency indices did not correspond to 4-8Hz, but instead 4,6531-6.8571Hz. In this case, I wonder how do I get the accurate indices for the 4-8Hz? according to the frequency parameters I defined, there was no integer 4 and 8, the proximal number was 4.1020 and 7.9592. Now I come up with the code as follows to define the frequency indices:
FrequencyWindow = [4 8];
Fqidx = (dsearchn(frex’,FrequencyWindow’))’;
Fqidx = [Fqidx(1):Fqidx(2)]; % the index of frequency

By implementing this code, it can only give me the value : 4.1020 - 7.9592 (Hz), is this ok?

Thanks in advance for your reply.

Yup, that looks good. With the amount of smoothing imposed by wavelet convolution (or any other time-frequency analysis method), you don’t need to take the frequency boundaries so precisely. “4-8 Hz” is the same thing as “4.1020 - 7.9592 Hz.”