From bug-request at octave dot org Thu Dec 16 07:32:18 2004 Subject: [octave-forge] two patches for Savintsky-Golay filters From: "Pascal A. Dupuis" To: bug at octave dot org Date: Thu, 16 Dec 2004 05:01:26 -0600 --k1lZvvs/B4yU6o8G Content-Type: text/plain; charset=us-ascii Content-Disposition: inline Hello, I wandered why the sgolay and sgolayfilt couldn't be used also to estimate time derivatives. The reason is quite simple: the filter() function is a convolution, so the filters coefficient are used in the reverse order. I made a small change to sgolayfilt in such a way that the SG coeffs are defined as in Numerical Recipes, but are used the right way in the computation. Now the coefficients returned by sgolay for m>0 can be used to compute the mth derivative. Best regards Pascal -- Dr. ir. Pascal Dupuis K. U. Leuven, ESAT/ELECTA (formerly ELEN): http://www.esat.kuleuven.ac.be/ Kasteelpark Arenberg, 10; B-3001 Leuven-Heverlee, Belgium Tel. +32-16-32 10 21 -- Fax +32-16-32 19 85 --k1lZvvs/B4yU6o8G Content-Type: image/x-coreldrawpattern Content-Disposition: attachment; filename="sgolayfilt.pat" Content-Transfer-Encoding: quoted-printable --- sgolayfilt.m.orig 2004-12-15 17:39:37.000000000 +0100=0A+++ sgolayfilt.= m 2004-12-15 17:54:27.000000000 +0100=0A at @ -14,7 +14,7 @@=0A ## along with = this program; if not, write to the Free Software=0A ## Foundation, Inc., 59= Temple Place, Suite 330, Boston, MA 02111-1307 USA=0A =0A-## y =3D sgola= yfilt (x, p, n)=0A+## y =3D sgolayfilt (x, p, n [, m [, ts]])=0A ## Smoo= th the data in x with a Savitsky-Golay smoothing filter of =0A ## polyno= mial order p and length n, n odd, n > p. By default, p=3D3=0A ## and n= =3Dp+2 or n=3Dp+3 if p is even dot =0A at @ -36,20 +36,24 @@=0A ##=0A ## See also:= sgolay=0A =0A-## TODO: Doesn't accept "weight vector" as fourth parameter.= =0A+## 15 Dec 2004 modified by Pascal Dupuis =0A+## Author: Paul Kienzle =0A+=0A #= # TODO: Patch filter.cc so that it accepts matrix arguments=0A =0A-function= y =3D sgolayfilt (x, p, n)=0A+function y =3D sgolayfilt (x, p, n, m, ts)= =0A =0A if nargin < 1 || nargin > 3 =0A- usage("y =3D sgolayfilt(x,p,n= ) or y =3D sgolayfilt(x,F)"); =0A+ usage("y =3D sgolayfilt(x,p,n [, m [,= ts]]) or y =3D sgolayfilt(x,F)"); =0A endif=0A =0A if (nargin < 2)=0A = p =3D 3;=0A endif=0A- if (nargin =3D=3D 3)=0A- F =3D sgolay(p, n)= ;=0A+ if nargin < 4, m =3D 0; endif=0A+ if nargin < 5, ts =3D 1; endif=0A= + if (nargin >=3D 3)=0A+ F =3D sgolay(p, n, m, ts);=0A elseif (prod(s= ize(p)) =3D=3D 1)=0A n =3D p+3-rem(p,2);=0A F =3D sgolay(p, n);=0A at = at -73,8 +77,12 @@=0A ## The last k rows of F are used to filter the last = k points=0A ## of the data set based on the last n points of the dataset.= =0A ## The remaining data is filtered using the central row of F.=0A+ ##= As the filter coefficients are used in the reverse order of what=0A+ ## s= eems the logical notation, reverse F(k+1, :) so that antisymmetric=0A+ ## = sequences are used with the right sign.=0A+=0A k =3D floor(n/2);=0A- z = =3D filter(F(k+1,:), 1, x);=0A+ z =3D filter(F(k+1,n:-1:1), 1, x);=0A y = =3D [ F(1:k,:)*x(1:n,:) ; z(n:len,:) ; F(k+2:n,:)*x(len-n+1:len,:) ];=0A = =0A if (transpose) y =3D y.'; endif=0A--- sgolay.m.orig 2004-12-15 17:39:= 24.000000000 +0100=0A+++ sgolay.m 2004-12-15 17:56:15.000000000 +0100=0A at @ = -14,9 +14,10 at @=0A ## along with this program; if not, write to the Free So= ftware=0A ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 0211= 1-1307 USA=0A =0A-## F =3D sgolay (p, n)=0A+## F =3D sgolay (p, n [, m [, = ts]])=0A ## Computes the filter coefficients for all Savitzsky-Golay smoo= thing=0A-## filters of order p for length n (odd).=0A+## filters of ord= er p for length n (odd). m can be used in order to=0A+## get directly the= mth derivative. In this case, ts is a scaling factor. =0A ##=0A ## The ear= ly rows of F smooth based on future values and later rows=0A ## smooth base= d on past values, with the middle row using half future=0A at @ -36,30 +37,29 = at @=0A ##=0A ## See also: sgolayfilt=0A =0A+## 15 Dec 2004 modified by Pasca= l Dupuis =0A ## Author: Paul Kienzle =0A ## Based on smooth.m by E. Farhi =0A =0A-## TODO: Doesn't accept "weight vector" as fourth= parameter since=0A-## TODO: Should be able to estimate derivatives using s= econd and=0A-## TODO: subsequent columns of A, but they seem to have the wr= ong=0A-## TODO: sign and the wrong scale, so I won't put that in.=0A+functi= on F =3D sgolay (p, n, m, ts)=0A =0A-function F =3D sgolay (p, n)=0A-=0A- = if (nargin < 2 || nargin > 3)=0A- usage ("F =3D sgolay (p, n)");=0A+ if= (nargin < 2 || nargin > 4)=0A+ usage ("F =3D sgolay (p, n [, m [, ts]])= ");=0A elseif rem(n,2) !=3D 1=0A error ("sgolay needs an odd filter l= ength n");=0A elseif p >=3D n=0A error ("sgolay needs filter length n= larger than polynomial order p");=0A else =0A+ if nargin < 3, m =3D 0= ; endif=0A+ if nargin < 4, ts =3D 1; endif=0A k =3D floor (n/2);=0A = F =3D zeros (n, n);=0A for row =3D 1:k+1=0A A =3D pinv( ( [(1= :n)-row]'*ones(1,p+1) ) .^ ( ones(n,1)*[0:p] ) );=0A- F(row,:) =3D A(1= ,:);=0A+ F(row,:) =3D A(1+m,:);=0A end=0A F(k+2:n,:) =3D F(k:-= 1:1,n:-1:1);=0A endif=0A-=0A+ if m > 1, F =3D F * factorial(m); endif= =0A+ if m > 0, F =3D F / (ts^m); endif=0A endfunction=0A --k1lZvvs/B4yU6o8G-- ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------