Mercurial > hg > octave-lyh
diff scripts/signal/arma_rnd.m @ 3191:e4f4b2d26ee9
[project @ 1998-10-23 05:43:59 by jwe]
author | jwe |
---|---|
date | Fri, 23 Oct 1998 05:44:01 +0000 |
parents | |
children | f8dde1807dee |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/signal/arma_rnd.m @@ -0,0 +1,79 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## 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, 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 file. If not, write to the Free Software Foundation, +## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +## usage: arma_rnd (a, b, v, t [, n]) +## +## Returns a simulation of the ARMA model +## x(n) = a(1) * x(n-1) + ... + a(k) * x(n-k) + +## + e(n) + b(1) * e(n-1) + ... + b(l) * e(n-l), +## where k is the length of vector a, l is the length of vector b and e +## is gaussian white noise with variance v. The function returns a +## vector of length t. +## +## The optional parameter n gives the number of dummy x(i) used for +## initialization, i.e., a sequence of length t+n is generated and +## x(n+1:t+n) is returned. If n is omitted, n=100 is used. + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Simulate an ARMA process + +function x = arma_rnd (a, b, v, t, n) + + unwind_protect + orig_listelemok = empty_list_elements_ok; + empty_list_elements_ok = "true"; + + if (nargin == 4) + n = 100; + elseif (nargin == 5) + if (!is_scalar (t)) + error ("arma_rnd: n must be a scalar"); + endif + else + usage ("arma_rnd (a, b, v, t [, n])"); + endif + + if ( (min (size (a)) > 1) || (min (size (b)) > 1) ) + error ("arma_rnd: a and b must not be matrices"); + endif + + if (!is_scalar (t)) + error ("arma_rnd: t must be a scalar"); + endif + + ar = length (a); + br = length (b); + + a = reshape (a, ar, 1); + b = reshape (b, br, 1); + + a = [1; -a]; # apply our notational convention + b = [1; b]; + + n = min (n, ar + br); + + e = sqrt (v) * randn (t + n, 1); + + x = filter (b, a, e); + x = x(n + 1 : t + n); + + unwind_protect_cleanup + + empty_list_elements_ok = orig_listelemok; + + end_unwind_protect + +endfunction