From: <jpi...@us...> - 2012-05-30 11:41:46
|
Revision: 10541 http://octave.svn.sourceforge.net/octave/?rev=10541&view=rev Author: jpicarbajal Date: 2012-05-30 11:41:35 +0000 (Wed, 30 May 2012) Log Message: ----------- system-identification: power spectrum tool. Wrapper Added Paths: ----------- trunk/octave-forge/main/system-identification/inst/tisean/pspec_bin.m trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m Added: trunk/octave-forge/main/system-identification/inst/tisean/pspec_bin.m =================================================================== --- trunk/octave-forge/main/system-identification/inst/tisean/pspec_bin.m (rev 0) +++ trunk/octave-forge/main/system-identification/inst/tisean/pspec_bin.m 2012-05-30 11:41:35 UTC (rev 10541) @@ -0,0 +1,70 @@ +%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...> +%% +%% 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 +%% 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} {@var{S} = } pspec_bin (@var{data}) +%% Computes a power spectrum by binning adjacent frequencies. This function @ +%% calls @code{spectrum} from the TISEAN package. +%% +%% @var{S} contains the power spectrum, @var{f} the frequency vector. +%% +%% @seealso{pspec_mep} +%% +%% @end deftypefn + +function s = pspec_bin (data, varargin) + + data = data(:); + [nT n] = size (data); + if n !=1 + warning ('Tisian:InvalidInputFormat', ... + 'Only first column of data will be used'); + end + + # --- Parse arguments --- # + parser = inputParser (); + parser.FunctionName = "pspec_bin"; + parser = addParamValue (parser,'Res', 1/nT, @(x)x>0); + parser = addParamValue (parser,'Sampling', 1, @(x)x>0); + parser = parse(parser,varargin{:}); + + res = parser.Results.Res; + sampl = parser.Results.Sampling; + + clear parser + # ------ # + + + flag.w = sprintf("-w%d", res); + flag.f = sprintf("-f%d", sampl); + + infile = tmpnam (); + outfile = tmpnam (); + + %% Write data to file + save ('-ascii',infile, 'data'); + + %% Prepare format of system call + func = file_in_loadpath ("spectrum"); + syscmd = sprintf("%s %s %s -o%s -V0 %s", ... + func, flag.f, flag.w, outfile, infile); + + %% Function call + system (syscmd); + + s = load (outfile); + s = complex (s(:,1), s(:,2)); + +endfunction Added: trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m =================================================================== --- trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m (rev 0) +++ trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m 2012-05-30 11:41:35 UTC (rev 10541) @@ -0,0 +1,104 @@ +%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...> +%% +%% 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 +%% 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} {[@var{S} @var{f} @var{AR} @var{ARerr}] = } pspec_mep (@var{data}) +%% @deftypefnx {Function File} {[@dots{}] = } pspec_mep (@dosts{},@var{property},@var{value}) +%% Estimates the power spectrum of a scalar data set on the basis of the maximum @ +%% entropy principle. This function calls @code{mem_spec} from the TISEAN package. +%% +%% @var{S} contains the power spectrum, @var{f} the frequency vector. +%% +%% Optional arguments are: +%% @table @samp +%% @item Poles +%% Number of poles of the AR model. Deafault 128. +%% @item Freqs +%% Number of frequencies to use. Defualt 2000 Hz. +%% @item Sampling +%% Sampling frequency of the data. Default 1 Hz +%% @item AR +%% This is just a swtich, it doesn't need to be followed by a value. If present +%% the coefficients of the AR model are returned in @var{AR} and the estimation +%% error in @var{ARerr}. +%% @end table +%% +%% @seealso{pspec_bin} +%% +%% @end deftypefn + +function [s f AR ARerr] = pspec_mep (data, varargin) + + # --- Parse arguments --- # + parser = inputParser (); + parser.FunctionName = "pspec_mep"; + parser = addParamValue (parser,'Poles', 128, @(x)x>0); + parser = addParamValue (parser,'Freqs', 2e3, @(x)x>0); + parser = addParamValue (parser,'Sampling', 1, @(x)x>0); + parser = addSwitch (parser,'AR'); + parser = parse(parser,varargin{:}); + + poles = parser.Results.Poles; + freqs = parser.Results.Freqs; + sampl = parser.Results.Sampling; + ar = parser.Results.AR; + + clear parser + # ------ # + + data = data(:); + [nT n] = size (data); + if n !=1 + warning ('Tisian:InvalidInputFormat', ... + 'Only first column of data will be used'); + end + + flag.p = sprintf("-p%d", poles); + flag.P = sprintf("-P%d", freqs); + flag.f = sprintf("-f%d", sampl); + flag.V = sprintf("-V%d", 2*ar); + + infile = tmpnam (); + outfile = tmpnam (); + + %% Write data to file + save ('-ascii',infile, 'data'); + + %% Prepare format of system call + func = file_in_loadpath ("mem_spec"); + syscmd = sprintf("%s %s %s %s %s -o%s %s", ... + func, flag.p, flag.P, flag.f, flag.V, outfile, infile); + + %% Function call + system (syscmd); + + %% Read results + if ar + ARerrtxt = textread (outfile, "%s", 1){1}; + [~, AR] = textread (outfile, "%s%f", poles, "headerlines", 1); + [f,s] = textread (outfile, "%f%f", freqs, "headerlines", 1+poles); + + %% Parse AR error + ARerr = str2num (strsplit (ARerrtxt, "="){2}); + else + s = load (outfile); + s = complex (s(:,1), s(:,2)); + AR = []; + ARerr = []; + end + + + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <jpi...@us...> - 2012-05-30 18:38:00
|
Revision: 10543 http://octave.svn.sourceforge.net/octave/?rev=10543&view=rev Author: jpicarbajal Date: 2012-05-30 18:37:53 +0000 (Wed, 30 May 2012) Log Message: ----------- system-identification: nonlinear prediction error. Tisean Wrapper Modified Paths: -------------- trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m Added Paths: ----------- trunk/octave-forge/main/system-identification/inst/tisean/nstatz.m trunk/octave-forge/main/system-identification/inst/tisean/recurrplt.m Added: trunk/octave-forge/main/system-identification/inst/tisean/nstatz.m =================================================================== --- trunk/octave-forge/main/system-identification/inst/tisean/nstatz.m (rev 0) +++ trunk/octave-forge/main/system-identification/inst/tisean/nstatz.m 2012-05-30 18:37:53 UTC (rev 10543) @@ -0,0 +1,109 @@ +%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...> +%% +%% 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 +%% 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} {@var{m} = } nstatz (@var{data}, @var{n}) +%% @deftypefn {Function File} {[@var{f} @var{p} @var{ferr}] = } nstatz (@dots{}) +%% @deftypefnx {Function File} {@dots{} = } nststz (@dots{},@var{property},@var{value}) +%% Seeks for nonstationarity in a time series by dividing it into a number of @ +%% segments and calculating the cross-forecast errors between the different @ +%% segments. The model used for the forecast is zeroth order model as proposed by @ +%% Schreiber. This function calls @code{nstat_z} from the TISEAN package. +%% +%% @end deftypefn + +function varargout = nstatz (data, n, varargin) + + data = data(:); + [nT nc] = size (data); + if nc !=1 + warning ('Tisian:InvalidInputFormat', ... + 'Only first column of data will be used'); + end + amp = abs((max(data)-min(data))); + + # --- Parse arguments --- # + parser = inputParser (); + parser.FunctionName = "nstatz"; + parser = addParamValue (parser,'Edim', 3, @(x)x>0); + parser = addParamValue (parser,'Delay', 1, @(x)x>0); + parser = addParamValue (parser,'Points', Inf, @(x)x>0); + parser = addParamValue (parser,'minNeigh', 30, @(x)x>0); + parser = addParamValue (parser,'Rini', amp*1e-3, @(x)x>0); + parser = addParamValue (parser,'Rslope', 1.2, @(x)x>0); + parser = addParamValue (parser,'Fstep', 1, @(x)x>0); + parser = addParamValue (parser,'CausW', 1, @(x)x>0); + parser = parse(parser,varargin{:}); + + edim = parser.Results.Edim; + delay = parser.Results.Delay; + points = parser.Results.Points; + minN = parser.Results.minNeigh; + rini = parser.Results.Rini; + rslope = parser.Results.Rslope; + fstep = parser.Results.Fstep; + causw = parser.Results.CausW; + if ismember("CausW", parser.UsingDefaults) + causw = fstep; + end + clear parser + # ------ # + + flag.num = sprintf("-#%d", n); + + flag.m = sprintf("-m%d", edim); + flag.d = sprintf("-d%d", delay); + if isinf (points) + flag.n = ""; + else + flag.n = sprintf("-d%d", points); + end + flag.k = sprintf("-k%d", minN); + flag.r = sprintf("-r%f", rini); + flag.f = sprintf("-f%f", rslope); + flag.s = sprintf("-s%d", fstep); + flag.C = sprintf("-C%d", causw); + + infile = tmpnam (); + outfile = tmpnam (); + + %% Write data to file + save ('-ascii',infile, 'data'); + + %% Prepare format of system call + func = file_in_loadpath ("nstat_z"); + syscmd = sprintf("%s %s %s %s %s %s %s %s %s %s -o%s -V0 %s", func, flag.num, ... + flag.m, flag.d, flag.n, flag.k, flag.r, flag.f, flag.s, flag.C, ... + outfile, infile); + + %% Function call + system (syscmd); + s = load (outfile); + f = s(:,1); + p = s(:,2); + ferr = s(:,3); + + if nargout == 1 + m = NA (n); + idx = sub2ind ([n,n], f, p); + m(idx) = ferr; + varargout{1} = m; + elseif nargout == 3 + varargout = {f,p,ferr}; + else + print_usage (); + end + +endfunction Modified: trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m =================================================================== --- trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m 2012-05-30 12:03:18 UTC (rev 10542) +++ trunk/octave-forge/main/system-identification/inst/tisean/pspec_mep.m 2012-05-30 18:37:53 UTC (rev 10543) @@ -15,7 +15,7 @@ %% -*- texinfo -*- %% @deftypefn {Function File} {[@var{S} @var{f} @var{AR} @var{ARerr}] = } pspec_mep (@var{data}) -%% @deftypefnx {Function File} {[@dots{}] = } pspec_mep (@dosts{},@var{property},@var{value}) +%% @deftypefnx {Function File} {[@dots{}] = } pspec_mep (@dots{},@var{property},@var{value}) %% Estimates the power spectrum of a scalar data set on the basis of the maximum @ %% entropy principle. This function calls @code{mem_spec} from the TISEAN package. %% Added: trunk/octave-forge/main/system-identification/inst/tisean/recurrplt.m =================================================================== --- trunk/octave-forge/main/system-identification/inst/tisean/recurrplt.m (rev 0) +++ trunk/octave-forge/main/system-identification/inst/tisean/recurrplt.m 2012-05-30 18:37:53 UTC (rev 10543) @@ -0,0 +1,91 @@ +%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...> +%% +%% 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 +%% 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} {[@var{M} @var{i} @var{j}] = } recurrplt (@var{data}) +%% @deftypefnx {Function File} {[@dots{}] = } recurrplt (@dots{},@var{property},@var{value}) +%% Produces a recurrence plot of the, possibly multivariate, data set. That +%% means, for each point in the data set it looks for all points, such that the +%% distance between these two points is smaller than a given size in a given +%% embedding space. This function calls @code{recurr} from the TISEAN package. +%% +%% Returns a pairs of integers representing the indexes of the pairs of points +%% having a distance smaller than R. +%% +%% Optional arguments are: +%% @table @samp +%% @item Edim +%% A 1x2 vector with number of components, embedding dimension. Deafult [1,2]. +%% @item Delay +%% Time delay. Default 1. +%% @item R +%% Size of the neighbourhood. Default data-interval/1000. +%% @item Decimate +%% Number between 0-1. Output only the percentage a of all points found +%% should help to keep the output size reasonably smal +%% @end table +%% +%% @seealso{delayvec} +%% +%% @end deftypefn + +function [i j] = recurrplt (data, varargin) + + [nT nc] = size (data); + amp = abs((max(data(:))-min(data(:)))); + + # --- Parse arguments --- # + parser = inputParser (); + parser.FunctionName = "recurrplt"; + parser = addParamValue (parser,'Edim', [1,2], @ismatrix); + parser = addParamValue (parser,'Delay', 1, @(x)x>0); + parser = addParamValue (parser,'R', amp*1e-3, @(x)x>0); + parser = addParamValue (parser,'Decimate', 1, @(x)x>0 && x <=1); + parser = parse(parser,varargin{:}); + + edim = parser.Results.Edim; + delay = parser.Results.Delay; + radius = parser.Results.R; + percent = parser.Results.Decimate; + + clear parser + # ------ # + + flag.m = sprintf("-m%d,%d", edim); + flag.d = sprintf("-d%d", delay); + flag.r = sprintf("-r%f", radius); + flag.P = sprintf("-%%%f", percent*100); + flag.c = sprintf("-c%d", nc); + + infile = tmpnam (); + outfile = tmpnam (); + + %% Write data to file + save ('-ascii',infile, 'data'); + + %% Prepare format of system call + func = file_in_loadpath ("recurr"); + syscmd = sprintf("%s %s %s %s %s %s -V0 -o%s %s", ... + func, flag.c, flag.m, flag.d, flag.r, flag.P, outfile, infile); + + %% Function call + system (syscmd); + + %% Read results + d = load(outfile); + i = d(:,1); + j = d(:,2); + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |