User defined wavelet for continuous time analysis

6 次查看（过去 30 天）
Catherine Davey
Catherine Davey 2018-11-21
评论： Catherine Davey 2020-5-11
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

请先登录，再回答此问题。

采纳的回答

Catherine Davey
Catherine Davey 2018-11-22
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);
  3 个评论
显示 1更早的评论
César Alcaide León
César Alcaide León 2020-5-10
I am also intereseted in that code
Please, should you send me to cesaralcaide@gmail.com?
It is for my PhD thesis, thanks
Catherine Davey
Catherine Davey 2020-5-11
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');

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Continuous Wavelet Transforms 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by