Wavelet analysis toolbox
The zip file contains the following functions (see their help text for input / output variables):
wavelets.m: returns the amplitude time-course of a signal, or the amplitudes and phase-difference time courses of two signals.
getERA.m, getPLV.m: take amplitude or phase-difference vectors and find the average trial time course (Event-Related Amplitude and the Phase-Locking Value, respectively). Input includes a vector of trial start-times and the samples per trial.
test_wavelets.m: example file using simulate_data.m for data.
Code of wavelets.m
function varargout = wavelets(varargin) % function amp = wavelets(Fs, freq_mid, freq_sd, vector1) % function [amp1, amp2, phase_difference] = wavelets(Fs, freq_mid, freq_sd, vector1, vector2) % % T.E.Gladwin, 2007 Fs = varargin{1}; freq_mid = varargin{2}; freq_sd = varargin{3}; vector1 = varargin{4}; if length(varargin) > 4, vector2 = varargin{5}; else vector2 = []; end; [realPart1, imPart1] = getComplexConvolutions(vector1, Fs, freq_mid, freq_sd); varargout{1} = CC2Amp(realPart1, imPart1); if ~isempty(vector2), [realPart2, imPart2] = getComplexConvolutions(vector2, Fs, freq_mid, freq_sd); varargout{2} = CC2Amp(realPart2, imPart2); varargout{3} = CC2PD(realPart1, imPart1, realPart2, imPart2); end; % Subfunctions function [realPart, imPart] = getComplexConvolutions(signal, Fs, freq_mid, freq_sd) a = 1; b = my_Morlet('cos', freq_mid, freq_sd, Fs); realPart = filter(b,a,signal) * (1 / Fs); b = my_Morlet('sin', freq_mid, freq_sd, Fs); imPart = filter(b,a,signal) * (1 / Fs); shifter = ceil(length(b) / 2); realPart(1:(end - shifter + 1)) = realPart(shifter:end); realPart(1:shifter) = 0; realPart((end - shifter):end) = 0; imPart(1:(end - shifter + 1)) = imPart(shifter:end); imPart(1:shifter) = 0; imPart((end - shifter):end) = 0; function wl = my_Morlet(type, freq_mid, freq_sd, Fs); par = 7; sigT = 1 / (2 * pi * freq_sd); t = (-par * sigT):(1 / Fs):(par * sigT); if strcmpi(type, 'cos') == 1, osc = cos(2 * pi * freq_mid * t); else, osc = sin(2 * pi * freq_mid * t); end; g = Gaussian(t, sigT); wl = g .* osc; wl = wl * 2 / (sigT * sqrt(2 * pi)); function g = Gaussian(t, sigT); g = exp((-t .^ 2) / (2 * sigT ^ 2)); function amp = CC2Amp(realPart, imPart) amp = sqrt(imPart .^ 2 + realPart .^ 2); function angle_diff = CC2PD(realPart_A, imPart_A, realPart_B, imPart_B) % get the phases angle_A = atan2(imPart_A, realPart_A); angle_B = atan2(imPart_B, realPart_B); % get the phase differences angle_diff = angle_A - angle_B;
Code of getERA.m
function ERA = getERA(amp, trial_starts, samples_per_trial) % ERA = getERA(amp, nTrials) % % T.E.Gladwin, 2007 % extract trials memory1 = amp(:); amp = []; for n = 1:length(trial_starts), a = trial_starts(n); b = a + samples_per_trial - 1; if a > 0 & b < length(memory1), amp = [amp memory1(a:b)]; end; end; % Reshape the vector into a trial-length x nTrials matrix, and average over % trials (== over columns) ERA = mean(amp, 2);
Code of getPLV.m
function PLV = getPLV(angle_diff, trial_starts, samples_per_trial) % PLV = getPLV(PD, trial_starts, samples_per_trial) % % T.E.Gladwin, 2007 % extract trials memory1 = angle_diff(:); angle_diff = []; for n = 1:length(trial_starts), a = trial_starts(n); b = a + samples_per_trial - 1; if a > 0 & b < length(memory1), angle_diff = [angle_diff; memory1(a:b)]; end; end; % Get the x and y coordinates of the unit phase-difference vectors x = cos(angle_diff); y = sin(angle_diff); % Get the vector sum (sum over columns), get the length, and finally % normalize it to get the instantaneous PLV. % Reshape the vectors into trial-length x nTrials matrices x = reshape(x, samples_per_trial, length(x) / samples_per_trial); y = reshape(y, samples_per_trial, length(y) / samples_per_trial); x_sum = sum(x, 2); y_sum = sum(y, 2); sum_length = sqrt(x_sum .^ 2 + y_sum .^ 2); PLV = sum_length / length(trial_starts);