Merge pull request #656 from aodn/spike_qc_nonburst

Spike qc nonburst - part 1

AutomaticQC/SpikeClassifiers/SavGol.m

@@ -0,0 +1,42 @@
+function g = SavGol (f, nl, nr, M)
+
+% SAVGOL SavGol smoothes the data in the vector f by means of a
+%        Savitzky-Golay smoothing filter.
+%
+%        Input: f : noisy data
+%        nl: number of points to the left of the reference point
+%        nr: number of points to the right of the reference point
+%        M : degree of the least squares polynomial
+%
+%        Output: g: smoothed data
+%
+%        W. H. Press and S. A. Teukolsky,
+%        Savitzky-Golay Smoothing Filters,
+%        Computers in Physics, 4 (1990), pp. 669-672.
+
+% matrix A
+A = ones (nl+nr+1, M+1);
+for j = M:-1:1,
+  A (:, j) = [-nl:nr]' .* A (:, j+1);
+end
+
+% filter coefficients c
+[Q, R] = qr (A);
+c = Q (:, M+1) / R (M+1, M+1);
+
+% smoothing of the noisy data
+% Note that there are two equivalent ways to apply the Savitzky-Golay
+% filter to the vector f.  In the first case we use a for-loop whereas
+% in the second case we use the faster built-in function filter.
+%
+% g = f;
+% n = size (f);
+% for i = 1+nl:n-nr,
+%   g (i) = c' * f (i-nl:i+nr);
+% end
+%
+n = length (f);
+g = filter (c (nl+nr+1:-1:1), 1, f);
+g (1:nl) = f (1:nl);
+g (nl+1:n-nr) = g (nl+nr+1:n);
+g (n-nr+1:n) = f (n-nr+1:n);

AutomaticQC/SpikeClassifiers/imosSpikeClassifierHampel.m

@@ -0,0 +1,57 @@
+function [spikes, fsignal] = imosSpikeClassifierHampel(signal, window, stdfactor)
+% function spikes, fsignal = imosSpikeClassifierHampel(signal, window, stdfactor)
+%
+% A wrapper to Hampel, a median windowed filter, for both
+% Burst and Equally spaced time series.
+%
+% The hampel method is a robust parametric denoising filter that uses
+% the MAD deviation (scaled) from the median in a pre-defined
+% indexed centred window.
+%
+% See also hampel.m.
+%
+% Inputs:
+%
+% signal - A 1-d signal.
+% window - A window index range. If signal is in burst mode, this is ignored.
+% stdfactor - A multiplying scale for the standard deviation.
+%
+% Outputs:
+%
+% spikes - An array with spikes indexes.
+% fsignal - A filtered signal where the spikes are subsituted
+%           by the median of the window.
+%
+% Example:
+%
+%
+% author: [email protected]
+%
+
+% 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>.
+%
+
+if nargin < 3
+    window = 3;
+    stdfactor = 10;
+elseif nargin < 2
+    stdfactor = 10;
+end
+
+[fsignal, bind, ~, ~] = hampel(signal, window, stdfactor);
+spikes = find(bind);
+end

AutomaticQC/SpikeClassifiers/imosSpikeClassifierNonBurstSavGolOTSU.m

@@ -0,0 +1,111 @@
+function [spikes] = imosSpikeClassifierNonBurstSavGolOTSU(signal, window, pdeg, nbins, oscale)
+% function spikes = imosSpikeClassifierNonBurstSavGolOTSU(signal, window, pdeg, nbins, oscale)
+%
+% Detect spikes in signal by using the Otsu threshold method
+% applied to the noise estimated by a Savitzy Golay Filter applied on the signal.
+%
+% The routine uses the SavGol.m implementation provided by
+% Ken Ridgway.
+%
+% Inputs:
+%
+% signal - the signal.
+%  window - the SavGol window argument (odd integer).
+%  pdeg - the SavGol lsq polynomial order (integer).
+%  nbins - the number of histogram bins in the OTSU threshold method (positive integer)
+%  oscale - a constant scale factor to weight the OTSU threshold method (positive number)
+%
+% Outputs:
+%
+% spikes - spikes indexes.
+%
+% Example:
+%
+% t = 0:3600:86400*5;
+% omega = 2*pi/86400;
+% signal = 10*sin(omega*t);
+% signal([20]) = max(signal).^3;
+% [spikes] = imosSpikeClassifierNonBurstSavGolOTSU(signal,5,2,100,5);
+% assert(spikes,20)
+%
+% author: Unknown
+% author: Ken Ridgway
+% author: [email protected] [rewrite/refactoring]
+
+%
+% 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>.
+%
+
+if nargin == 1
+    window = 5; %minimum filtering
+    pdeg = 2; %quadratic polynomial
+    nbins = 256; % default otsu bins
+    oscale = 1; %arbitrary scale used to reduce the otsu threshold
+elseif nargin == 2
+    pdeg = 2; %quadratic polynomial
+    nbins = 256; % default bins in despiking1
+    oscale = 1; %arbitrary scale used to reduce the otsu threshold
+elseif nargin == 3
+    nbins = 256; % default bins in despiking1
+    oscale = 1; %arbitrary scale used to reduce the otsu threshold
+elseif nargin == 4
+    oscale = 1; %arbitrary scale used to reduce the otsu threshold
+end
+
+if window <= 3
+    error('Window need to be even and larger than 3')
+else
+
+    if rem(window, 2) == 0
+        window = window - 1;
+    end
+
+end
+
+fsignal = SavGol(signal, (window - 1) / 2, (window - 1) / 2, pdeg);
+noise = signal - fsignal;
+abs_noise = abs(noise);
+otsu = otsu_threshold(abs_noise, nbins);
+threshold = otsu / oscale;
+ring_spikes = find(abs_noise > threshold | abs_noise <- threshold);
+spikes = centred_indexes(ring_spikes);
+
+end
+
+function [cgind] = centred_indexes(sind)
+% Return the central value of a group of
+% consecutive integers in a sorted array of integers.
+
+s_start = sind(1);
+s_end = sind(end);
+cgind = sind*NaN;
+c=1;
+for k = 1:length(sind) - 1
+    adjacent = sind(k + 1) - sind(k) <= 1;
+    if adjacent
+        continue
+    end
+
+    s_end = sind(k);
+    cgind(c) = floor((s_end + s_start) / 2);
+    s_start = sind(k + 1);
+    c = c+1;
+end
+cgind(end) = floor((s_end + s_start) / 2);
+
+cgind = cgind(~isnan(cgind));
+end

AutomaticQC/SpikeClassifiers/imosSpikeClassifierOTSU.m

@@ -0,0 +1,114 @@
+function [spikes,threshold] = imosSpikeClassifierOTSU(signal, nbins, oscale, centralize)
+% function [spikes,threshold] = imosSpikeClassifierOTSU(signal, nbins, oscale, centralize)
+%
+% Detect spikes in by using the Otsu threshold method.
+%
+% Inputs:
+%
+%  signal - the signal.
+%  nbins - the number of histogram bins in the OTSU threshold method (positive integer)
+%  oscale - a constant scale factor to weight the OTSU threshold method (positive number)
+%  centralize - a boolean to average/centralize the detected indexes (default: true)
+%
+% Outputs:
+%
+% spikes - spikes indexes.
+% threshold - the threshold used for spike detection.
+%
+% Example:
+%
+% signal = [0,1,0];
+% [spikes] = imosSpikeClassifierOTSU(signal,3,1,false);
+% assert(spikes==2)
+%
+% %centralized
+% signal = [0,1,1,1,0]
+% [spikes] = imosSpikeClassifierOTSU(signal,3,1,false);
+% assert(spikes==[4,5,6])
+% [spikes] = imosSpikeClassifierOTSU(signal,3,1,true);
+% assert(spikes==[5])
+%
+% %harmonic
+% omega = 2*pi/86400;
+% t = 0:3600:86400*5;
+% signal = 10*sin(omega*t)
+% signal([10,20]) = max(signal)^3;
+% [spikes] = imosSpikeClassifierOTSU(signal,100,1,true);
+% assert(spikes==[10,20])
+%
+% %
+% omega = 2*pi/86400;
+% t = 0:3600:86400*5;
+% signal = 10*sin(omega*t)
+% signal(20) = signal(20)*2;
+% signal(21:end) = signal(21:end)+signal(20);
+% signal([33,55]) = sign(signal([33,55])).*signal([33,55]).^2;
+% [spikes] = imosSpikeClassifierOTSU(signal,100,1,true);
+% assert(spikes==[33,55])
+%
+%
+% author: [email protected] [rewrite/refactoring]
+
+%
+% 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>.
+%
+
+if nargin == 1
+    nbins = 256;
+    oscale = 1;
+    centralize = true;
+elseif nargin == 2
+    oscale = 1;
+    centralize = true;
+elseif nargin == 3
+    centralize = true;
+end
+
+otsu = otsu_threshold(signal, nbins);
+threshold = otsu / oscale;
+spikes = find(signal > abs(threshold) | signal < -1*abs(threshold));
+if centralize
+    spikes = centred_indexes(spikes);
+end
+
+end
+
+function [cgind] = centred_indexes(sind)
+% Return the central value of a group of
+% consecutive integers in a sorted array of integers.
+
+s_start = sind(1);
+cgind = sind * NaN;
+c = 1;
+
+for k = 1:length(sind) - 1
+    adjacent = sind(k + 1) - sind(k) <= 1;
+
+    if adjacent
+        continue
+    end
+
+    s_end = sind(k);
+    cgind(c) = floor((s_end + s_start) / 2);
+    s_start = sind(k + 1);
+    c = c + 1;
+end
+
+s_end = sind(end);
+cgind(end) = floor((s_end + s_start) / 2);
+cgind = cgind(~isnan(cgind));
+end