Mercurial > forge
changeset 10299:4f08f44a2705 octave-forge
Added cubicwgt and the wavelet coefficient functions.
author | benjf5 |
---|---|
date | Tue, 22 May 2012 02:02:32 +0000 |
parents | a187606d8d02 |
children | 8ffb34d12d44 |
files | extra/lssa/cubicwgt.m extra/lssa/nucorrcoeff.m extra/lssa/nuwaveletcoeff.m |
diffstat | 3 files changed, 87 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/lssa/cubicwgt.m Tue May 22 02:02:32 2012 +0000 @@ -0,0 +1,22 @@ +## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com> +## +## 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/>. + +## This function implements the windowing function on page 10 of the doc. +## if t is in [-1,1] then the windowed term is a = 1 + ( |t|^2 * ( 2|t| - 3 ) +## else the windowed term is 0. +function a = cubicwgt(t) ## where t is the set of time values + a = abs(t); + a = ifelse( ( a < 1 ), 1 + ( ( a .^ 2 ) .* ( 2 .* a - 3 ) ), a = 0); +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/lssa/nucorrcoeff.m Tue May 22 02:02:32 2012 +0000 @@ -0,0 +1,38 @@ +## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com> +## +## 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/>. + +## 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; + ## 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? + length(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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/lssa/nuwaveletcoeff.m Tue May 22 02:02:32 2012 +0000 @@ -0,0 +1,27 @@ +## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com> +## +## 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/>. + +function coeff = nuwaveletcoeff( X , Y , t , o , wgt = cubicwgt , wgtrad = 1 ) + so = 0.05 .* o; + 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