User defined wavelet for continuous time analysis

Hi. I have generated my own user defined wavelet, which I've added using wavemngr. How can I then use this wavelet in a continuous time analysis of my signal? The new cwt will only allow a subset of predefined wavelets. The new cwt function is apparently supposed to run the old cwt if you use the old syntax but I can't get that to work at all. In which case how do you run a continuous time analysis if your mother wavelet is user defined? I feel like I'm missing something here because it seems too obvious a gap!

Thanks for your help,

Katie

Answer 1

[Catherine Davey]

Ok, I finally got it to work, which I'll include here in case other people are having problems. I think there's a few different ways to trigger the use of the old cwt, but one way is to have the first input a cell array stipulating the signal time (either lower & upper bounds, or a vector) and the signal itself;

>> sig = {voltage, time}; % my example has extracellular recording of neuron voltage

I put my wavelet values in to 'wav', rather than calling on the wavemngr using its shortname, and put my required scales into a positive vector called 'scales'. There may be other ways to do it, but I got it to work by then calling:

>> coeff = cwt({voltage, time}, scales, wav);

[Md. Atiqur Rahman]

Could you show me the code at [e-mail] please ? It would help me a lot because I am working on it in my PhD study.

[César Alcaide León]

I am also intereseted in that code

Please, should you send me to [e-mail]?

It is for my PhD thesis, thanks

[Catherine Davey]

Hey. I'm really sorry Atiqur Rahman, I didn't notice your request for the code. I have many PhD students and know how much they would appreciate this sort of support. So I'm sorry about that.

César - I'm not sure that my code will help that much, but I've copied a snippet here. Happy to email you more - the code was written for a GUI so there's a lot of it! I've just added the important bits here.

To run it:

% make mother wavelets for each template

APwavs = generateNewSETWavelets( rescaleAPs, 'continuous', 'polynomial', 'aicc', true, true );

addSETUserMotherWavelet( APwavs );

wav = APwavs{1}.mother;

% get centre freq of wavelet - seems to be determined by the freq of the

% best fit sin wave, & the sampling resolution

% fc = num sin waves that fit mother wavelet / ( time(end) - time(1) )

fcentre = centfrq('ap1',2,'plot'); % centre freq of wavelet

% scal2freq returns the pseudo-frequencies corresponding to the scales

% and wavelet, where fscale = fcentre / (scales * dt)

% - therefore: scale = fcentre / (fscales * dt)

% - we're interested in freqs less than 1000

figure;

scales = round( fcentre ./ ((1000:-100:100) * dt ), 'plot' );

fscales = scal2frq(scales, 'ap1', dt);

fscales = fcentre ./ (scales * dt);

coeff = cwt({voltage, time}, scales, wav);

% % addSETUserMotherWavelet( APwavs )

%

% Generate user mother wavelets to fit AP templates using

% generateNewSETWavelets, and add them to the wavelet toolbox here

function addSETUserMotherWavelet( APwavs )

% APwavs has the following fields for each AP template

% APwavs{ti}.mother

% APwavs{ti}.pred

% APwavs{ti}.degree

numwav = length( APwavs );

for wi=1:numwav

wav = APwavs{wi}.mother;

deg = APwavs{wi}.degree;

time = APwavs{wi}.time;

% wavelet has to be even else it's zero padded

if ~iseven( length(wav) )

% repeat last value rather than adding 0 just in case it causes a

% big jump

wav = [wav(:); wav(end)];

end

% longname can be whatevs I think

longname = ['APwave' num2str(wi)];

% shortname must be 4 characters or less

shortname = ['ap' num2str(wi)];

fname = [shortname '.mat'];

% need lower & upper bound of time vector

lb = length(wav)/2 - length(wav) + 1;

ub = length(wav)/2;

% must create a wavelet .mat file with the same name as the wavelet's

% shortname, and a variable in the file named the same as shortname

% - for multiple wavelets in the family it's a bit different i think

eval( [shortname '= wav;'] );

% time = ( lb : ub )';

eval( ['X= time;'] );

eval( ['Y= wav;'] );

save( fname, shortname );

save( fname, 'X', 'Y' );

out = wavemngr('read'); % get list of current wavelets

% if wavelet already exists gotta delete before adding new one

% - compare with shortname + space so we don't mistake ap10 for ap1,

% for example

if contains( toVec(out')', [shortname ' '] )

wavemngr('del', shortname);

end

wavemngr('add', longname, shortname, 4, '', fname, [lb ub]);

% List wavelets families.

wavemngr('read');

end

end

% % [APwavs, APpred] = generateNewSETWavelets( APtemplates, [regularity], [fitmethod], [goftest], [zeromean], [doplot] );

%

% Construct mother wavelets from AP templates so we can later do continuous

% time wavelet analysis. Note that all wavelets are required to integrate

% to zero so a zero meaned version of the AP templates may fit better.

%

% Inputs:

% APtemplates - matrix of time x template

% regularity - ensure mother wavelets are 'differentiable' or 'continuous'

% fitmethod - fit using 'polynomial' or 'orthconst'

% goftest - evaluate gof of mother wavelet to AP template using

% 'aic', 'aicc', 'bic', 'fpe', 'hqc', 'mse' (no param penalty)

% zeromean - mother wavelets seem to fit better if AP is zero meaned,

% but I'm not sure how this will later impact its use

% doplot - plot best fit user defined wavelet based on goftest for each AP template

% Outputs:

% APwavs{i}.mother - cell array of best fit mother wavelet for each template

% APwavs{i}.pred - cell array of prediction of mother wavelet for each template

% APwavs{i}.degree - degree of polynomial or orthogonal constants

%

% See also: pat2cwav calcGOF

function APwavs = generateNewSETWavelets( APtemplates, varargin )

if nargin==0, help generateNewSETWavelets; return; end

regularity = 'differentiable'; % 'continuous'; %

fitmethod = 'polynomial'; % 'orthconst'; %

goftest = 'mse'; % 'aic'; % 'aicc'; % 'bic'; % 'hqc'; % 'fpe'; % 'mse'; %

zeromean = true;

doplot = true;

optargs = {regularity, fitmethod, goftest, zeromean, doplot};

narg = length( varargin );

optargs(1:narg) = varargin;

[regularity, fitmethod, goftest, zeromean, doplot] = optargs{:};

[nS, nT] = size( APtemplates.data );

t = linspace(0, 1, nS);

% APtemplates: struct with fields:

% type: 'ap'

% data: [28×8 single]

% time: [28×1 double]

% dt: 2.0000e-04

% params: [1×1 struct]

% APfamily: {1×8 cell}

% name: 'Voltage_APs'

APwavs = cell(0);

if ~any( strcmpi( regularity, {'differentiable', 'continuous'} ) )

str = sprintf('\tgenerateNewSETWavelets: unknown regularity (set to differentiable or continuous)\n');

cprintf( 'Errors*', str );

return;

end

if ~any( strcmpi( goftest, {'aic', 'bic', 'aicc', 'fpe', 'hqc', 'mse'} ) )

str = sprintf('\tgenerateNewSETWavelets: unknown gof fit (set to aic, aicc, bic, fpe, hqc or mse)\n');

cprintf( 'Errors*', str );

return;

end

if ~any( strcmpi( fitmethod, {'polynomial', 'orthconst', 'discpoly', 'lagrange', 'cos'} ) )

str = sprintf('\tgenerateNewSETWavelets: unknown fit method (set to polynomial or orthconst)\n');

cprintf( 'Errors*', str );

return;

end

maxpoly = floor( nS / 2 );

minpoly = ternaryOp( strcmpi( regularity, 'continuous' ), 3, 5 );

npoly = maxpoly - minpoly + 1;

APwavs = cell(nT, 1); APpred = cell(nT, 1);

malloc = @(ndim) zeros([ndim nT]); % dims x num_templates

mse = malloc(npoly);

gof = malloc(npoly);

wav = malloc([nS npoly]);

pred = malloc([nS npoly]);

for ti=1:nT

y = APtemplates.data(:,ti);

if zeromean

y = y - mean(y);

end

for pp=1:npoly

polydeg = minpoly + pp - 1;

[Y,~,nc] = pat2cwav(y(:)', fitmethod, polydeg, regularity) ;

wav(:,pp,ti) = Y(:);

pred(:,pp,ti) = Y(:)*nc;

% calc mse & gof to see which polynomial degree we should use

mse(pp,ti) = sum((y(:)-Y(:)*nc).^2);

% AP template built from avging so make length a little longer

gof(pp,ti) = calcGOF( goftest, mse(pp,ti), polydeg, length(y)*4, false );

end

% I'm finding that gof is too influence by number of samples coz

% there's so few (about 28), so perhaps use mse instead?

[~, best] = min( gof(:,ti) );

APwavs{ti}.mother = wav( :,best,ti);

APwavs{ti}.pred = pred(:,best,ti);

APwavs{ti}.degree = best + minpoly - 1;

APwavs{ti}.time = APtemplates.time; % shouldn't copy for each but easiest for now

end

if doplot

figure;

nr = ceil(sqrt(nT)); nc = ceil(nT/nr);

for ti=1:nT

y = APtemplates.data(:,ti);

if zeromean

y = y - mean(y);

end

subplot(nr,nc,ti)

hold on;

plot(t, APwavs{ti}.pred);

plot(t, y, 'k-', 'linewidth', 2);

str = sprintf('AP %d (deg %d)', ti, APwavs{ti}.degree);

title( str );

end

end

end

% save('mother.mat', 'Y', 'X');