feat(qc): Include new spikeQC Burst functions

Provide 2 spikeQC functions that works on burst data,
along with their default configurations.

a. BurstHampel is Hampel considering burst sampling
b. BurstRunningStats is a wrapper to apply any statistical
function of scale and dispersion along bursts.

AutomaticQC/SpikeClassifiers/imosSpikeClassifierBurstHampel.m

@@ -0,0 +1,176 @@
+function [spikes] = imosSpikeClassifierBurstHampel(bindexes, signal, use_burst_window, half_window_width, madfactor, lower_mad_limit, repeated_only)
+% function [spikes] = imosSpikeClassifierBurstHampel(bindexes, signal,  use_burst_window,  half_window_width, madfactor, lower_mad_limit, repeated_only)
+%
+% A Hampel classifier for burst data. The half_window_width is applied at burst level instead of
+% sample level.
+%
+% Inputs:
+%
+% bindexes - a cell with the invidivudal burst range indexes.
+% signal - A Burst 1-d signal.
+% use_burst_window  - a boolean to consider the half_window_width applied at burst scale. If false, will consider the burst series as continuous.
+% half_window_width - The half_window_width burst range (the full window size will be 1+2*half_window_width)
+% madfactor -  A multipling scale for the MAD.
+% repeated_only - a boolean to mark spikes in burst only if they are detected more than one time (only for half_window_width>0).
+% lower_mad_limit - a lower threshold for the MAD values, which values below will be ignored.
+%
+% Outputs:
+%
+% spikes - An array with spikes indexes.
+% fsignal - A filtered signal where the spikes are substituted by median values of the window.
+%
+% Example:
+%
+% % simple spikes
+% x = randn(1,100)*1e-2;
+% spikes = [3,7,33,92,99]
+% x(spikes) = 1000;
+% bduration = 6;
+% v = 1:bduration:length(x)+bduration;
+% for k=1:length(v)-1;
+% bind{k} = [v(k),min(v(k+1)-1,length(x))];
+% end
+% [dspikes] = imosSpikeClassifierBurstHampel(bind,x);
+% assert(isequal(dspikes,spikes));
+% % equal to Hampel
+% [dspikes] = imosSpikeClassifierBurstHampel(bind,x,length(bind));
+% [dspikes2] = imosSpikeClassifierHampel(x,length(bind));
+% assert(isequal(dspikes,spikes));
+% assert(isequal(dspikes,dspikes2));
+%
+% % detecting entire burst as spike
+% x = randn(1,100)*1e-2;
+% fullburst_spiked = [3,7,13,14,15,16,17,18,33,92,99]
+% x(fullburst_spiked) = 1000;
+% [dspikes] = imosSpikeClassifierBurstHampel(bind,x);
+% assert(isequal(dspikes,[3,7,33,92,99]));
+% [dspikes2] = imosSpikeClassifierBurstHampel(bind,x,true,2);
+% assert(isequal(dspikes2,fullburst_spiked))
+%
+%
+% author: hugo.oliveira@utas.edu.au
+%
+
+% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated
+% Marine Observing System (IMOS).
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation version 3 of the License.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program.
+% If not, see <https://www.gnu.org/licenses/gpl-3.0.en.html>.
+%
+
+narginchk(2, 7)
+
+if nargin == 2
+    use_burst_window = 1;
+    half_window_width = 0;
+    madfactor = 10;
+    lower_mad_limit = 0;
+    repeated_only = false;
+elseif nargin == 3
+    half_window_width = 0;
+    madfactor = 10;
+    lower_mad_limit = 0;
+    repeated_only = false;
+elseif nargin == 4
+    madfactor = 10;
+    lower_mad_limit = 0;
+    repeated_only = false;
+elseif nargin == 5
+    lower_mad_limit = 0;
+    repeated_only = false;
+elseif nargin == 6
+    repeated_only = false;
+end
+
+if ~use_burst_window
+    if half_window_width == 0
+        half_window_width = 1;
+    end
+    spikes = imosSpikeClassifierHampel(signal,half_window_width,madfactor,lower_mad_limit);
+    return
+end
+
+nbursts = double(half_window_width) .* 2 + 1;
+burst_is_series = nbursts >= length(bindexes);
+if burst_is_series
+    spikes = imosSpikeClassifierHampel(signal, ceil(length(signal)/2+1), madfactor,lower_mad_limit);
+    return
+end
+
+spikes = NaN(size(signal));
+%left bry
+bs = 1;
+range = 1:bindexes{nbursts}(end);
+[bspikes] = imosSpikeClassifierHampel(signal(range), length(range), madfactor,lower_mad_limit);
+
+if ~isempty(bspikes)
+    blen = length(bspikes);
+    spikes(bs:bs + blen - 1) = bspikes;
+    bs = bs + blen;
+end
+
+
+series_end_at_next_window = nbursts+1 >= length(bindexes)-nbursts;
+if series_end_at_next_window
+    bie = bindexes{nbursts+1}(1);
+else
+    sburst = 1;
+    ni = nbursts+1;
+    ne = length(bindexes)-nbursts;
+    h=waitbar(0,'Computing Spikes');
+    wbval=0;
+    wbint=0.05;
+    for eburst = ni:ne
+        if eburst/ne > (wbval+wbint)
+            wbval = wbval+wbint;
+            waitbar(eburst/ne,h,'Computing Spikes');
+        end
+        sburst = sburst + 1;
+        bis = bindexes{sburst}(1);
+        bie = bindexes{eburst}(end);
+        csignal = signal(bis:bie);
+        [bspikes] = imosSpikeClassifierHampel(csignal, length(csignal), madfactor,lower_mad_limit);
+
+        if ~isempty(bspikes)
+            blen = length(bspikes);
+            spikes(bs:bs + blen - 1) = bspikes + bis - 1;
+            bs = bs + blen;
+        end
+
+    end
+end
+if ~isempty(h)
+    h.delete()
+end
+%right bry
+range = bie + 1:bindexes{end}(end);
+[bspikes] = imosSpikeClassifierHampel(signal(range), length(range), madfactor,lower_mad_limit);
+
+if ~isempty(bspikes)
+    blen = length(bspikes);
+    spikes(bs:bs + blen - 1) = bspikes + bie;
+end
+
+spikes = spikes(~isnan(spikes));
+if repeated_only
+    srange = min(spikes)-1:max(spikes)+1;
+    candidates = find(histcounts(spikes,srange)>1);
+    if ~isempty(candidates)
+        spikes = candidates;
+    end
+else
+    spikes = unique(spikes);
+end
+
+
+end

AutomaticQC/SpikeClassifiers/imosSpikeClassifierBurstRunningStats.m

@@ -0,0 +1,74 @@
+function [spikes] = imosSpikeClassifierBurstRunningStats(bindexes, signal, scalefun, dispersionfun, dispersion_factor)
+% function [spikes] = imosSpikeClassifierBurstRunningStats(bindexes, signal, scalefun, dispersionfun, dispersion_factor)
+%
+% A simple thresholding method based on user defined deviations around a central scale. Values above
+% the dispersion will be tagged.
+%
+% Inputs:
+%
+% bindexes - a cell with the invidivudal burst range indexes.
+% signal - A Burst 1-d signal.
+% scalefun - the function to estimate the centre of the distribution. (e.g. mean or median)
+%          Default: @mean
+% dispersionfun - the function to estimate the dispersion (e.g. std, mean_absolute_deviation, median_absolute_deviation)
+%               Default: @std
+% dispersion_factor - the multiplying factor for the dispersion.
+%                   Default: 2
+%
+% Outputs:
+%
+% spikes - An array with spikes indexes.
+% cscale - the computed scale
+% cdispersion - the computed dispersion (scaled).
+%
+% Example:
+% z = randn(1,10);
+% z(10) = 10000;
+% [spikes] = imosSpikeClassifierRunningStats(z,@mean,@std,2);
+% assert(spikes(10)==1)
+% assert(spikes(1:9)==0)
+%
+% author: hugo.oliveira@utas.edu.au
+%
+
+% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated
+% Marine Observing System (IMOS).
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation version 3 of the License.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program.
+% If not, see <https://www.gnu.org/licenses/gpl-3.0.en.html>.
+%
+
+%detect bursts
+
+spikes = NaN(size(signal));
+ne = length(bindexes);
+
+h=waitbar(0,'Computing Spikes');
+wbval=0;
+wbint=0.05;
+for k=1:ne
+	if k/ne > (wbval+wbint)
+		wbval=wbval+wbint;
+		waitbar(k/ne,h,'Computing Spikes');
+	end
+	brange = bindexes{k};
+	xs = brange(1);
+	xe = brange(2);
+	bspikes = imosSpikeClassifierRunningStats(signal(xs:xe),scalefun,dispersionfun,dispersion_factor);
+	if ~isempty(bspikes)
+		spikes(xs:xs+length(bspikes)-1) = bspikes+xs-1;
+	end
+end
+h.delete();
+spikes = unique(spikes(~isnan(spikes)));
+end

AutomaticQC/SpikeClassifiers/imosSpikeClassifierRunningStats.m

@@ -0,0 +1,53 @@
+function [spikes] = imosSpikeClassifierRunningStats(signal, scalefun, dispersionfun, dispersion_factor)
+% function [spikes] = imosSpikeClassifierRunningStats(signal, scalefun, dispersionfun, dispersion_factor)
+%
+% A simple thresholding method based on user defined deviations around a central scale. Values above
+% the dispersion will be tagged.
+%
+% Inputs:
+%
+% signal - A 1-d signal.
+% scalefun - the function to estimate the centre of the distribution. (e.g. mean or median)
+%          Default: @mean
+% dispersionfun - the function to estimate the dispersion (e.g. std, mean_absolute_deviation, median_absolute_deviation)
+%               Default: @std
+% dispersion_factor - the multiplying factor for the dispersion.
+%                   Default: 2
+%
+% Outputs:
+%
+% spikes - An array with spikes indexes.
+% cscale - the computed scale
+% cdispersion - the computed dispersion (scaled).
+%
+% Example:
+% z = randn(1,10);
+% z(10) = 10000;
+% [spikes] = imosSpikeClassifierRunningStats(z,@mean,@std,2);
+% assert(spikes(10)==1)
+% assert(spikes(1:9)==0)
+%
+% author: hugo.oliveira@utas.edu.au
+%
+
+% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated
+% Marine Observing System (IMOS).
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation version 3 of the License.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program.
+% If not, see <https://www.gnu.org/licenses/gpl-3.0.en.html>.
+%
+cscale = scalefun(signal);
+cdispersion = dispersion_factor * dispersionfun(signal);
+above_and_below = (signal > (cscale + cdispersion)) + (signal < (cscale - cdispersion));
+spikes = find(above_and_below);
+end

AutomaticQC/imosTimeSeriesSpikeQCBurst.txt

@@ -0,0 +1,26 @@
+%If auto is on, use the default options from named method
+auto_function = burst_hampel
+
+% burst Hampel filter
+burst_hampel_function = imosSpikeClassifierBurstHampel
+% boolean use windows as bursts instead of samples
+burst_hampel_use_burst_window = 0
+% The number of bursts to consider in a Hampel/Median window
+burst_hampel_half_window_width = 60
+%multiplicative factor for the mad within the window
+burst_hampel_madfactor = 5
+%multiplicative factor for the mad within the window
+burst_hampel_lower_mad_limit = 0
+%number of burst to aggregate as a window in the filter
+%detect all not only repeated
+burst_hampel_repeated_only = 0
+
+%Running Statistics
+burst_runningstats_function = imosSpikeClassifierBurstRunningStats
+
+% the scale statistic function
+burst_runningstats_scalefun = nanmedian
+% the dispersion statistic function
+burst_runningstats_dispersionfun = mad
+% the dispersion scale
+burst_runningstats_dispersion_factor = 2