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);