X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Fsgolay.m;fp=octave_packages%2Fsignal-1.1.3%2Fsgolay.m;h=a35e76872fc1738efdae30633fd08dfb4ac62cf9;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/signal-1.1.3/sgolay.m b/octave_packages/signal-1.1.3/sgolay.m
new file mode 100644
index 0000000..a35e768
--- /dev/null
+++ b/octave_packages/signal-1.1.3/sgolay.m
@@ -0,0 +1,104 @@
+## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
+## Copyright (C) 2004 Pascal Dupuis <Pascal.Dupuis@esat.kuleuven.ac.be>
+##
+## 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; either version 3 of the License, or (at your option) any later
+## version.
+##
+## 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 <http://www.gnu.org/licenses/>.
+
+## F = sgolay (p, n [, m [, ts]])
+##   Computes the filter coefficients for all Savitzsky-Golay smoothing
+##   filters of order p for length n (odd). m can be used in order to
+##   get directly the mth derivative. In this case, ts is a scaling factor. 
+##
+## The early rows of F smooth based on future values and later rows
+## smooth based on past values, with the middle row using half future
+## and half past.  In particular, you can use row i to estimate x(k)
+## based on the i-1 preceding values and the n-i following values of x
+## values as y(k) = F(i,:) * x(k-i+1:k+n-i).
+##
+## Normally, you would apply the first (n-1)/2 rows to the first k
+## points of the vector, the last k rows to the last k points of the
+## vector and middle row to the remainder, but for example if you were
+## running on a realtime system where you wanted to smooth based on the
+## all the data collected up to the current time, with a lag of five
+## samples, you could apply just the filter on row n-5 to your window
+## of length n each time you added a new sample.
+##
+## Reference: Numerical recipes in C. p 650
+##
+## See also: sgolayfilt
+
+## Based on smooth.m by E. Farhi <manuf@ldv.univ-montp2.fr>
+
+function F = sgolay (p, n, m = 0, ts = 1)
+
+  if (nargin < 2 || nargin > 4)
+    print_usage;
+  elseif rem(n,2) != 1
+    error ("sgolay needs an odd filter length n");
+  elseif p >= n
+    error ("sgolay needs filter length n larger than polynomial order p");
+  else
+    if length(m) > 1, error("weight vector unimplemented"); endif
+
+    ## Construct a set of filters from complete causal to completely
+    ## noncausal, one filter per row.  For the bulk of your data you
+    ## will use the central filter, but towards the ends you will need
+    ## a filter that doesn't go beyond the end points.
+    F = zeros (n, n);
+    k = floor (n/2);
+    for row = 1:k+1
+      ## Construct a matrix of weights Cij = xi ^ j.  The points xi are
+      ## equally spaced on the unit grid, with past points using negative
+      ## values and future points using positive values.
+      C = ( [(1:n)-row]'*ones(1,p+1) ) .^ ( ones(n,1)*[0:p] );
+      ## A = pseudo-inverse (C), so C*A = I; this is constructed from the SVD 
+      A = pinv(C);
+      ## Take the row of the matrix corresponding to the derivative
+      ## you want to compute.
+      F(row,:) = A(1+m,:);
+    end
+    ## The filters shifted to the right are symmetric with those to the left.
+    F(k+2:n,:) = (-1)^m*F(k:-1:1,n:-1:1);
+
+  endif
+  F =  F * ( prod(1:m) / (ts^m) );
+endfunction
+
+%!test
+%! N=2^12;
+%! t=[0:N-1]'/N;
+%! dt=t(2)-t(1);
+%! w = 2*pi*50;
+%! offset = 0.5; # 50 Hz carrier
+%! # exponential modulation and its derivatives
+%! d = 1+exp(-3*(t-offset));
+%! dd = -3*exp(-3*(t-offset));
+%! d2d = 9*exp(-3*(t-offset));
+%! d3d = -27*exp(-3*(t-offset));
+%! # modulated carrier and its derivatives
+%! x = d.*sin(w*t);
+%! dx = dd.*sin(w*t) + w*d.*cos(w*t);
+%! d2x = (d2d-w^2*d).*sin(w*t) + 2*w*dd.*cos(w*t);
+%! d3x = (d3d-3*w^2*dd).*sin(w*t) + (3*w*d2d-w^3*d).*cos(w*t);
+%!
+%! y = sgolayfilt(x,sgolay(8,41,0,dt));
+%! assert(norm(y-x)/norm(x),0,5e-6);
+%!
+%! y = sgolayfilt(x,sgolay(8,41,1,dt));
+%! assert(norm(y-dx)/norm(dx),0,5e-6);
+%! 
+%! y = sgolayfilt(x,sgolay(8,41,2,dt));
+%! assert(norm(y-d2x)/norm(d2x),0,1e-5);
+%! 
+%! y = sgolayfilt(x,sgolay(8,41,3,dt));
+%! assert(norm(y-d3x)/norm(d3x),0,1e-4);