diff options
Diffstat (limited to 'detrend.sci')
| -rw-r--r-- | detrend.sci | 100 |
1 files changed, 100 insertions, 0 deletions
diff --git a/detrend.sci b/detrend.sci new file mode 100644 index 0000000..993857e --- /dev/null +++ b/detrend.sci @@ -0,0 +1,100 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) Bruno Pincon +// +// Copyright (C) 2012 - 2016 - Scilab Enterprises +// +// This file is hereby licensed under the terms of the GNU GPL v2.0, +// pursuant to article 5.3.4 of the CeCILL v.2.1. +// This file was originally licensed under the terms of the CeCILL v2.1, +// and continues to be available under such terms. +// For more information, see the COPYING file which you should have received +// along with this program. +// +function [y] = detrend(x, flag, bp) + // + // this function removes the constant or linear or + // piecewise linear trend from a vector x. If x is + // matrix this function removes the trend of each + // column of x. + // + // flag = "constant" or "c" to removes the constant trend + // (simply the mean of the signal) + // flag = "linear" or "l" to remove the linear or piecewise + // linear trend. + // + // To define the piecewise linear trend the function needs the + // breakpoints and these must be given as the third argument bp. + // + // The "instants" of the signal x are in fact from 0 to m-1 + // (m = length(x) if x is a vector and m = size(x,1) in case + // x is a matrix). So bp must be reals in [0 m-1]. + // + + rhs = argn(2) + if rhs < 1 | rhs > 3 then + error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"detrend",1,3)); + elseif rhs == 1 + flag = "linear"; bp = [] + elseif rhs == 2 + bp = [] + end + + if type(x)~=1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",1)); + end + if type(flag)~=10 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",2)); + end + if ~(type(bp)==1 & isreal(bp)) then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",3)); + end + + [mx,nx] = size(x) + if mx==1 | nx==1 then + x_is_vector = %t; x = x(:); m = mx*nx; n = 1 + elseif mx*nx == 0 then + y = []; return + else + x_is_vector = %f; m = mx; n = nx + end + + + if flag == "constant" | flag == "c" then + y = x - ones(m,1)*mean(x,"r") + elseif flag == "linear" | flag == "l" + bp = unique([0 ; bp(:) ; m-1]) + // delete all the breakpoints outside [0,m-1] + while bp(1) < 0, bp(1)=[], end + while bp($) > m-1, bp($)=[], end + // breakpoints are 0-based so add one to + // compare them with signal vector indices (1-based) + bp = bp + 1; + nbp = length(bp); + // build the least square matrix with hat functions + // (as a basis for continuous piecewise linear functions) + A = zeros(m, nbp) + k1 = 1 + delta_bp = diff(bp) + for j = 2:nbp-1 + k2 = ceil(bp(j))-1 + ind = (k1:k2)' + A(ind,j-1) = (bp(j) - ind)/delta_bp(j-1) + A(ind,j) = (ind - bp(j-1))/delta_bp(j-1) + k1 = k2+1 + end + ind = (k1:m)' + A(ind,nbp-1) = (m - ind)/delta_bp(nbp-1) + A(ind,nbp) = (ind - bp(nbp-1))/delta_bp(nbp-1) + // solve the least square pb and retrieve the fitted + // piecewise linear func off the signal + y = x - A*(A\x) + else + error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n") ,.. + "detrend",2,"''constant'',''c'',''linear'',''l''")); + end + + if x_is_vector then + y = matrix(y,mx,nx) + end + +endfunction |