Contents
Complex Morlet wavelet
% *대한뇌파연구회 집담회: 뇌파분석워크샵*
% This tutorial is published online at http://enslab.pbworks.com
% made by Ki-Young Jung, MD
% ENSLab. Neuroscience Research Institute
% Seoul National University College of Medicine
% code made based on 'Analysis neural time series data' by Mike X Cohen.
% This tutorial is made using a ERP data published in following paper.
% Cha KS, Sunwoo J-S, Byun J-I, et al.
% Working memory deficits in patients with idiopathic restless legs syndrome
% are associated with abnormal theta-band neural synchrony.
% Journal of Sleep Research. 2021;30:e13287.
% last updated May. 28 2022
understanding dot product and convolution in matlab
% make matrix
A = [1, 2, 3; 4, 5, 6] % 2 by 3 matrix
B = [7,8; 9, 0; 1, 2] % 3 by 2 matrix
% dot product between A and B matrix
AB_dot = A * B % result is 2 by 2 matrix
% convolution between A row vector and B col vector
AB_conv = conv(A(1,:), B(:, 1))
% length of conv result = length A + length B - 1
A =
1 2 3
4 5 6
B =
7 8
9 0
1 2
AB_dot =
28 14
79 44
AB_conv =
7
23
40
29
3
Part 1: analytic signal: single frequency
step 1. make complex wavelet kernel
clear all; close all;
srate = 200;
wtime = -2:1/srate:2;
hz = linspace(0,srate/2,floor(length(wtime)/2)+1);
freq = 5; % Hz center frequency to convolve %%%
s = 6/(2*pi*freq); % wavelet width
csine = exp(2*1i*pi*freq*wtime); % complex sine wave
gaus = exp( -(wtime.^2) / (2*s^2) ); % gaussian kernel
cmw = csine .* gaus; % complex morlet wavelet
cmwX = fft(cmw)/length(wtime); % fft of cmw
cmw_pow = cmwX .* conj(cmwX); % power of cmw
% plot of sine wave, gaussian kernel, cmw
figure(1), clf
subplot(131), plot(wtime, csine)
title('sine wave')
subplot(132), plot(wtime, gaus)
title('gaussian kernel')
subplot(133),
plot(wtime, real(cmw), 'b'),
hold on
plot(wtime, imag(cmw), 'b:'),
title('complex morlet wavelet')
% plot wavelet in 3D
figure(2), clf
subplot(121)
plot3(wtime,real(cmw),imag(cmw))
xlabel('Time'), ylabel('Real part'), zlabel('Imaginary part')
axis image
title('Complex Morlet wavelet in the time domain')
rotate3d on
subplot(122)
plot(hz, cmw_pow(1:length(hz)))
% plot(hz,abs(cmwX(1:length(hz))).^2)
set(gca,'xlim',[0 freq*3])
xlabel('Frequencies (Hz)'), ylabel('Amplitude')
title('Power spectrum of Morlet wavelet')
step 2: make fake EEG data to be analyzed
srate = 200;
datatime = -5:1/srate:5;
tim1 = [ones(1,200*2) zeros(1, 200*3) ones(1,200*2) zeros(1, 200*3) 1];
tim2 = [zeros(1,200*2) ones(1, 200*3) zeros(1,200*2) ones(1, 200*3) 1];
sine1 = exp(1i*2*pi*5.*(datatime.*tim1)); % 5 Hz
sine2 = exp(1i*2*pi*12.*(datatime.*tim2)); % 12 Hz
data = sine1 + sine2;
step 3: convolution between cmw and EEG data
% set parameter
nData = length(data);
nKern = length(wtime);
nConv = nData+nKern-1;
nHfkn = floor(nKern/2+1);
% convolution in frequency domain
fft_cmw = fft(cmw, nConv);
fft_data = fft(data,nConv);
wave_conv = ifft(fft_cmw .* fft_data);
wave_conv = wave_conv(nHfkn:end-nHfkn+1);
use matlab conv function in time domain
% wave_conv = conv(cmw, data);
% wave_conv = wave_conv(nHfkn:end-nHfkn+1);
step 4: get wavelet power
wave_pow = abs(wave_conv) .^2;
wave_pow = wave_pow/max(wave_pow);
% plot the time-domain signal
figure(3), clf
subplot(211)
plot(datatime, real(data))
xlabel('Time (s)'), ylabel('Amplitude')
title('Time-domain signal')
% the power time series at one frequency
subplot(212)
plot(datatime, wave_pow)
set(gca, 'YLim', [-1 2])
xlabel('Time (s)'), ylabel('Amplitude')
title([ 'Time-frequency power (' num2str(freq) ' Hz)' ])
% plot multiple frequency in time domain
figure(4), clf
for freq = 1:15
s = 5 /(2*pi*freq); % wavelet width
csine = exp(2*1i*pi*freq*wtime); % complex sine wave
gaus = exp( -(wtime.^2) / (2*s^2) ); % gaussian kernel
cmw = csine .* gaus; % complex morlet wavelet
fft_cmw = fft(cmw, nConv);
fft_data = fft(data,nConv);
wave_conv = ifft(fft_cmw .* fft_data);
wave_conv = wave_conv(nHfkn:end-nHfkn+1);
wave_pow = abs(wave_conv) .^2;
wave_pow = wave_pow/max(wave_pow);
subplot(3, 5, freq)
plot(datatime, wave_pow)
set(gca, 'YLim', [-1 2])
xlabel('Time (s)'), ylabel('Amplitude')
title([ 'Time-frequency power (' num2str(freq) ' Hz)' ])
end
part 2: analytic signal: multiple frequency
% set parameter
srate = 200;
wtime = -2:1/srate:2;
% make EEG time series
datatime = -5:1/srate:5;
tim1 = [ones(1,200*2) zeros(1, 200*3) ones(1,200*2) zeros(1, 200*3) 1];
tim2 = [zeros(1,200*2) ones(1, 200*3) zeros(1,200*2) ones(1, 200*3) 1];
sine1 = exp(1i*2*pi*5.*(datatime.*tim1)); % 5 Hz
sine2 = exp(1i*2*pi*12.*(datatime.*tim2)); % 12 Hz
data = sine1 + sine2;
% set convolution parameter
nData = length(data);
nKern = length(wtime);
nConv = nData+nKern-1;
nHfkn = floor(nKern/2+1);
% fft of signal
nFreq = 30;
freq = linspace(1, 30, nFreq);
s = 5 ./ (2*pi.*freq);
% s = linspace(4,12,nFreq) ./ (2*pi.*freq);
fft_data = fft(data,nConv);
tf = zeros(nFreq, nData);
% convolution loop at multiple center frequency
for fi=1:nFreq
% create Morlet wavelet
csine = exp(2*1i*pi*freq(fi)*wtime);
gaus = exp( -(wtime.^2) / (2*s(fi)^2) );
cmw = csine .* gaus;
% compute FFT
cmwX = fft(cmw,nConv);
cmwX = cmwX./max(cmwX); %normalization
% convolution
conv_res = ifft( fft_data .* cmwX );
% trim wings
conv_res = conv_res(nHfkn:end-nHfkn+1);
% power
tf(fi,:) = conv_res.*conj(conv_res); % power
end
% fft of whole data
fft_data = fft(data);
fft_data = fft_data / nData;
fft_pow = fft_data .* conj(fft_data);
hz = linspace(0,srate/2,floor(nData/2)+1);
figure(4), clf
subplot(311)
plot(datatime, data)
xlabel('Time (s)'), ylabel('Amplitude')
title('Time-domain signal')
subplot(312)
plot(hz, fft_pow(1:length(hz)))
set(gca, 'XLim', [0 30], 'YLim', [0 1]*0.5)
subplot(313)
contourf(datatime,freq,tf,40,'linecolor','none')
set(gca,'clim',[0 1], 'ylim', [1 20], 'xlim', [-5 5])
xlabel('Time (s)'), ylabel('Amplitude')
title([ 'Time-frequency power (Hz) '])
% end of part 2
경고: 복소수 X 및/또는 Y 인수의 허수부는 무시됩니다.
part 3: wavelet analysis on real EEG data
% reference: Cha KS, Sunwoo J-S, Byun J-I, et al.
% Working memory deficits in patients with idiopathic restless legs syndrome
% are associated with abnormal theta-band neural synchrony.
% Journal of Sleep Research. 2021;30:e13287.
% paper title
figure(5), clf
imagesc(imread('RLS_sternberg.png'))
part 3-1: convolution method
clear all; close all;
% load data
load C24.mat
srate = EEG.srate;
data = squeeze(EEG.data(48,:,:)); % Pz channel
data1 = reshape(data,1,[]);
nData = length(data1);
% create wavelet
wavtime = -2:1/EEG.srate:2;
nFreq = 50;
frex = linspace(1,50,nFreq);
s = linspace(4,12,nFreq) ./ (2*pi.*frex);
nKern = length(wavtime);
nConv = nData+nKern-1;
nHfkn = floor(nKern/2+1);
fft_eeg = fft(data1,nConv);
tf = zeros(nFreq, length(data));
for fi=1:nFreq
% create Morlet wavelet
csine = exp(2*1i*pi*frex(fi)*wavtime);
gaus = exp( -(wavtime.^2) / (2*s(fi)^2) );
cmw = csine .* gaus;
% compute FFT
cmwX = fft(cmw,nConv);
cmwX = cmwX./max(cmwX);
% convolution
conv_res = ifft( fft_eeg .* cmwX );
% trim wings
conv_res = conv_res(nHfkn:end-nHfkn+1);
% reshape
conv_res = reshape(conv_res, [], EEG.trials);
% power
wav_pow = conv_res.*conj(conv_res);
% power
tf(fi,:) = mean(wav_pow, 2);
end
% baseline correction of tf
base_period = [-500 -200];
baseidx = dsearchn(EEG.times', base_period');
% compute the average per signal
basePow = mean(tf(:,baseidx(1):baseidx(2)),2);
% dB-normalization
tf_db = 10*log10(tf ./ basePow);
% baseline correction of ERP
base_period = [-200 0];
base_erp = mean(data(baseidx(1):baseidx(2), :), 1);
norm_erp = data - base_erp;
erp = mean(norm_erp, 2);
% plot
figure(6), clf
subplot(211)
plot(EEG.times, erp, 'r')
set(gca,'ylim',[-1 1.5]*3, 'xlim', [-200 800])
subplot(212)
contourf(EEG.times,frex,tf_db,40,'linecolor','none')
set(gca,'clim',[0 1]*2, 'ylim', [1 20], 'xlim', [-200 800])
% end of part 3-1
part 3-2: matlab cwt function method
clear all; close all;
% data
load C24.mat
srate = EEG.srate;
data = squeeze(EEG.data(48,:,:)); % Pz channel
data1 = reshape(data,1,[]);
nData = length(data1);
% set cwt function parameter
Ts = 1/srate; % temporal resolution
freq=[1:50]; % frequency resolution
nFreq = length(freq);
wname='cmor1-1'; % complex morlet wavelet as a basis function
Fcenter = centfrq(wname)/Ts; % set center frequency
scale = Fcenter*(1./freq);
coef=cwt(data1,scale,wname); % compute continuous wavelet transform
F = scal2frq(scale,'cmor1-1', Ts);
% get power
pow =abs(coef); % magnitude
pow = reshape(pow, [nFreq, EEG.pnts, EEG.trials]);
tf = squeeze(mean(pow, 3));
% baseline correction of tf
base_period = [-500 -200];
baseidx = dsearchn(EEG.times', base_period');
% compute the average per signal
basePow = mean(tf(:,baseidx(1):baseidx(2)),2);
% dB-normalization
tf_db = 10*log10(tf ./ basePow);
% baseline correction of ERP
base_period = [-200 0];
base_erp = mean(data(baseidx(1):baseidx(2), :), 1);
norm_erp = data - base_erp;
erp = mean(norm_erp, 2);
% plot
figure(7); clf
subplot(211);
plot(EEG.times, erp, 'r');
set(gca,'ylim',[-1 1.5]*3, 'xlim', [-200 800])
title('EEG time series')
subplot(212)
contourf(EEG.times, F, tf_db,40,'linecolor','none')
xlabel('Time (s)'), ylabel('Frequencies (Hz)')
set(gca,'clim',[0 1]*1, 'ylim', [1 20], 'xlim', [-200 800])
title('continuous Morlet wavelet')
% end of part 3-2
% written by Jung, Ki-Young %
Comments (0)
You don't have permission to comment on this page.