You can subscribe to this list here.
2001 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}
(10) 
_{Aug}
(5) 
_{Sep}
(3) 
_{Oct}
(41) 
_{Nov}
(41) 
_{Dec}
(33) 

2002 
_{Jan}
(75) 
_{Feb}
(10) 
_{Mar}
(170) 
_{Apr}
(174) 
_{May}
(66) 
_{Jun}
(11) 
_{Jul}
(10) 
_{Aug}
(44) 
_{Sep}
(73) 
_{Oct}
(28) 
_{Nov}
(139) 
_{Dec}
(52) 
2003 
_{Jan}
(35) 
_{Feb}
(93) 
_{Mar}
(62) 
_{Apr}
(10) 
_{May}
(55) 
_{Jun}
(70) 
_{Jul}
(37) 
_{Aug}
(16) 
_{Sep}
(56) 
_{Oct}
(31) 
_{Nov}
(57) 
_{Dec}
(83) 
2004 
_{Jan}
(85) 
_{Feb}
(67) 
_{Mar}
(27) 
_{Apr}
(37) 
_{May}
(75) 
_{Jun}
(85) 
_{Jul}
(160) 
_{Aug}
(68) 
_{Sep}
(104) 
_{Oct}
(25) 
_{Nov}
(39) 
_{Dec}
(23) 
2005 
_{Jan}
(10) 
_{Feb}
(45) 
_{Mar}
(43) 
_{Apr}
(19) 
_{May}
(108) 
_{Jun}
(31) 
_{Jul}
(41) 
_{Aug}
(23) 
_{Sep}
(65) 
_{Oct}
(58) 
_{Nov}
(44) 
_{Dec}
(54) 
2006 
_{Jan}
(96) 
_{Feb}
(27) 
_{Mar}
(69) 
_{Apr}
(59) 
_{May}
(67) 
_{Jun}
(35) 
_{Jul}
(13) 
_{Aug}
(461) 
_{Sep}
(160) 
_{Oct}
(399) 
_{Nov}
(32) 
_{Dec}
(72) 
2007 
_{Jan}
(316) 
_{Feb}
(305) 
_{Mar}
(318) 
_{Apr}
(54) 
_{May}
(194) 
_{Jun}
(173) 
_{Jul}
(282) 
_{Aug}
(91) 
_{Sep}
(227) 
_{Oct}
(365) 
_{Nov}
(168) 
_{Dec}
(18) 
2008 
_{Jan}
(71) 
_{Feb}
(111) 
_{Mar}
(155) 
_{Apr}
(173) 
_{May}
(70) 
_{Jun}
(67) 
_{Jul}
(55) 
_{Aug}
(83) 
_{Sep}
(32) 
_{Oct}
(68) 
_{Nov}
(80) 
_{Dec}
(29) 
2009 
_{Jan}
(46) 
_{Feb}
(18) 
_{Mar}
(95) 
_{Apr}
(76) 
_{May}
(140) 
_{Jun}
(98) 
_{Jul}
(84) 
_{Aug}
(123) 
_{Sep}
(94) 
_{Oct}
(131) 
_{Nov}
(142) 
_{Dec}
(125) 
2010 
_{Jan}
(128) 
_{Feb}
(158) 
_{Mar}
(172) 
_{Apr}
(134) 
_{May}
(94) 
_{Jun}
(84) 
_{Jul}
(32) 
_{Aug}
(127) 
_{Sep}
(167) 
_{Oct}
(109) 
_{Nov}
(69) 
_{Dec}
(78) 
2011 
_{Jan}
(39) 
_{Feb}
(58) 
_{Mar}
(52) 
_{Apr}
(47) 
_{May}
(56) 
_{Jun}
(76) 
_{Jul}
(55) 
_{Aug}
(54) 
_{Sep}
(165) 
_{Oct}
(255) 
_{Nov}
(328) 
_{Dec}
(263) 
2012 
_{Jan}
(82) 
_{Feb}
(147) 
_{Mar}
(400) 
_{Apr}
(216) 
_{May}
(209) 
_{Jun}
(160) 
_{Jul}
(86) 
_{Aug}
(141) 
_{Sep}
(156) 
_{Oct}
(6) 
_{Nov}

_{Dec}

S  M  T  W  T  F  S 



1
(12) 
2
(1) 
3
(4) 
4
(4) 
5
(4) 
6
(5) 
7
(5) 
8
(9) 
9
(9) 
10
(11) 
11
(11) 
12
(14) 
13

14
(10) 
15
(8) 
16
(5) 
17
(1) 
18
(10) 
19
(2) 
20
(5) 
21
(14) 
22
(15) 
23
(9) 
24
(12) 
25
(2) 
26
(1) 
27
(5) 
28
(3) 
29
(3) 
30
(10) 
31
(5) 


From: <paramaniac@us...>  20120525 23:14:20

Revision: 10522 http://octave.svn.sourceforge.net/octave/?rev=10522&view=rev Author: paramaniac Date: 20120525 23:14:13 +0000 (Fri, 25 May 2012) Log Message:  controldevel: save some early draft code Added Paths:  trunk/octaveforge/extra/controldevel/devel/rarx.m Added: trunk/octaveforge/extra/controldevel/devel/rarx.m ===================================================================  trunk/octaveforge/extra/controldevel/devel/rarx.m (rev 0) +++ trunk/octaveforge/extra/controldevel/devel/rarx.m 20120525 23:14:13 UTC (rev 10522) @@ 0,0 +1,202 @@ +## Copyright (C) 2012 Lukas F. Reichlin +## +## This file is part of LTI Syncope. +## +## LTI Syncope 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. +## +## LTI Syncope 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>;. + +## * texinfo * +## @deftypefn {Function File} {@var{sys} =} arx (@var{dat}, @var{na}, @var{nb}) +## ARX +## @end deftypefn + +## Author: Lukas Reichlin <lukas.reichlin@...> +## Created: April 2012 +## Version: 0.1 + +function [sys, varargout] = rarx (dat, na, nb) + + ## TODO: delays + + if (nargin != 3) + print_usage (); + endif + + if (! isa (dat, "iddata")) + error ("arx: first argument must be an iddata dataset"); + endif + + ## p: outputs, m: inputs, ex: experiments + [~, p, m, ex] = size (dat); + + ## extract data + Y = dat.y; + U = dat.u; + tsam = dat.tsam; + + ## multiexperiment data requires equal sampling times + if (ex > 1 && ! isequal (tsam{:})) + error ("arx: require equally sampled experiments"); + else + tsam = tsam{1}; + endif + + + if (is_real_scalar (na, nb)) + na = repmat (na, p, 1); # na(pby1) + nb = repmat (nb, p, m); # nb(pbym) + elseif (! (is_real_vector (na) && is_real_matrix (nb) \ + && rows (na) == p && rows (nb) == p && columns (nb) == m)) + error ("arx: require na(%dx1) instead of (%dx%d) and nb(%dx%d) instead of (%dx%d)", \ + p, rows (na), columns (na), p, m, rows (nb), columns (nb)); + endif + + max_nb = max (nb, [], 2); # one maximum for each row/output, max_nb(pby1) + n = max (na, max_nb); # n(pby1) + + ## create empty cells for numerator and denominator polynomials + num = cell (p, m+p); + den = cell (p, m+p); + + ## MIMO (pbym) models are identified as p MISO (1bym) models + ## For multiexperiment data, minimize the trace of the error + for i = 1 : p # for every output + Phi = cell (ex, 1); # one regression matrix per experiment + for e = 1 : ex # for every experiment + ## avoid warning: toeplitz: column wins antidiagonal conflict + ## therefore set first row element equal to y(1) + PhiY = toeplitz (Y{e}(1:end1, i), [Y{e}(1, i); zeros(na(i)1, 1)]); + ## create MISO Phi for every experiment + PhiU = arrayfun (@(x) toeplitz (U{e}(1:end1, x), [U{e}(1, x); zeros(nb(i,x)1, 1)]), 1:m, "uniformoutput", false); + Phi{e} = (horzcat (PhiY, PhiU{:}))(n(i):end, :); + endfor + + ## compute parameter vector Theta + Theta = __theta__ (Phi, Y, i, n); + + ## extract polynomial matrices A and B from Theta + ## A is a scalar polynomial for output i, i=1:p + ## B is polynomial row vector (1bym) for output i + A = [1; Theta(1:na(i))]; # a0 = 1, a1 = Theta(1), an = Theta(n) + ThetaB = Theta(na(i)+1:end); # all polynomials from B are in one column vector + B = mat2cell (ThetaB, nb(i,:)); # now separate the polynomials, one for each input + B = reshape (B, 1, []); # make B a row cell (1bym) + B = cellfun (@(x) [0; x], B, "uniformoutput", false); # b0 = 0 (leading zero required by filt) + + ## add error inputs + Be = repmat ({0}, 1, p); # there are as many error inputs as system outputs (p) + Be(i) = 1; # inputs m+1:m+p are zero, except m+i which is one + num(i, :) = [B, Be]; # numerator polynomials for output i, individual for each input + den(i, :) = repmat ({A}, 1, m+p); # in a row (output i), all inputs have the same denominator polynomial + endfor + + ## A(q) y(t) = B(q) u(t) + e(t) + ## there is only one A per row + ## B(z) and A(z) are a Matrix Fraction Description (MFD) + ## y = A^1(q) B(q) u(t) + A^1(q) e(t) + ## since A(q) is a diagonal polynomial matrix, its inverse is trivial: + ## the corresponding transfer function has common row denominators. + + sys = filt (num, den, tsam); # filt creates a transfer function in z^1 + + ## compute initial state vector x0 if requested + ## this makes only sense for statespace models, therefore convert TF to SS + if (nargout > 1) + sys = prescale (ss (sys(:,1:m))); + x0 = slib01cd (Y, U, sys.a, sys.b, sys.c, sys.d, 0.0) + ## return x0 as vector for singleexperiment data + ## instead of a cell containing one vector + if (numel (x0) == 1) + x0 = x0{1}; + endif + varargout{1} = x0; + endif + +endfunction + + +%function theta = __theta__ (phi, y, i, n) +function Theta = __theta__ (Phi, Y, i, n) + + + if (numel (Phi) == 1) # singleexperiment dataset + % recursive leastsquares with efficient matrix inversion + lambda = default 1 + Theta = ? + P = ? + for t = + phi = Phi{1}(t,:); + y = Y{1}(t?offset n(i), :); + ## note phi != phi.' + den = lambda + phi*P*phi.'; + L = P * phi / den; + P = (P  (P * phi.' * phi * P) / den) / lambda + Theta += L * (y  phi*Theta); + endfor +%{ + ## use "squareroot algorithm" + A = horzcat (phi{1}, y{1}(n(i)+1:end, i)); # [Phi, Y] + R0 = triu (qr (A, 0)); # 0 for economysize R (without zero rows) + R1 = R0(1:end1, 1:end1); # R1 is triangular  can we exploit this in R1\R2? + R2 = R0(1:end1, end); + theta = __ls_svd__ (R1, R2); # R1 \ R2 + + ## Theta = Phi \ Y(n+1:end, :); # naive formula + ## theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i)); +%} + else # multiexperiment dataset + ## TODO: find more sophisticated formula than + ## Theta = (Phi1' Phi + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...) + + ## covariance matrix C = (Phi1' Phi + Phi2' Phi2 + ...) + tmp = cellfun (@(Phi) Phi.' * Phi, phi, "uniformoutput", false); + rc = cellfun (@rcond, tmp); # C auch noch testen? QR oder SVD? + C = plus (tmp{:}); + + ## PhiTY = (Phi1' Y1 + Phi2' Y2 + ...) + tmp = cellfun (@(Phi, Y) Phi.' * Y(n(i)+1:end, i), phi, y, "uniformoutput", false); + PhiTY = plus (tmp{:}); + + ## pseudoinverse Theta = C \ Phi'Y + theta = __ls_svd__ (C, PhiTY); + endif + +endfunction + + +function x = __ls_svd__ (A, b) + + ## solve the problem Ax=b + ## x = A\b would also work, + ## but this way we have better control and warnings + + ## solve linear least squares problem by pseudoinverse + ## the pseudoinverse is computed by singular value decomposition + ## M = U S V* > M+ = V S+ U* + ## Th = Ph \ Y = Ph+ Y + ## Th = V S+ U* Y, S+ = 1 ./ diag (S) + + [U, S, V] = svd (A, 0); # 0 for "economy size" decomposition + S = diag (S); # extract main diagonal + r = sum (S > eps*S(1)); + if (r < length (S)) + warning ("arx: rankdeficient coefficient matrix"); + warning ("sampling time too small"); + warning ("persistence of excitation"); + endif + V = V(:, 1:r); + S = S(1:r); + U = U(:, 1:r); + x = V * (S .\ (U' * b)); # U' is the conjugate transpose + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <benjf5@us...>  20120525 14:49:59

Revision: 10521 http://octave.svn.sourceforge.net/octave/?rev=10521&view=rev Author: benjf5 Date: 20120525 14:49:45 +0000 (Fri, 25 May 2012) Log Message:  Fixed lssa filenames, blatant error in lsreal. Added Paths:  trunk/octaveforge/extra/lssa/lscorrcoeff.m trunk/octaveforge/extra/lssa/lsreal.m trunk/octaveforge/extra/lssa/lswaveletcoeff.m Removed Paths:  trunk/octaveforge/extra/lssa/nucorrcoeff.m trunk/octaveforge/extra/lssa/nureal.m trunk/octaveforge/extra/lssa/nuwaveletcoeff.m Copied: trunk/octaveforge/extra/lssa/lscorrcoeff.m (from rev 10520, trunk/octaveforge/extra/lssa/nucorrcoeff.m) ===================================================================  trunk/octaveforge/extra/lssa/lscorrcoeff.m (rev 0) +++ trunk/octaveforge/extra/lssa/lscorrcoeff.m 20120525 14:49:45 UTC (rev 10521) @@ 0,0 +1,67 @@ +## Copyright (C) 2012 Benjamin Lewis <benjf5@...> +## +## 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/>;. + +## * texinfo * +## @deftypefn {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq) +## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window) +## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window, winradius) +## +## Return the coefficient of the wavelet correlation of time +## series (@var{abc1}, @var{ord1}) and (@var{abc2}, @var{ord2}). +## @var{window} is used to apply a windowing function, its +## default is cubicwgt if left blank, and its radius is 1, +## as defined as the default for @var{winradius}. +## +## @end deftypefn + +## Demo with sin, cos as Nir suggested. +%!demo +%! x = 1:10; +%! y = sin(x); +%! z = cos(x); +%! a = nucorrcoeff(x,y,x,z,0.5,0.9) +%! ## This generates the correlation coefficient at time 0.5 and circular freq. 0.9 + + +## nucorrcoeff, computes a coefficient of the wavelet correlation of two time series + +function coeff = nucorrcoeff(x1, y1, x2, y2, t, o, wgt = @cubicwgt, wgtrad = 1) + so = 0.05 * o; + ## This code can only, as of currently, work on vectors; I haven't figured out a way to make it work on a matrix. + if( ( ndims(x1) == 2 ) && ! ( rows(x1) == 1 ) ) + x1 = reshape(x1,1,length(x1)); + y1 = reshape(y1,1,length(y1)); + x2 = reshape(x2,1,length(x2)); + y2 = reshape(y2,1,length(y2)); + endif + ## The first solution that comes to mind is admittedly slightly ugly and has a data footprint of O(2n) + ## but it is vectorised. + mask = ( abs( x1  t ) * so ) < wgtrad; + mask = mask .* [ 1 : length(mask) ]; + rx1 = x1(mask); ## I've kept the variable names from the R function here + ry1 = y1(mask); ## Needs to have a noisy error if length(y1) != length(x1)  add this! + mask = ( abs( x2  t ) * so ) < wgtrad; + mask = mask .* [ 1 : length(mask) ]; + rx2 = x2(mask); + ry2 = y2(mask); + ## I've used the same mask for all of these as it's an otherwise unimportant variable ... can this leak memory? + lnength(rx1) ##printing this length is probably used as a warning if 0 is returned; I inculded it + ## in particular to maintain an exact duplicate of the R function. + s = sum( wgt( ( rx1  t ) .* so ) ) * sum( wgt( ( rx2  t ) .* so ) ); + coeff = ifelse( s != 0 , ( sum( wgt( ( rx1  t ) .* so ) .* exp( i .* o .* rx1 ) .* ry1 ) + * sum( wgt( ( rx2  t ) .* so ) .* exp( i .* o .* rx2 ) .* ry2 ) ) / s, 0 ); + + +endfunction \ No newline at end of file Copied: trunk/octaveforge/extra/lssa/lsreal.m (from rev 10520, trunk/octaveforge/extra/lssa/nureal.m) ===================================================================  trunk/octaveforge/extra/lssa/lsreal.m (rev 0) +++ trunk/octaveforge/extra/lssa/lsreal.m 20120525 14:49:45 UTC (rev 10521) @@ 0,0 +1,60 @@ +## Copyright (C) 2012 Benjamin Lewis +## +## 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 2 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/>;. + +## * texinfo * +## @deftypefn {Function File} {transform =} nureal ( time, mag, maxfreq, numcoeff, numoctaves) +## +## Return the real leastsquares transform of the time series +## defined, based on the maximal frequency @var{maxfreq}, the +## number of coefficients @var{numcoeff}, and the number of +## octaves @var{numoctaves}. It works only for vectors currently. +## +## @end deftypefn + +function transform = nureal( t, x, omegamax, ncoeff, noctave) + ## the R function runs the following command: + ## nureal( double X, double Y, int min(X,Y), int ncoeff, int noctave, double omegamax, complex rp) + ## this means that in the C, *tptr is X and *xptr is Y. Yes, I know. I think I'll rename them. + ## n is the min of X and Y ... (as is k) and ncoeff is ... ncoeff, while noctave is noctave and + ## o is omegamax. + ## where rp = complex(noctave*ncoeff) so ... I can just store that as noctave*ncoeff and have no + ## problems, I guess. + ## Possibly throw an error if ncoeff <= 0. + k = n = min(length(x),length(t)); ## THIS IS VECTORONLY. I'd need to add another bit of code to + ## make it arraysafe, and that's not knowing right now what else will be necessary. + transform = zeros(1,(noctave * ncoeff)); ## In the C code, this is rendered as a Complex, but Octave doesn't care. + od = 2 ^ (  1 / ncoeff ); ## this will cause a crash if ncoeff=0; prefer error & quit? + o = omegamax; + ## ot is just defined as a Real here, I'm leaving it until it needs to be called. + n1 = 1 / n; ## written n_1 in the C, but I'd prefer to not get into underscores here. + ## zeta and iota are defined as Complex here, leaving them until they need to be defined. + ## I'm not convinced I won't want ncoeff again, so ... + ncoeffp = ncoeff; + ncoeffp *= noctave; + for iter = 1:ncoeffp + ## This method is an application of Eq. 8 on page 6 of the text, as well as Eq. 7 + ot = o .* t; + zeta = sum( ( cos(ot) .* x )  ( sin(ot) .* x .* i ) ); + ot = ot .* 2; + iota = sum( cos(ot)  ( sin(ot) .* i ) ); + zeta = zeta .* n1; + iota = iota .* n1; + transform(iter) = 2 / ( 1  ( real(iota) ^ 2 )  ( imag(iota) ^ 2 ) ) * ( conj(zeta)  (conj(iota) * zeta )); + o = o .* od; + endfor + + ## transform = rp; + +endfunction \ No newline at end of file Copied: trunk/octaveforge/extra/lssa/lswaveletcoeff.m (from rev 10520, trunk/octaveforge/extra/lssa/nuwaveletcoeff.m) ===================================================================  trunk/octaveforge/extra/lssa/lswaveletcoeff.m (rev 0) +++ trunk/octaveforge/extra/lssa/lswaveletcoeff.m 20120525 14:49:45 UTC (rev 10521) @@ 0,0 +1,55 @@ +## Copyright (C) 2012 Benjamin Lewis <benjf5@...> +## +## 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 2 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/>;. + +## * texinfo * +## @deftypefn {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq) +## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window) +## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window, winradius) +## +## Return the coefficient of the wavelet transform of the +## time series (@var{abc}, @var{ord}) at time @var{time} +## and frequency @var{freq}; optional variable @var{window} +## provides a windowing function and defaults to cubicwgt, +## while @var{winradius} is the windowing radius, and defaults +## to 1 (the radius of cubicwgt.) +## +## @end deftypefn + +%!demo +%! x = 1:10; +%! y = sin(x); +%! xt = x'; +%! yt = y'; +%! a = nuwaveletcoeff(x,y,0.5,0.9) +%! b = nuwaveletcoeff(xt,yt,0.5,0.9) +%! ## Generates the wavelet transform coefficient for time 0.5 and circ. freq. 0.9, for row & column vectors. + + +function coeff = nuwaveletcoeff( x , y , t , o , wgt = @cubicwgt , wgtrad = 1 ) + so = 0.05 .* o; + if ( ( ndims(x) == 2 ) && ! ( rows(x) == 1 ) ) + x = reshape(x,1,length(x)); + y = reshape(y,1,length(y)); + endif + mask = abs( x  t ) * so < wgtrad; + mask = mask .* [ 1 : length(mask) ]; + rx = x(mask); + ## This was the fastest way to extract a matching subset that I could think of, but it has a complexity O(2n). + ry = y(mask); + ## Going by the R code, this can use the same mask. + s = sum( wgt( ( x  t ) .* so ) ); + coeff = ifelse( s != 0 , sum( wgt( ( rx  t ) .* so) .* exp( i .* o .* ( rx  t ) ) .* ry ) ./ s , 0 ); + +endfunction \ No newline at end of file Deleted: trunk/octaveforge/extra/lssa/nucorrcoeff.m ===================================================================  trunk/octaveforge/extra/lssa/nucorrcoeff.m 20120524 20:41:47 UTC (rev 10520) +++ trunk/octaveforge/extra/lssa/nucorrcoeff.m 20120525 14:49:45 UTC (rev 10521) @@ 1,67 +0,0 @@ ## Copyright (C) 2012 Benjamin Lewis <benjf5@...> ## ## 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/>;.  ## * texinfo * ## @deftypefn {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq) ## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window) ## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window, winradius) ## ## Return the coefficient of the wavelet correlation of time ## series (@var{abc1}, @var{ord1}) and (@var{abc2}, @var{ord2}). ## @var{window} is used to apply a windowing function, its ## default is cubicwgt if left blank, and its radius is 1, ## as defined as the default for @var{winradius}. ## ## @end deftypefn  ## Demo with sin, cos as Nir suggested. %!demo %! x = 1:10; %! y = sin(x); %! z = cos(x); %! a = nucorrcoeff(x,y,x,z,0.5,0.9) %! ## This generates the correlation coefficient at time 0.5 and circular freq. 0.9   ## nucorrcoeff, computes a coefficient of the wavelet correlation of two time series  function coeff = nucorrcoeff(x1, y1, x2, y2, t, o, wgt = @cubicwgt, wgtrad = 1)  so = 0.05 * o;  ## This code can only, as of currently, work on vectors; I haven't figured out a way to make it work on a matrix.  if( ( ndims(x1) == 2 ) && ! ( rows(x1) == 1 ) )  x1 = reshape(x1,1,length(x1));  y1 = reshape(y1,1,length(y1));  x2 = reshape(x2,1,length(x2));  y2 = reshape(y2,1,length(y2));  endif  ## The first solution that comes to mind is admittedly slightly ugly and has a data footprint of O(2n)  ## but it is vectorised.  mask = ( abs( x1  t ) * so ) < wgtrad;  mask = mask .* [ 1 : length(mask) ];  rx1 = x1(mask); ## I've kept the variable names from the R function here  ry1 = y1(mask); ## Needs to have a noisy error if length(y1) != length(x1)  add this!  mask = ( abs( x2  t ) * so ) < wgtrad;  mask = mask .* [ 1 : length(mask) ];  rx2 = x2(mask);  ry2 = y2(mask);  ## I've used the same mask for all of these as it's an otherwise unimportant variable ... can this leak memory?  lnength(rx1) ##printing this length is probably used as a warning if 0 is returned; I inculded it  ## in particular to maintain an exact duplicate of the R function.  s = sum( wgt( ( rx1  t ) .* so ) ) * sum( wgt( ( rx2  t ) .* so ) );  coeff = ifelse( s != 0 , ( sum( wgt( ( rx1  t ) .* so ) .* exp( i .* o .* rx1 ) .* ry1 )  * sum( wgt( ( rx2  t ) .* so ) .* exp( i .* o .* rx2 ) .* ry2 ) ) / s, 0 );   endfunction \ No newline at end of file Deleted: trunk/octaveforge/extra/lssa/nureal.m ===================================================================  trunk/octaveforge/extra/lssa/nureal.m 20120524 20:41:47 UTC (rev 10520) +++ trunk/octaveforge/extra/lssa/nureal.m 20120525 14:49:45 UTC (rev 10521) @@ 1,72 +0,0 @@ ## Copyright (C) 2012 Benjamin Lewis ## ## 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 2 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/>;.  ## * texinfo * ## @deftypefn {Function File} {transform =} nureal ( mag, time, maxfreq, numcoeff, numoctaves) ## ## Return the real leastsquares transform of the time series ## defined, based on the maximal frequency @var{maxfreq}, the ## number of coefficients @var{numcoeff}, and the number of ## octaves @var{numoctaves}. It works only for vectors currently. ## ## @end deftypefn  function transform = nureal( x, t, omegamax, ncoeff, noctave)  ## the R function runs the following command:  ## nureal( double X, double Y, int min(X,Y), int ncoeff, int noctave, double omegamax, complex rp)  ## this means that in the C, *tptr is X and *xptr is Y. Yes, I know. I think I'll rename them.  ## n is the min of X and Y ... (as is k) and ncoeff is ... ncoeff, while noctave is noctave and  ## o is omegamax.  ## where rp = complex(noctave*ncoeff) so ... I can just store that as noctave*ncoeff and have no  ## problems, I guess.  ## Possibly throw an error if ncoeff <= 0.  k = n = min ( min ( x , t ) ); ## THIS IS VECTORONLY. I'd need to add another bit of code to  ## make it arraysafe, and that's not knowing right now what else will be necessary.  rp = zeros(1,(noctave * ncoeff)); ## In the C code, this is rendered as a Complex, but Octave doesn't care.  od = 2 ^ (  1 / ncoeff ); ## this will cause a crash if ncoeff=0; prefer error & quit?  o = omegamax;  ## ot is just defined as a Real here, I'm leaving it until it needs to be called.  n1 = 1 / n; ## written n_1 in the C, but I'd prefer to not get into underscores here.  ## zeta and iota are defined as Complex here, leaving them until they need to be defined.  ## I'm not convinced I won't want ncoeff again, so ...  ncoeffp = ncoeff;  for ( ncoeffp *= noctave , iter = 1 ; sign(ncoeffp) ; o *= od )  #{  zeta = iota = 0;  for ( SRCFIRST ; SRCAVAIL ; SRCNEXT ) ##This is going to be vectorised shortly.  ## This code can't run yet ... I'm going to work out what SRCFIRST, SRCAVAIL, SRCNEXT are and  ## replace them with what they should be in this context. I've kept them as reminders.  ot = o * SRCT; ## Same with SRCT ... I think it means what was originally set as Y. Macros?  zeta += cos(ot) * SRCX;  zeta = sin(ot) * SRCX * i; ## More sure now. I don't think I can vectorise this ... who am I kidding?  ot *= 2;  iota += cos(ot);  iota = sin(ot) * i;  endfor  }#  ## Commented out the convertedfromC code because it's replaced by the four lines below.  ## This method is an application of Eq. 8 on page 6 of the text, as well as Eq. 7  ot = o .* t;  zeta = sum( ( cos(ot) .* x )  ( sin(ot) .* x .* i ) );  ot = ot .* 2;  iota = sum( cos(ot)  ( sin(ot) .* i ) );  zeta *= n1;  iota *= n1;  rp(iter++) = 2 / ( 1  ( real(iota) ^ 2 )  ( imag(iota) ^ 2 ) ) * ( conj(zeta)  (conj(iota) * zeta ));  endfor   transform = rp;  endfunction \ No newline at end of file Deleted: trunk/octaveforge/extra/lssa/nuwaveletcoeff.m ===================================================================  trunk/octaveforge/extra/lssa/nuwaveletcoeff.m 20120524 20:41:47 UTC (rev 10520) +++ trunk/octaveforge/extra/lssa/nuwaveletcoeff.m 20120525 14:49:45 UTC (rev 10521) @@ 1,55 +0,0 @@ ## Copyright (C) 2012 Benjamin Lewis <benjf5@...> ## ## 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 2 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/>;.  ## * texinfo * ## @deftypefn {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq) ## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window) ## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window, winradius) ## ## Return the coefficient of the wavelet transform of the ## time series (@var{abc}, @var{ord}) at time @var{time} ## and frequency @var{freq}; optional variable @var{window} ## provides a windowing function and defaults to cubicwgt, ## while @var{winradius} is the windowing radius, and defaults ## to 1 (the radius of cubicwgt.) ## ## @end deftypefn  %!demo %! x = 1:10; %! y = sin(x); %! xt = x'; %! yt = y'; %! a = nuwaveletcoeff(x,y,0.5,0.9) %! b = nuwaveletcoeff(xt,yt,0.5,0.9) %! ## Generates the wavelet transform coefficient for time 0.5 and circ. freq. 0.9, for row & column vectors.   function coeff = nuwaveletcoeff( x , y , t , o , wgt = @cubicwgt , wgtrad = 1 )  so = 0.05 .* o;  if ( ( ndims(x) == 2 ) && ! ( rows(x) == 1 ) )  x = reshape(x,1,length(x));  y = reshape(y,1,length(y));  endif  mask = abs( x  t ) * so < wgtrad;  mask = mask .* [ 1 : length(mask) ];  rx = x(mask);  ## This was the fastest way to extract a matching subset that I could think of, but it has a complexity O(2n).  ry = y(mask);  ## Going by the R code, this can use the same mask.  s = sum( wgt( ( x  t ) .* so ) );  coeff = ifelse( s != 0 , sum( wgt( ( rx  t ) .* so) .* exp( i .* o .* ( rx  t ) ) .* ry ) ./ s , 0 );  endfunction \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 