1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) Bruno Pincon
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
//
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
|
