Mercurial > hg > octave-lyh
annotate scripts/signal/arch_rnd.m @ 17295:011e5d67baf1
interpft.m: Widen tolerance on tests to allow for FFT error
* scripts/general/interpft.m: Widen tolerance on tests to allow for
error in FFT routines, for example when not linking with FFTW.
| author | Mike Miller <mtmiller@ieee.org> |
|---|---|
| date | Tue, 20 Aug 2013 09:46:17 -0400 |
| parents | 5d3a684236b0 |
| children | 1c89599167a6 |
| rev | line source |
|---|---|
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
1 ## Copyright (C) 1995-2012 Kurt Hornik |
| 3426 | 2 ## |
| 3922 | 3 ## This file is part of Octave. |
| 4 ## | |
| 5 ## Octave is free software; you can redistribute it and/or modify it | |
| 6 ## under the terms of the GNU General Public License as published by | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 3426 | 9 ## |
| 3922 | 10 ## Octave is distributed in the hope that it will be useful, but |
| 3191 | 11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
| 12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 3426 | 13 ## General Public License for more details. |
| 14 ## | |
| 3191 | 15 ## You should have received a copy of the GNU General Public License |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 3191 | 18 |
| 3449 | 19 ## -*- texinfo -*- |
| 3500 | 20 ## @deftypefn {Function File} {} arch_rnd (@var{a}, @var{b}, @var{t}) |
| 21 ## Simulate an ARCH sequence of length @var{t} with AR | |
| 3449 | 22 ## coefficients @var{b} and CH coefficients @var{a}. I.e., the result |
| 3500 | 23 ## @math{y(t)} follows the model |
|
10846
a4f482e66b65
Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents:
9245
diff
changeset
|
24 ## @c Set example in small font to prevent overfull line |
| 3191 | 25 ## |
| 7031 | 26 ## @smallexample |
| 3499 | 27 ## y(t) = b(1) + b(2) * y(t-1) + @dots{} + b(lb) * y(t-lb+1) + e(t), |
| 7031 | 28 ## @end smallexample |
| 3449 | 29 ## |
| 30 ## @noindent | |
| 3499 | 31 ## where @math{e(t)}, given @var{y} up to time @math{t-1}, is |
| 32 ## @math{N(0, h(t))}, with | |
|
10846
a4f482e66b65
Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents:
9245
diff
changeset
|
33 ## @c Set example in small font to prevent overfull line |
| 3449 | 34 ## |
| 7031 | 35 ## @smallexample |
| 3499 | 36 ## h(t) = a(1) + a(2) * e(t-1)^2 + @dots{} + a(la) * e(t-la+1)^2 |
| 7031 | 37 ## @end smallexample |
| 3449 | 38 ## @end deftypefn |
| 3191 | 39 |
| 5428 | 40 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> |
| 3457 | 41 ## Description: Simulate an ARCH process |
| 3191 | 42 |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
43 function y = arch_rnd (a, b, t) |
| 3426 | 44 |
| 3191 | 45 if (nargin != 3) |
| 6046 | 46 print_usage (); |
| 3191 | 47 endif |
| 3426 | 48 |
| 3457 | 49 if (! ((min (size (a)) == 1) && (min (size (b)) == 1))) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11469
diff
changeset
|
50 error ("arch_rnd: A and B must both be scalars or vectors"); |
| 3191 | 51 endif |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
52 if (! (isscalar (t) && (t > 0) && (rem (t, 1) == 0))) |
| 3457 | 53 error ("arch_rnd: T must be a positive integer"); |
| 3191 | 54 endif |
| 3426 | 55 |
| 3457 | 56 if (! (a(1) > 0)) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11469
diff
changeset
|
57 error ("arch_rnd: A(1) must be positive"); |
| 3191 | 58 endif |
| 59 ## perhaps add a test for the roots of a(z) here ... | |
| 3426 | 60 |
| 3191 | 61 la = length (a); |
| 62 a = reshape (a, 1, la); | |
| 63 if (la == 1) | |
| 64 a = [a, 0]; | |
| 65 la = la + 1; | |
| 66 endif | |
| 7436 | 67 |
| 3191 | 68 lb = length (b); |
| 69 b = reshape (b, 1, lb); | |
| 70 if (lb == 1) | |
| 71 b = [b, 0]; | |
| 72 lb = lb + 1; | |
| 73 endif | |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
74 m = max ([la, lb]); |
| 3426 | 75 |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
76 e = zeros (t, 1); |
|
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
77 h = zeros (t, 1); |
|
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
78 y = zeros (t, 1); |
| 3426 | 79 |
| 3191 | 80 h(1) = a(1); |
| 81 e(1) = sqrt (h(1)) * randn; | |
| 82 y(1) = b(1) + e(1); | |
| 3426 | 83 |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
84 for t = 2:m |
| 3238 | 85 ta = min ([t, la]); |
| 7436 | 86 h(t) = a(1) + a(2:ta) * e(t-ta+1:t-1).^2; |
| 3191 | 87 e(t) = sqrt (h(t)) * randn; |
| 3238 | 88 tb = min ([t, lb]); |
| 7436 | 89 y(t) = b(1) + b(2:tb) * y(t-tb+1:t-1) + e(t); |
| 3191 | 90 endfor |
| 7436 | 91 |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
92 if (t > m) |
|
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
93 for t = m+1:t |
| 7436 | 94 h(t) = a(1) + a(2:la) * e(t-la+1:t-1).^2; |
| 3191 | 95 e(t) = sqrt (h(t)) * randn; |
| 7436 | 96 y(t) = b(1) + b(2:lb) * y(t-tb+1:t-1) + e(t); |
| 3191 | 97 endfor |
| 98 endif | |
| 3426 | 99 |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
100 y = y(1:t); |
| 3426 | 101 |
| 3191 | 102 endfunction |