usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/000755 000765 000120 00000000000 12402246361 023717 5ustar00jburnettadmin000000 000000 usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/arch_fit.m000644 000765 000120 00000006712 12402246361 025662 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Kurt Hornik ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{a}, @var{b}] =} arch_fit (@var{y}, @var{x}, @var{p}, @var{iter}, @var{gamma}, @var{a0}, @var{b0}) ## Fit an ARCH regression model to the time series @var{y} using the ## scoring algorithm in Engle's original ARCH paper. The model is ## ## @example ## @group ## y(t) = b(1) * x(t,1) + @dots{} + b(k) * x(t,k) + e(t), ## h(t) = a(1) + a(2) * e(t-1)^2 + @dots{} + a(p+1) * e(t-p)^2 ## @end group ## @end example ## ## @noindent ## in which @math{e(t)} is @math{N(0, h(t))}, given a time-series vector ## @var{y} up to time @math{t-1} and a matrix of (ordinary) regressors ## @var{x} up to @math{t}. The order of the regression of the residual ## variance is specified by @var{p}. ## ## If invoked as @code{arch_fit (@var{y}, @var{k}, @var{p})} with a ## positive integer @var{k}, fit an ARCH(@var{k}, @var{p}) process, ## i.e., do the above with the @math{t}-th row of @var{x} given by ## ## @example ## [1, y(t-1), @dots{}, y(t-k)] ## @end example ## ## Optionally, one can specify the number of iterations @var{iter}, the ## updating factor @var{gamma}, and initial values @math{a0} and ## @math{b0} for the scoring algorithm. ## @end deftypefn ## Author: KH ## Description: Fit an ARCH regression model function [a, b] = arch_fit (y, x, p, iter, gamma, a0, b0) if ((nargin < 3) || (nargin == 6) || (nargin > 7)) print_usage (); endif if (! (isvector (y))) error ("arch_fit: Y must be a vector"); endif T = length (y); y = reshape (y, T, 1); [rx, cx] = size (x); if ((rx == 1) && (cx == 1)) x = autoreg_matrix (y, x); elseif (! (rx == T)) error ("arch_fit: either rows (X) == length (Y), or X is a scalar"); endif [T, k] = size (x); if (nargin == 7) a = a0; b = b0; e = y - x * b; else [b, v_b, e] = ols (y, x); a = [v_b, (zeros (1, p))]'; if (nargin < 5) gamma = 0.1; if (nargin < 4) iter = 50; endif endif endif esq = e.^2; Z = autoreg_matrix (esq, p); for i = 1 : iter; h = Z * a; tmp = esq ./ h.^2 - 1 ./ h; s = 1 ./ h(1:T-p); for j = 1 : p; s = s - a(j+1) * tmp(j+1:T-p+j); endfor r = 1 ./ h(1:T-p); for j = 1:p; r = r + 2 * h(j+1:T-p+j).^2 .* esq(1:T-p); endfor r = sqrt (r); X_tilde = x(1:T-p, :) .* (r * ones (1,k)); e_tilde = e(1:T-p) .*s ./ r; delta_b = inv (X_tilde' * X_tilde) * X_tilde' * e_tilde; b = b + gamma * delta_b; e = y - x * b; esq = e .^ 2; Z = autoreg_matrix (esq, p); h = Z * a; f = esq ./ h - ones (T,1); Z_tilde = Z ./ (h * ones (1, p+1)); delta_a = inv (Z_tilde' * Z_tilde) * Z_tilde' * f; a = a + gamma * delta_a; endfor endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/arch_rnd.m000644 000765 000120 00000005307 12402246361 025662 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Kurt Hornik ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} arch_rnd (@var{a}, @var{b}, @var{t}) ## Simulate an ARCH sequence of length @var{t} with AR ## coefficients @var{b} and CH coefficients @var{a}. I.e., the result ## @math{y(t)} follows the model ## @c Set example in small font to prevent overfull line ## ## @smallexample ## y(t) = b(1) + b(2) * y(t-1) + @dots{} + b(lb) * y(t-lb+1) + e(t), ## @end smallexample ## ## @noindent ## where @math{e(t)}, given @var{y} up to time @math{t-1}, is ## @math{N(0, h(t))}, with ## @c Set example in small font to prevent overfull line ## ## @smallexample ## h(t) = a(1) + a(2) * e(t-1)^2 + @dots{} + a(la) * e(t-la+1)^2 ## @end smallexample ## @end deftypefn ## Author: KH ## Description: Simulate an ARCH process function y = arch_rnd (a, b, t) if (nargin != 3) print_usage (); endif if (! ((min (size (a)) == 1) && (min (size (b)) == 1))) error ("arch_rnd: A and B must both be scalars or vectors"); endif if (! (isscalar (t) && (t > 0) && (rem (t, 1) == 0))) error ("arch_rnd: T must be a positive integer"); endif if (! (a(1) > 0)) error ("arch_rnd: A(1) must be positive"); endif ## perhaps add a test for the roots of a(z) here ... la = length (a); a = reshape (a, 1, la); if (la == 1) a = [a, 0]; la = la + 1; endif lb = length (b); b = reshape (b, 1, lb); if (lb == 1) b = [b, 0]; lb = lb + 1; endif m = max ([la, lb]); e = zeros (t, 1); h = zeros (t, 1); y = zeros (t, 1); h(1) = a(1); e(1) = sqrt (h(1)) * randn; y(1) = b(1) + e(1); for t = 2:m ta = min ([t, la]); h(t) = a(1) + a(2:ta) * e(t-ta+1:t-1).^2; e(t) = sqrt (h(t)) * randn; tb = min ([t, lb]); y(t) = b(1) + b(2:tb) * y(t-tb+1:t-1) + e(t); endfor if (t > m) for t = m+1:t h(t) = a(1) + a(2:la) * e(t-la+1:t-1).^2; e(t) = sqrt (h(t)) * randn; y(t) = b(1) + b(2:lb) * y(t-tb+1:t-1) + e(t); endfor endif y = y(1:t); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/arch_test.m000644 000765 000120 00000005521 12402246361 026054 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Kurt Hornik ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{pval}, @var{lm}] =} arch_test (@var{y}, @var{x}, @var{p}) ## For a linear regression model ## ## @example ## y = x * b + e ## @end example ## ## @noindent ## perform a Lagrange Multiplier (LM) test of the null hypothesis of no ## conditional heteroscedascity against the alternative of CH(@var{p}). ## ## I.e., the model is ## ## @example ## y(t) = b(1) * x(t,1) + @dots{} + b(k) * x(t,k) + e(t), ## @end example ## ## @noindent ## given @var{y} up to @math{t-1} and @var{x} up to @math{t}, ## @math{e}(t) is @math{N(0, h(t))} with ## ## @example ## h(t) = v + a(1) * e(t-1)^2 + @dots{} + a(p) * e(t-p)^2, ## @end example ## ## @noindent ## and the null is @math{a(1)} == @dots{} == @math{a(p)} == 0. ## ## If the second argument is a scalar integer, @math{k}, perform the same ## test in a linear autoregression model of order @math{k}, i.e., with ## ## @example ## [1, y(t-1), @dots{}, y(t-@var{k})] ## @end example ## ## @noindent ## as the @math{t}-th row of @var{x}. ## ## Under the null, LM approximately has a chisquare distribution with ## @var{p} degrees of freedom and @var{pval} is the @math{p}-value (1 ## minus the CDF of this distribution at LM) of the test. ## ## If no output argument is given, the @math{p}-value is displayed. ## @end deftypefn ## Author: KH ## Description: Test for conditional heteroscedascity function [pval, lm] = arch_test (y, x, p) if (nargin != 3) error ("arch_test: 3 input arguments required"); endif if (! (isvector (y))) error ("arch_test: Y must be a vector"); endif T = length (y); y = reshape (y, T, 1); [rx, cx] = size (x); if ((rx == 1) && (cx == 1)) x = autoreg_matrix (y, x); elseif (! (rx == T)) error ("arch_test: either rows (X) == length (Y), or X is a scalar"); endif if (! (isscalar (p) && (rem (p, 1) == 0) && (p > 0))) error ("arch_test: P must be a positive integer"); endif [b, v_b, e] = ols (y, x); Z = autoreg_matrix (e.^2, p); f = e.^2 / v_b - ones (T, 1); f = Z' * f; lm = f' * inv (Z'*Z) * f / 2; pval = 1 - chi2cdf (lm, p); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/arma_rnd.m000644 000765 000120 00000004432 12402246361 025663 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} arma_rnd (@var{a}, @var{b}, @var{v}, @var{t}, @var{n}) ## Return a simulation of the ARMA model ## ## @example ## @group ## x(n) = a(1) * x(n-1) + @dots{} + a(k) * x(n-k) ## + e(n) + b(1) * e(n-1) + @dots{} + b(l) * e(n-l) ## @end group ## @end example ## ## @noindent ## in which @var{k} is the length of vector @var{a}, @var{l} is the ## length of vector @var{b} and @var{e} is Gaussian white noise with ## variance @var{v}. The function returns a vector of length @var{t}. ## ## The optional parameter @var{n} gives the number of dummy ## @var{x}(@var{i}) used for initialization, i.e., a sequence of length ## @var{t}+@var{n} is generated and @var{x}(@var{n}+1:@var{t}+@var{n}) ## is returned. If @var{n} is omitted, @var{n} = 100 is used. ## @end deftypefn ## Author: FL ## Description: Simulate an ARMA process function x = arma_rnd (a, b, v, t, n) if (nargin == 4) n = 100; elseif (nargin == 5) if (!isscalar (n)) error ("arma_rnd: N must be a scalar"); endif else print_usage (); endif if ((min (size (a)) > 1) || (min (size (b)) > 1)) error ("arma_rnd: A and B must not be matrices"); endif if (!isscalar (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); ## Apply our notational convention. a = [1; -a]; 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); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/autoreg_matrix.m000644 000765 000120 00000003453 12402246361 027134 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Kurt Hornik ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} autoreg_matrix (@var{y}, @var{k}) ## Given a time series (vector) @var{y}, return a matrix with ones in the ## first column and the first @var{k} lagged values of @var{y} in the ## other columns. I.e., for @var{t} > @var{k}, @code{[1, ## @var{y}(@var{t}-1), @dots{}, @var{y}(@var{t}-@var{k})]} is the t-th row ## of the result. The resulting matrix may be used as a regressor matrix ## in autoregressions. ## @end deftypefn ## Author: KH ## Description: Design matrix for autoregressions function X = autoreg_matrix (y, k) if (nargin != 2) print_usage (); endif if (! (isvector (y))) error ("autoreg_matrix: Y must be a vector"); endif T = length (y); y = reshape (y, T, 1); X = ones (T, k+1); for j = 1 : k; X(:, j+1) = [(zeros (j, 1)); y(1:T-j)]; endfor endfunction %!test %! K = 4; %! A = zeros (1,K+1); %! A(1) = 1; %! B = eye (K+1); %! B(:,1) = 1; %! assert (autoreg_matrix (A,K), B); %!error autoreg_matrix () %!error autoreg_matrix (1) %!error autoreg_matrix (ones (4,1), 5) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/bartlett.m000644 000765 000120 00000003424 12402246361 025721 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Andreas Weingessel ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} bartlett (@var{m}) ## Return the filter coefficients of a Bartlett (triangular) window of ## length @var{m}. ## ## For a definition of the Bartlett window, see e.g., A. V. Oppenheim & ## R. W. Schafer, @cite{Discrete-Time Signal Processing}. ## @end deftypefn ## Author: AW ## Description: Coefficients of the Bartlett (triangular) window function c = bartlett (m) if (nargin != 1) print_usage (); endif if (! (isscalar (m) && (m == fix (m)) && (m > 0))) error ("bartlett: M has to be an integer > 0"); endif if (m == 1) c = 1; else m = m - 1; n = fix (m / 2); c = [2*(0:n)/m, 2-2*(n+1:m)/m]'; endif endfunction %!assert (bartlett (1), 1) %!assert (bartlett (2), zeros (2,1)) %!assert (bartlett (16), fliplr (bartlett (16))) %!assert (bartlett (15), fliplr (bartlett (15))) %!test %! N = 9; %! A = bartlett (N); %! assert (A(ceil (N/2)), 1); %!error bartlett () %!error bartlett (0.5) %!error bartlett (-1) %!error bartlett (ones (1,4)) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/blackman.m000644 000765 000120 00000003530 12402246361 025646 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Andreas Weingessel ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} blackman (@var{m}) ## Return the filter coefficients of a Blackman window of length @var{m}. ## ## For a definition of the Blackman window, see e.g., A. V. Oppenheim & ## R. W. Schafer, @cite{Discrete-Time Signal Processing}. ## @end deftypefn ## Author: AW ## Description: Coefficients of the Blackman window function c = blackman (m) if (nargin != 1) print_usage (); endif if (! (isscalar (m) && (m == fix (m)) && (m > 0))) error ("blackman: M has to be an integer > 0"); endif if (m == 1) c = 1; else m = m - 1; k = (0 : m)' / m; c = 0.42 - 0.5 * cos (2 * pi * k) + 0.08 * cos (4 * pi * k); endif endfunction %!assert (blackman (1), 1) %!assert (blackman (2), zeros (2,1), 1e-6) %!assert (blackman (16), fliplr (blackman (16))) %!assert (blackman (15), fliplr (blackman (15))) %!test %! N = 9; %! A = blackman (N); %! assert (A (ceil (N/2)), 1, 1e-6); %! assert ([A(1), A(length (A))], zeros (1, 2), 1e-6); %!error blackman () %!error blackman (0.5) %!error blackman (-1) %!error blackman (ones (1,4)) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/detrend.m000644 000765 000120 00000005005 12402246361 025522 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Kurt Hornik ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} detrend (@var{x}, @var{p}) ## If @var{x} is a vector, @code{detrend (@var{x}, @var{p})} removes the ## best fit of a polynomial of order @var{p} from the data @var{x}. ## ## If @var{x} is a matrix, @code{detrend (@var{x}, @var{p})} does the same ## for each column in @var{x}. ## ## The second argument is optional. If it is not specified, a value of 1 ## is assumed. This corresponds to removing a linear trend. ## ## The order of the polynomial can also be given as a string, in which case ## @var{p} must be either @qcode{"constant"} (corresponds to ## @code{@var{p}=0}) or ## @qcode{"linear"} (corresponds to @code{@var{p}=1}). ## @seealso{polyfit} ## @end deftypefn ## Author: KH ## Created: 11 October 1994 ## Adapted-By: jwe function y = detrend (x, p = 1) ## Check input if (nargin > 0 && isreal (x) && ndims (x) <= 2) ## Check p if (ischar (p) && strcmpi (p, "constant")) p = 0; elseif (ischar (p) && strcmpi (p, "linear")) p = 1; elseif (!isscalar (p) || p < 0 || p != fix (p)) error ("detrend: second input argument must be 'constant', 'linear' or a positive integer"); endif else error ("detrend: first input argument must be a real vector or matrix"); endif [m, n] = size (x); if (m == 1) x = x'; endif r = rows (x); b = ((1 : r)' * ones (1, p + 1)) .^ (ones (r, 1) * (0 : p)); y = x - b * (b \ x); if (m == 1) y = y'; endif endfunction %!test %! N = 32; %! x = (0:1:N-1)/N + 2; %! y = detrend (x); %! assert (abs (y(:)) < 20*eps); %!test %! N = 32; %! t = (0:1:N-1)/N; %! x = t .* t + 2; %! y = detrend (x,2); %! assert (abs (y(:)) < 30*eps); %!test %! N = 32; %! t = (0:1:N-1)/N; %! x = [t;4*t-3]'; %! y = detrend (x); %! assert (abs (y(:)) < 20*eps); usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/diffpara.m000644 000765 000120 00000004741 12402246361 025657 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{d}, @var{dd}] =} diffpara (@var{x}, @var{a}, @var{b}) ## Return the estimator @var{d} for the differencing parameter of an ## integrated time series. ## ## The frequencies from @math{[2*pi*a/t, 2*pi*b/T]} are used for the ## estimation. If @var{b} is omitted, the interval ## @math{[2*pi/T, 2*pi*a/T]} is used. If both @var{b} and @var{a} are ## omitted then @math{a = 0.5 * sqrt (T)} and @math{b = 1.5 * sqrt (T)} ## is used, where @math{T} is the sample size. If @var{x} is a matrix, ## the differencing parameter of each column is estimated. ## ## The estimators for all frequencies in the intervals ## described above is returned in @var{dd}. The value of @var{d} is ## simply the mean of @var{dd}. ## ## Reference: P.J. Brockwell & R.A. Davis. @cite{Time Series: ## Theory and Methods}. Springer 1987. ## @end deftypefn ## Author: FL ## Description: Estimate the fractional differencing parameter function [d, dd] = diffpara (x, a, b) if ((nargin < 1) || (nargin > 3)) print_usage (); else if (isvector (x)) n = length (x); k = 1; x = reshape (x, n, 1); else [n, k] = size (x); endif if (nargin == 1) a = 0.5 * sqrt (n); b = 1.5 * sqrt (n); elseif (nargin == 2) b = a; a = 1; endif endif if (! (isscalar (a) && isscalar (b))) error ("diffpara: A and B must be scalars"); endif dd = zeros (b - a + 1, k); for l = 1:k w = 2 * pi * (1 : n-1) / n; x = 2 * log (abs (1 - exp (-i*w))); y = log (periodogram (x(2:n,l))); x = center (x); y = center (y); for m = a:b dd(m-a+1) = - x(1:m) * y(1:m) / sumsq (x(1:m)); endfor endfor d = mean (dd); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/durbinlevinson.m000644 000765 000120 00000004651 12402246361 027144 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} durbinlevinson (@var{c}, @var{oldphi}, @var{oldv}) ## Perform one step of the Durbin-Levinson algorithm. ## ## The vector @var{c} specifies the autocovariances @code{[gamma_0, @dots{}, ## gamma_t]} from lag 0 to @var{t}, @var{oldphi} specifies the ## coefficients based on @var{c}(@var{t}-1) and @var{oldv} specifies the ## corresponding error. ## ## If @var{oldphi} and @var{oldv} are omitted, all steps from 1 to ## @var{t} of the algorithm are performed. ## @end deftypefn ## Author: FL ## Description: Perform one step of the Durbin-Levinson algorithm function [newphi, newv] = durbinlevinson (c, oldphi, oldv) if (! ((nargin == 1) || (nargin == 3))) print_usage (); endif if (columns (c) > 1) c = c'; endif newphi = 0; newv = 0; if (nargin == 3) t = length (oldphi) + 1; if (length (c) < t+1) error ("durbinlevinson: C too small"); endif if (oldv == 0) error ("durbinlevinson: OLDV = 0"); endif if (rows (oldphi) > 1) oldphi = oldphi'; endif newphi = zeros (1, t); newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv; for i = 2 : t newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1); endfor newv = (1 - newphi(1)^2) * oldv; elseif(nargin == 1) tt = length (c)-1; oldphi = c(2) / c(1); oldv = (1 - oldphi^2) * c(1); for t = 2 : tt newphi = zeros (1, t); newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv; for i = 2 : t newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1); endfor newv = (1 - newphi(1)^2) * oldv; oldv = newv; oldphi = newphi; endfor endif endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/fftconv.m000644 000765 000120 00000006251 12402246361 025546 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1994-2013 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} fftconv (@var{x}, @var{y}) ## @deftypefnx {Function File} {} fftconv (@var{x}, @var{y}, @var{n}) ## Convolve two vectors using the FFT for computation. ## ## @code{c = fftconv (@var{x}, @var{y})} returns a vector of length equal to ## @code{length (@var{x}) + length (@var{y}) - 1}. ## If @var{x} and @var{y} are the coefficient vectors of two polynomials, the ## returned value is the coefficient vector of the product polynomial. ## ## The computation uses the FFT by calling the function @code{fftfilt}. If ## the optional argument @var{n} is specified, an N-point FFT is used. ## @seealso{deconv, conv, conv2} ## @end deftypefn ## Author: KH ## Created: 3 September 1994 ## Adapted-By: jwe function c = fftconv (x, y, n) if (nargin < 2 || nargin > 3) print_usage (); endif if (! (isvector (x) && isvector (y))) error ("fftconv: both A and B must be vectors"); endif la = length (x); lb = length (y); if ((la == 1) || (lb == 1)) c = x * y; else lc = la + lb - 1; x(lc) = 0; y(lc) = 0; if (nargin == 2) c = fftfilt (x, y); else if (! isscalar (n)) error ("fftconv: N must be a scalar"); endif c = fftfilt (x, y, n); endif endif endfunction %% FIXME: Borrow tests from conv.m. May need a tolerance on the assert comparison %!test %! x = ones (3,1); %! y = ones (1,3); %! b = 2; %! c = 3; %! assert (fftconv (x, x), [1; 2; 3; 2; 1], 5*eps); %! assert (fftconv (y, y), [1, 2, 3, 2, 1], 5*eps); %! assert (fftconv (x, y), [1, 2, 3, 2, 1], 5*eps); %! assert (fftconv (y, x), [1; 2; 3; 2; 1], 5*eps); %! assert (fftconv (c, x), [3; 3; 3], 5*eps); %! assert (fftconv (c, y), [3, 3, 3], 5*eps); %! assert (fftconv (x, c), [3; 3; 3], 5*eps); %! assert (fftconv (y, c), [3, 3, 3], 5*eps); %! assert (fftconv (b, c), 6, 5*eps); %!test %! a = 1:10; %! b = 1:3; %! assert (size (conv (a,b)), [1, numel(a)+numel(b)-1]); %! assert (size (conv (b,a)), [1, numel(a)+numel(b)-1]); %! a = (1:10).'; %! b = 1:3; %! assert (size (conv (a,b)), [numel(a)+numel(b)-1, 1]); %! assert (size (conv (b,a)), [numel(a)+numel(b)-1, 1]); %!test %! a = 1:10; %! b = (1:3).'; %! assert (size (conv (a,b)), [1, numel(a)+numel(b)-1]); %! assert (size (conv (b,a)), [1, numel(a)+numel(b)-1]); %% Test input validation %!error fftconv (1) %!error fftconv (1,2,3,4) %!error fftconv ([1, 2; 3, 4], 3) %!error fftconv (2, []) %!error fftconv ([1,1], [2,2] , [3, 4]) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/fftfilt.m000644 000765 000120 00000013222 12402246361 025533 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1994-2013 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} fftfilt (@var{b}, @var{x}) ## @deftypefnx {Function File} {} fftfilt (@var{b}, @var{x}, @var{n}) ## ## With two arguments, @code{fftfilt} filters @var{x} with the FIR filter ## @var{b} using the FFT. ## ## Given the optional third argument, @var{n}, @code{fftfilt} uses the ## overlap-add method to filter @var{x} with @var{b} using an N-point ## FFT@. The FFT size must be an even power of 2 and must be greater than ## or equal to the length of @var{b}. If the specified @var{n} does not ## meet these criteria, it is automatically adjusted to the nearest value ## that does. ## ## If @var{x} is a matrix, filter each column of the matrix. ## @seealso{filter, filter2} ## @end deftypefn ## Author: Kurt Hornik ## Created: 3 September 1994 ## Adapted-By: jwe function y = fftfilt (b, x, n) ## If N is not specified explicitly, we do not use the overlap-add ## method at all because loops are really slow. Otherwise, we only ## ensure that the number of points in the FFT is the smallest power ## of two larger than N and length(b). This could result in length ## one blocks, but if the user knows better ... if (nargin < 2 || nargin > 3) print_usage (); endif transpose = (rows (x) == 1); if (transpose) x = x.'; endif [r_x, c_x] = size (x); [r_b, c_b] = size (b); if (! isvector (b)) error ("fftfilt: B must be a vector"); endif if (ndims (x) != 2) error ("fftfilt: X must be a 1-D or 2-D array"); endif l_b = r_b * c_b; b = reshape (b, l_b, 1); if (nargin == 2) ## Use FFT with the smallest power of 2 which is >= length (x) + ## length (b) - 1 as number of points ... n = 2 ^ nextpow2 (r_x + l_b - 1); B = fft (b, n); y = ifft (fft (x, n) .* B(:, ones (1, c_x))); else ## Use overlap-add method ... if (! (isscalar (n))) error ("fftfilt: N has to be a scalar"); endif n = 2 ^ nextpow2 (max ([n, l_b])); L = n - l_b + 1; B = fft (b, n); B = B(:, ones (c_x,1)); R = ceil (r_x / L); y = zeros (r_x, c_x); for r = 1:R; lo = (r - 1) * L + 1; hi = min (r * L, r_x); tmp = zeros (n, c_x); tmp(1:(hi-lo+1),:) = x(lo:hi,:); tmp = ifft (fft (tmp) .* B); hi = min (lo+n-1, r_x); y(lo:hi,:) = y(lo:hi,:) + tmp(1:(hi-lo+1),:); endfor endif y = y(1:r_x, :); ## Final cleanups: ## - If both b and x are real, y should be real. ## - If b is real and x is imaginary, y should be imaginary. ## - If b is imaginary and x is real, y should be imaginary. ## - If both b and x are imaginary, y should be real. xisreal = all (imag (x) == 0); xisimag = all (real (x) == 0); if (all (imag (b) == 0)) y (:,xisreal) = real (y (:,xisreal)); y (:,xisimag) = complex (real (y (:,xisimag)) * 0, imag (y (:,xisimag))); elseif (all (real (b) == 0)) y (:,xisreal) = complex (real (y (:,xisreal)) * 0, imag (y (:,xisreal))); y (:,xisimag) = real (y (:,xisimag)); endif ## - If both x and b are integer in both real and imaginary ## components, y should be integer. if (! any (b - fix (b))) idx = find (! any (x - fix (x))); y (:, idx) = round (y (:, idx)); endif ## Transpose after cleanup, otherwise rounding fails. if (transpose) y = y.'; endif endfunction %!shared b, x, r %!test %! b = [1 1]; %! x = [1, zeros(1,9)]; %! assert (fftfilt (b, x ), [1 1 0 0 0 0 0 0 0 0] ); %! assert (fftfilt (b, x.'), [1 1 0 0 0 0 0 0 0 0].'); %! assert (fftfilt (b.',x ), [1 1 0 0 0 0 0 0 0 0] ); %! assert (fftfilt (b.',x.'), [1 1 0 0 0 0 0 0 0 0].'); %! assert (fftfilt (b, [x.' x.']), [1 1 0 0 0 0 0 0 0 0].'*[1 1]); %! assert (fftfilt (b, [x.'+2*eps x.']) == [1 1 0 0 0 0 0 0 0 0].'*[1 1], [false(10, 1) true(10, 1)]); %!test %! r = sqrt (1/2) * (1+i); %! b = b*r; %! assert (fftfilt (b, x ), r*[1 1 0 0 0 0 0 0 0 0] , eps ); %! assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], 2*eps); %! assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps ); %!test %! b = [1 1]; %! x = zeros (10,3); x(1,1)=-1; x(1,2)=1; %! y0 = zeros (10,3); y0(1:2,1)=-1; y0(1:2,2)=1; %! y = fftfilt (b, x); %! assert (y0, y); %! y = fftfilt (b*i, x); %! assert (y0*i, y); %! y = fftfilt (b, x*i); %! assert (y0*i, y); %! y = fftfilt (b*i, x*i); %! assert (-y0, y); %! x = rand (10, 1); %! y = fftfilt (b, [x x*i]); %! assert (true, isreal (y(:,1))); %! assert (false, any (real (y(:,2)))); %!test %! b = rand (10, 1); %! x = rand (10, 1); %! y0 = filter (b, 1, x); %! y = fftfilt (b, x); %! assert (y0, y, 16*eps); %! y0 = filter (b*i, 1, x*i); %! y = fftfilt (b*i, x*i); %! assert (y0, y, 16*eps); %!test %! b = rand (10, 1) + i*rand (10, 1); %! x = rand (10, 1) + i*rand (10, 1); %! y0 = filter (b, 1, x); %! y = fftfilt (b, x); %! assert (y0, y, 55*eps); %% Test input validation %!error fftfilt (1) %!error fftfilt (1, 2, 3, 4) %!error fftfilt (ones (2), 1) %!error fftfilt (2, ones (3,3,3)) %!error fftfilt (2, 1, ones (2)) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/fftshift.m000644 000765 000120 00000007001 12402246361 025710 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1997-2013 Vincent Cautaerts ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} fftshift (@var{x}) ## @deftypefnx {Function File} {} fftshift (@var{x}, @var{dim}) ## Perform a shift of the vector @var{x}, for use with the @code{fft} ## and @code{ifft} functions, in order the move the frequency 0 to the ## center of the vector or matrix. ## ## If @var{x} is a vector of @math{N} elements corresponding to @math{N} ## time samples spaced by @nospell{@math{dt}}, then ## @code{fftshift (fft (@var{x}))} corresponds to frequencies ## ## @example ## f = [ -(ceil((N-1)/2):-1:1)*df 0 (1:floor((N-1)/2))*df ] ## @end example ## ## @noindent ## where @nospell{@math{df} = 1 / @math{dt}}. ## ## If @var{x} is a matrix, the same holds for rows and columns. If ## @var{x} is an array, then the same holds along each dimension. ## ## The optional @var{dim} argument can be used to limit the dimension ## along which the permutation occurs. ## @end deftypefn ## Author: Vincent Cautaerts ## Created: July 1997 ## Adapted-By: jwe function retval = fftshift (x, dim) if (nargin != 1 && nargin != 2) print_usage (); endif if (nargin == 2) if (! (isscalar (dim) && dim > 0 && dim == fix (dim))) error ("fftshift: dimension DIM must be a positive integer"); endif nd = ndims (x); sz = size (x); sz2 = ceil (sz(dim) / 2); idx = cell (); idx = repmat ({':'}, nd, 1); idx{dim} = [sz2+1:sz(dim), 1:sz2]; retval = x(idx{:}); else if (isvector (x)) xl = length (x); xx = ceil (xl/2); retval = x([xx+1:xl, 1:xx]); elseif (ismatrix (x)) nd = ndims (x); sz = size (x); sz2 = ceil (sz ./ 2); idx = cell (); for i = 1:nd idx{i} = [sz2(i)+1:sz(i), 1:sz2(i)]; endfor retval = x(idx{:}); else error ("fftshift: expecting vector or matrix argument"); endif endif endfunction %!test %! x = [0:7]; %! y = fftshift (x); %! assert (y, [4 5 6 7 0 1 2 3]); %! assert (fftshift (y), x); %!test %! x = [0:6]; %! y = fftshift (x); %! assert (y, [4 5 6 0 1 2 3]); %! assert (fftshift (y), [1 2 3 4 5 6 0]); %!test %! x = [0:7]'; %! y = fftshift (x); %! assert (y, [4;5;6;7;0;1;2;3]); %! assert (fftshift (y), x); %!test %! x = [0:6]'; %! y = fftshift (x); %! assert (y, [4;5;6;0;1;2;3]); %! assert (fftshift (y), [1;2;3;4;5;6;0]); %!test %! x = [0:3]; %! x = [x;2*x;3*x+1;4*x+1]; %! y = fftshift (x); %! assert (y, [[7 10 1 4];[9 13 1 5];[2 3 0 1];[4 6 0 2]]); %! assert (fftshift (y), x); %!test %! x = [0:3]; %! x = [x;2*x;3*x+1;4*x+1]; %! y = fftshift (x,1); %! assert (y, [[1 4 7 10];[1 5 9 13];[0 1 2 3];[0 2 4 6]]); %! assert (fftshift (y,1), x); %!test %! x = [0:3]; %! x = [x;2*x;3*x+1;4*x+1]; %! y = fftshift (x,2); %! assert (y, [[2 3 0 1];[4 6 0 2];[7 10 1 4];[9 13 1 5]]); %! assert (fftshift (y,2), x); usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/filter2.m000644 000765 000120 00000003376 12402246361 025455 0ustar00jburnettadmin000000 000000 ## Copyright (C) 2001-2013 Paul Kienzle ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {@var{y} =} filter2 (@var{b}, @var{x}) ## @deftypefnx {Function File} {@var{y} =} filter2 (@var{b}, @var{x}, @var{shape}) ## Apply the 2-D FIR filter @var{b} to @var{x}. If the argument ## @var{shape} is specified, return an array of the desired shape. ## Possible values are: ## ## @table @asis ## @item @qcode{"full"} ## pad @var{x} with zeros on all sides before filtering. ## ## @item @qcode{"same"} ## unpadded @var{x} (default) ## ## @item @qcode{"valid"} ## trim @var{x} after filtering so edge effects are no included. ## @end table ## ## Note this is just a variation on convolution, with the parameters ## reversed and @var{b} rotated 180 degrees. ## @seealso{conv2} ## @end deftypefn ## Author: Paul Kienzle ## 2001-02-08 ## * initial release function y = filter2 (b, x, shape) if (nargin < 2 || nargin > 3) print_usage (); endif if (nargin < 3) shape = "same"; endif [nr, nc] = size (b); y = conv2 (x, b(nr:-1:1, nc:-1:1), shape); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/fractdiff.m000644 000765 000120 00000003301 12402246361 026022 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} fractdiff (@var{x}, @var{d}) ## Compute the fractional differences @math{(1-L)^d x} where @math{L} ## denotes the lag-operator and @math{d} is greater than -1. ## @end deftypefn ## Author: FL ## Description: Compute fractional differences function retval = fractdiff (x, d) if (nargin != 2) print_usage (); endif N = 100; if (! isvector (x)) error ("fractdiff: X must be a vector"); endif if (! isscalar (d)) error ("fractdiff: D must be a scalar"); endif if (d >= 1) for k = 1 : d x = x(2 : length (x)) - x(1 : length (x) - 1); endfor endif if (d > -1) d = rem (d, 1); if (d != 0) n = (0 : N)'; w = real (gamma (-d+n) ./ gamma (-d) ./ gamma (n+1)); retval = fftfilt (w, x); retval = retval(1 : length (x)); else retval = x; endif else error ("fractdiff: D must be > -1"); endif endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/freqz.m000644 000765 000120 00000013467 12402246361 025237 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1994-2013 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}, @var{n}, "whole") ## @deftypefnx {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}) ## @deftypefnx {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}) ## @deftypefnx {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}, @var{n}) ## @deftypefnx {Function File} {@var{h} =} freqz (@var{b}, @var{a}, @var{w}) ## @deftypefnx {Function File} {[@dots{}] =} freqz (@dots{}, @var{Fs}) ## @deftypefnx {Function File} {} freqz (@dots{}) ## ## Return the complex frequency response @var{h} of the rational IIR filter ## whose numerator and denominator coefficients are @var{b} and @var{a}, ## respectively. The response is evaluated at @var{n} angular frequencies ## between 0 and ## @ifnottex ## 2*pi. ## @end ifnottex ## @tex ## $2\pi$. ## @end tex ## ## @noindent ## The output value @var{w} is a vector of the frequencies. ## ## If @var{a} is omitted, the denominator is assumed to be 1 (this ## corresponds to a simple FIR filter). ## ## If @var{n} is omitted, a value of 512 is assumed. ## For fastest computation, @var{n} should factor into a small number of ## small primes. ## ## If the fourth argument, @qcode{"whole"}, is omitted the response is ## evaluated at frequencies between 0 and ## @ifnottex ## pi. ## @end ifnottex ## @tex ## $\pi$. ## @end tex ## ## @code{freqz (@var{b}, @var{a}, @var{w})} ## ## Evaluate the response at the specific frequencies in the vector @var{w}. ## The values for @var{w} are measured in radians. ## ## @code{[@dots{}] = freqz (@dots{}, @var{Fs})} ## ## Return frequencies in Hz instead of radians assuming a sampling rate ## @var{Fs}. If you are evaluating the response at specific frequencies ## @var{w}, those frequencies should be requested in Hz rather than radians. ## ## @code{freqz (@dots{})} ## ## Plot the magnitude and phase response of @var{h} rather than returning them. ## ## @seealso{freqz_plot} ## @end deftypefn ## Author: jwe ??? function [h_r, f_r] = freqz (b, a, n, region, Fs) if (nargin < 1 || nargin > 5) print_usage (); elseif (nargin == 1) ## Response of an FIR filter. a = n = region = Fs = []; elseif (nargin == 2) ## Response of an IIR filter n = region = Fs = []; elseif (nargin == 3) region = Fs = []; elseif (nargin == 4) Fs = []; if (! ischar (region) && ! isempty (region)) Fs = region; region = []; endif endif if (isempty (b)) b = 1; endif if (isempty (a)) a = 1; endif if (isempty (n)) n = 512; endif if (isempty (region)) if (isreal (b) && isreal (a)) region = "half"; else region = "whole"; endif endif if (isempty (Fs)) freq_norm = true; if (nargout == 0) Fs = 2; else Fs = 2*pi; endif else freq_norm = false; endif a = a(:); b = b(:); if (! isscalar (n)) ## Explicit frequency vector given w = f = n; if (nargin == 4) ## Sampling rate Fs was specified w = 2*pi*f/Fs; endif k = max (length (b), length (a)); hb = polyval (postpad (b, k), exp (j*w)); ha = polyval (postpad (a, k), exp (j*w)); else ## polyval(fliplr(P),exp(jw)) is O(p n) and fft(x) is O(n log(n)), ## where p is the order of the polynomial P. For small p it ## would be faster to use polyval but in practice the overhead for ## polyval is much higher and the little bit of time saved isn't ## worth the extra code. k = max (length (b), length (a)); if (k > n/2 && nargout == 0) ## Ensure a causal phase response. n = n * 2 .^ ceil (log2 (2*k/n)); endif if (strcmp (region, "whole")) N = n; else N = 2*n; endif f = Fs * (0:n-1).' / N; pad_sz = N*ceil (k/N); b = postpad (b, pad_sz); a = postpad (a, pad_sz); hb = zeros (n, 1); ha = zeros (n, 1); for i = 1:N:pad_sz hb = hb + fft (postpad (b(i:i+N-1), N))(1:n); ha = ha + fft (postpad (a(i:i+N-1), N))(1:n); endfor endif h = hb ./ ha; if (nargout != 0) ## Return values and don't plot. h_r = h; f_r = f; else ## Plot and don't return values. freqz_plot (f, h, freq_norm); endif endfunction %!test # correct values and fft-polyval consistency %! ## butterworth filter, order 2, cutoff pi/2 radians %! b = [0.292893218813452 0.585786437626905 0.292893218813452]; %! a = [1 0 0.171572875253810]; %! [h,w] = freqz (b,a,32); %! assert (h(1),1,10*eps); %! assert (abs (h(17)).^2,0.5,10*eps); %! assert (h,freqz (b,a,w),10*eps); # fft should be consistent with polyval %!test # whole-half consistency %! b = [1 1 1]/3; # 3-sample average %! [h,w] = freqz (b,1,32,"whole"); %! assert (h(2:16),conj (h(32:-1:18)),20*eps); %! [h2,w2] = freqz (b,1,16,"half"); %! assert (h(1:16),h2,20*eps); %! assert (w(1:16),w2,20*eps); %!test # Sampling frequency properly interpreted %! b = [1 1 1]/3; a = [1 0.2]; %! [h,f] = freqz (b,a,16,320); %! assert (f,[0:15]'*10,10*eps); %! [h2,f2] = freqz (b,a,[0:15]*10,320); %! assert (f2,[0:15]*10,10*eps); %! assert (h,h2.',20*eps); %! [h3,f3] = freqz (b,a,32,"whole",320); %! assert (f3,[0:31]'*10,10*eps); usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/freqz_plot.m000644 000765 000120 00000003643 12402246361 026270 0ustar00jburnettadmin000000 000000 ## Copyright (C) 2002-2013 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} freqz_plot (@var{w}, @var{h}) ## @deftypefnx {Function File} {} freqz_plot (@var{w}, @var{h}, @var{freq_norm}) ## Plot the magnitude and phase response of @var{h}. ## ## If the optional @var{freq_norm} argument is true, the frequency vector ## @var{w} is in units of normalized radians. If @var{freq_norm} is false, or ## not given, then @var{w} is measured in Hertz. ## @seealso{freqz} ## @end deftypefn ## Author: Paul Kienzle function freqz_plot (w, h, freq_norm = false) if (nargin < 2 || nargin > 3) print_usage (); endif n = length (w); ## ## exclude zero-frequency ## h = h (2 : length (h)); ## w = w (2 : length (w)); ## n = n-1; mag = 20 * log10 (abs (h)); phase = unwrap (arg (h)); if (freq_norm) x_label = 'Normalized Frequency (\times\pi rad/sample)'; else x_label = "Frequency (Hz)"; endif subplot (2, 1, 1); plot (w, mag); grid ("on"); axis ([w(1), w(n)], "autoy"); xlabel (x_label); ylabel ("Magnitude (dB)"); subplot (2, 1, 2); plot (w, phase*360/(2*pi)); grid ("on"); axis ([w(1), w(n)], "autoy"); xlabel (x_label); ylabel ("Phase (degrees)"); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/hamming.m000644 000765 000120 00000003354 12402246361 025522 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Andreas Weingessel ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} hamming (@var{m}) ## Return the filter coefficients of a Hamming window of length @var{m}. ## ## For a definition of the Hamming window, see e.g., A. V. Oppenheim & ## R. W. Schafer, @cite{Discrete-Time Signal Processing}. ## @end deftypefn ## Author: AW ## Description: Coefficients of the Hamming window function c = hamming (m) if (nargin != 1) print_usage (); endif if (! (isscalar (m) && (m == fix (m)) && (m > 0))) error ("hamming: M has to be an integer > 0"); endif if (m == 1) c = 1; else m = m - 1; c = 0.54 - 0.46 * cos (2 * pi * (0:m)' / m); endif endfunction %!assert (hamming (1), 1) %!assert (hamming (2), (0.54 - 0.46)*ones (2,1)) %!assert (hamming (16), fliplr (hamming (16))) %!assert (hamming (15), fliplr (hamming (15))) %!test %! N = 15; %! A = hamming (N); %! assert (A (ceil (N/2)), 1); %!error hamming () %!error hamming (0.5) %!error hamming (-1) %!error hamming (ones (1,4)) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/hanning.m000644 000765 000120 00000003340 12402246361 025517 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Andreas Weingessel ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} hanning (@var{m}) ## Return the filter coefficients of a Hanning window of length @var{m}. ## ## For a definition of this window type, see e.g., A. V. Oppenheim & ## R. W. Schafer, @cite{Discrete-Time Signal Processing}. ## @end deftypefn ## Author: AW ## Description: Coefficients of the Hanning window function c = hanning (m) if (nargin != 1) print_usage (); endif if (! (isscalar (m) && (m == fix (m)) && (m > 0))) error ("hanning: M has to be an integer > 0"); endif if (m == 1) c = 1; else m = m - 1; c = 0.5 - 0.5 * cos (2 * pi * (0 : m)' / m); endif endfunction %!assert (hanning (1), 1); %!assert (hanning (2), zeros (2,1)); %!assert (hanning (16), fliplr (hanning (16))); %!assert (hanning (15), fliplr (hanning (15))); %!test %! N = 15; %! A = hanning (N); %! assert (A(ceil (N/2)), 1); %!error hanning () %!error hanning (0.5) %!error hanning (-1) %!error hanning (ones (1,4)) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/hurst.m000644 000765 000120 00000002631 12402246361 025244 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} hurst (@var{x}) ## Estimate the Hurst parameter of sample @var{x} via the rescaled range ## statistic. If @var{x} is a matrix, the parameter is estimated for ## every single column. ## @end deftypefn ## Author: FL ## Description: Estimate the Hurst parameter function H = hurst (x) if (nargin != 1) print_usage (); endif if (isscalar (x)) error ("hurst: X must not be a scalar"); elseif (isvector (x)) x = reshape (x, length (x), 1); endif [xr, xc] = size (x); s = std (x); w = cumsum (x - mean (x)); RS = (max (w) - min (w)) ./ s; H = log (RS) / log (xr); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/ifftshift.m000644 000765 000120 00000006005 12402246361 026064 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1997-2013 Vincent Cautaerts ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} ifftshift (@var{x}) ## @deftypefnx {Function File} {} ifftshift (@var{x}, @var{dim}) ## Undo the action of the @code{fftshift} function. For even length ## @var{x}, @code{fftshift} is its own inverse, but odd lengths differ ## slightly. ## @end deftypefn ## Author: Vincent Cautaerts ## Created: July 1997 ## Adapted-By: jwe ## Modified-By: Paul Kienzle, converted from fftshift ## Modified-By: David Bateman, add NDArray capability and option dim arg function retval = ifftshift (x, dim) retval = 0; if (nargin != 1 && nargin != 2) print_usage (); endif if (nargin == 2) if (! isscalar (dim)) error ("ifftshift: dimension must be an integer scalar"); endif nd = ndims (x); sz = size (x); sz2 = floor (sz(dim) / 2); idx = repmat ({':'}, nd, 1); idx{dim} = [sz2+1:sz(dim), 1:sz2]; retval = x(idx{:}); else if (isvector (x)) xl = length (x); xx = floor (xl/2); retval = x([xx+1:xl, 1:xx]); elseif (ismatrix (x)) nd = ndims (x); sz = size (x); sz2 = floor (sz ./ 2); idx = cell (); for i = 1:nd idx{i} = [sz2(i)+1:sz(i), 1:sz2(i)]; endfor retval = x(idx{:}); else error ("ifftshift: expecting vector or matrix argument"); endif endif endfunction %!test %! x = [0:7]; %! y = ifftshift (x); %! assert (y, [4 5 6 7 0 1 2 3]); %! assert (ifftshift (y), x); %!test %! x = [0:6]; %! y = ifftshift (x); %! assert (y, [3 4 5 6 0 1 2]); %! assert (ifftshift (y), [6 0 1 2 3 4 5]); %!test %! x = [0:7]'; %! y = ifftshift (x); %! assert (y, [4;5;6;7;0;1;2;3]); %! assert (ifftshift (y), x); %!test %! x = [0:6]'; %! y = ifftshift (x); %! assert (y, [3;4;5;6;0;1;2]); %! assert (ifftshift (y), [6;0;1;2;3;4;5]); %!test %! x = [0:3]; %! x = [x;2*x;3*x+1;4*x+1]; %! y = ifftshift (x); %! assert (y, [[7 10 1 4];[9 13 1 5];[2 3 0 1];[4 6 0 2]]); %! assert (ifftshift (y), x); %!test %! x = [0:3]; %! x = [x;2*x;3*x+1;4*x+1]; %! y = ifftshift (x,1); %! assert (y, [[1 4 7 10];[1 5 9 13];[0 1 2 3];[0 2 4 6]]); %! assert (ifftshift (y,1), x); %!test %! x = [0:3]; %! x = [x;2*x;3*x+1;4*x+1]; %! y = ifftshift (x,2); %! assert (y, [[2 3 0 1];[4 6 0 2];[7 10 1 4];[9 13 1 5]]); %! assert (ifftshift (y,2), x); usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/periodogram.m000644 000765 000120 00000011175 12402246361 026412 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## Copyright (C) 2010 Alois Schloegl ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[Pxx, @var{w}] =} periodogram (@var{x}) ## For a data matrix @var{x} from a sample of size @var{n}, return the ## periodogram. The angular frequency is returned in @var{w}. ## ## [Pxx,w] = periodogram (@var{x}). ## ## [Pxx,w] = periodogram (@var{x},win). ## ## [Pxx,w] = periodogram (@var{x},win,nfft). ## ## [Pxx,f] = periodogram (@var{x},win,nfft,Fs). ## ## [Pxx,f] = periodogram (@var{x},win,nfft,Fs,"range"). ## ## @itemize ## @item x: data; if real-valued a one-sided spectrum is estimated, ## if complex-valued or range indicates @qcode{"@nospell{twosided}"}, the full ## spectrum is estimated. ## ## @item win: weight data with window, x.*win is used for further computation, ## if window is empty, a rectangular window is used. ## ## @item nfft: number of frequency bins, default max (256, 2.^ceil (log2 (length (x)))). ## ## @item Fs: sampling rate, default 1. ## ## @item range: @qcode{"@nospell{onesided}"} computes spectrum from [0..nfft/2+1]. ## @qcode{"@nospell{twosided}"} computes spectrum from [0..nfft-1]. These ## strings can appear at any position in the list input arguments after ## window. ## ## @item @nospell{Pxx}: one-, or two-sided power spectrum. ## ## @item w: angular frequency [0..2*pi) (two-sided) or [0..pi] one-sided. ## ## @item f: frequency [0..Fs) (two-sided) or [0..Fs/2] one-sided. ## @end itemize ## @end deftypefn ## Author: FL ## Description: Compute the periodogram function [pxx, f] = periodogram (x, varargin) ## check input arguments if (nargin < 1 || nargin > 5) print_usage (); endif nfft = []; fs = []; range = []; window = []; j = 1; for k = 1:length (varargin) if (ischar (varargin{k})) range = varargin{k}; else switch (j) case 1 window = varargin{k}; case 2 nfft = varargin{k}; case 3 fs = varargin{k}; case 4 range = varargin{k}; endswitch j++; endif endfor [r, c] = size (x); if (r == 1) r = c; endif if (ischar (window)) range = window; window = []; endif; if (ischar (nfft)) range = nfft; nfft = []; endif; if (ischar (fs)) range = fs; fs = []; endif; if (! isempty (window)) if (all (size (x) == size (window))) x .*= window; elseif (rows (x) == rows (window) && columns (window) == 1) x .*= window (:,ones (1,c)); endif; endif if (numel (nfft)>1) error ("nfft must be scalar"); endif if (isempty (nfft)) nfft = max (256, 2.^ceil (log2 (r))); endif if (strcmp (range, "onesided")) range = 1; elseif (strcmp (range, "twosided")) range = 2; else range = 2-isreal (x); endif ## compute periodogram if (r>nfft) Pxx = 0; rr = rem (length (x), nfft); if (rr) x = [x(:); (zeros (nfft-rr, 1))]; endif x = sum (reshape (x, nfft, []), 2); endif if (isempty (window)) n = r; else n = sumsq (window); end; Pxx = (abs (fft (x, nfft))) .^ 2 / n ; if (nargin<4) Pxx /= 2*pi; elseif (! isempty (fs)) Pxx /= fs; endif ## generate output arguments if (range == 1) # onesided Pxx = Pxx(1:nfft/2+1) + [0; Pxx(end:-1:(nfft/2+2)); 0]; endif if (nargout != 1) if (range == 1) f = (0:nfft/2)'/nfft; elseif (range == 2) f = (0:nfft-1)'/nfft; endif if (nargin<4) f *= 2*pi; # generate w=2*pi*f elseif (! isempty (fs)) f *= fs; endif endif if (nargout == 0) if (nargin<4) plot (f/(2*pi), 10*log10 (Pxx)); xlabel ("normalized frequency [x pi rad]"); ylabel ("Power density [dB/rad/sample]"); else plot (f, 10*log10 (Pxx)); xlabel ("frequency [Hz]"); ylabel ("Power density [dB/Hz]"); endif grid on; title ("Periodogram Power Spectral Density Estimate"); else pxx = Pxx; endif endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/private/000755 000765 000120 00000000000 12402246361 025371 5ustar00jburnettadmin000000 000000 usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/sinc.m000644 000765 000120 00000002362 12402246361 025034 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1994-2013 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} sinc (@var{x}) ## Return ## @tex ## $\sin (\pi x)/(\pi x)$. ## @end tex ## @ifnottex ## sin (pi*x) / (pi*x). ## @end ifnottex ## @end deftypefn ## Author: jwe ??? function result = sinc (x) if (nargin != 1) print_usage (); endif result = ones (size (x)); i = (x != 0); if (any (i(:))) t = pi * x(i); result(i) = sin (t) ./ t; endif endfunction %!assert (sinc (0), 1) %!assert (sinc (1), 0,1e-6) %!assert (sinc (1/2), 2/pi, 1e-6) %!error sinc() usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/sinetone.m000644 000765 000120 00000003760 12402246361 025727 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} sinetone (@var{freq}, @var{rate}, @var{sec}, @var{ampl}) ## Return a sinetone of frequency @var{freq} with length of @var{sec} ## seconds at sampling rate @var{rate} and with amplitude @var{ampl}. ## The arguments @var{freq} and @var{ampl} may be vectors of common size. ## ## Defaults are @var{rate} = 8000, @var{sec} = 1 and @var{ampl} = 64. ## @end deftypefn ## Author: FL ## Description: Compute a sine tone function retval = sinetone (freq, rate, sec, ampl) if (nargin == 1) rate = 8000; sec = 1; ampl = 64; elseif (nargin == 2) sec = 1; ampl = 64; elseif (nargin == 3) ampl = 64; elseif ((nargin < 1) || (nargin > 4)) print_usage (); endif [err, freq, ampl] = common_size (freq, ampl); if (err || ! isvector (freq)) error ("sinetone: FREQ and AMPL must be vectors of common size"); endif if (! (isscalar (rate) && isscalar (sec))) error ("sinetone: RATE and SEC must be scalars"); endif n = length (freq); ns = round (rate * sec); retval = zeros (ns, n); for k = 1:n retval(:, k) = ampl(k) * sin (2 * pi * (1:ns) / rate * freq(k))'; endfor endfunction %!assert (size (sinetone (18e6, 150e6, 19550/150e6, 1)), [19550, 1]) usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/sinewave.m000644 000765 000120 00000003260 12402246361 025717 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Andreas Weingessel ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} sinewave (@var{m}, @var{n}, @var{d}) ## Return an @var{m}-element vector with @var{i}-th element given by ## @code{sin (2 * pi * (@var{i}+@var{d}-1) / @var{n})}. ## ## The default value for @var{d} is 0 and the default value for @var{n} ## is @var{m}. ## @end deftypefn ## Author: AW ## Description: Compute a sine wave function x = sinewave (m, n, d) if (nargin > 0 && nargin < 4) if (nargin < 3) d = 0; endif if (nargin < 2) n = m; endif x = sin (((1 : m) + d - 1) * 2 * pi / n); else print_usage (); endif endfunction %!assert (sinewave (1), 0) %!assert (sinewave (1, 4, 1), 1) %!assert (sinewave (1, 12, 1), 1/2, 1e-6) %!assert (sinewave (1, 12, 2), sqrt (3)/2, 1e-6) %!assert (sinewave (1, 20, 1), (sqrt (5)-1)/4, 1e-6) %!assert (sinewave (1), sinewave (1, 1,0)) %!assert (sinewave (3, 4), sinewave (3, 4, 0)) %!error sinewave () usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/spectral_adf.m000644 000765 000120 00000004254 12402246361 026531 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} spectral_adf (@var{c}) ## @deftypefnx {Function File} {} spectral_adf (@var{c}, @var{win}) ## @deftypefnx {Function File} {} spectral_adf (@var{c}, @var{win}, @var{b}) ## Return the spectral density estimator given a vector of ## autocovariances @var{c}, window name @var{win}, and bandwidth, ## @var{b}. ## ## The window name, e.g., @qcode{"triangle"} or @qcode{"rectangle"} is ## used to search for a function called @code{@var{win}_lw}. ## ## If @var{win} is omitted, the triangle window is used. If @var{b} is ## omitted, @code{1 / sqrt (length (@var{x}))} is used. ## @seealso{spectral_xdf} ## @end deftypefn ## Author: FL ## Description: Spectral density estimation function retval = spectral_adf (c, win, b) if (nargin < 1 || nargin > 3) print_usage (); endif cr = length (c); if (columns (c) > 1) c = c'; endif if (nargin < 3) b = 1 / ceil (sqrt (cr)); endif if (nargin == 1) w = triangle_lw (cr, b); elseif (! ischar (win)) error ("spectral_adf: WIN must be a string"); else win = str2func ([win "_lw"]); w = feval (win, cr, b); endif c = c .* w; retval = 2 * real (fft (c)) - c(1); retval = [(zeros (cr, 1)), retval]; retval(:, 1) = (0 : cr-1)' / cr; endfunction %% Test input validation %!error spectral_adf (); %!error spectral_adf (1, 2, 3, 4); %!error spectral_adf (1, 2); %!error spectral_adf (1, "invalid"); usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/spectral_xdf.m000644 000765 000120 00000004326 12402246361 026560 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} spectral_xdf (@var{x}) ## @deftypefnx {Function File} {} spectral_xdf (@var{x}, @var{win}) ## @deftypefnx {Function File} {} spectral_xdf (@var{x}, @var{win}, @var{b}) ## Return the spectral density estimator given a data vector @var{x}, ## window name @var{win}, and bandwidth, @var{b}. ## ## The window name, e.g., @qcode{"triangle"} or @qcode{"rectangle"} is ## used to search for a function called @code{@var{win}_sw}. ## ## If @var{win} is omitted, the triangle window is used. If @var{b} is ## omitted, @code{1 / sqrt (length (@var{x}))} is used. ## @seealso{spectral_adf} ## @end deftypefn ## Author: FL ## Description: Spectral density estimation function retval = spectral_xdf (x, win, b) if (nargin < 1 || nargin > 3) print_usage (); endif xr = length (x); if (columns (x) > 1) x = x'; endif if (nargin < 3) b = 1 / ceil (sqrt (xr)); endif if (nargin == 1) w = triangle_sw (xr, b); elseif (! ischar (win)) error ("spectral_xdf: WIN must be a string"); else win = str2func ([win "_sw"]); w = feval (win, xr, b); endif x = x - sum (x) / xr; retval = (abs (fft (x)) / xr).^2; retval = real (ifft (fft (retval) .* fft (w))); retval = [(zeros (xr, 1)), retval]; retval(:, 1) = (0 : xr-1)' / xr; endfunction %% Test input validation %!error spectral_xdf (); %!error spectral_xdf (1, 2, 3, 4); %!error spectral_xdf (1, 2); %!error spectral_xdf (1, "invalid"); usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/spencer.m000644 000765 000120 00000002606 12402246361 025540 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} spencer (@var{x}) ## Return Spencer's 15 point moving average of each column of ## @var{x}. ## @end deftypefn ## Author: FL ## Description: Apply Spencer's 15-point MA filter function retval = spencer (x) if (nargin != 1) print_usage (); endif [xr, xc] = size (x); n = xr; c = xc; if (isvector (x)) n = length (x); c = 1; x = reshape (x, n, 1); endif w = [-3, -6, -5, 3, 21, 46, 67, 74, 67, 46, 21, 3, -5, -6, -3] / 320; retval = fftfilt (w, x); retval = [zeros(7,c); retval(15:n,:); zeros(7,c);]; retval = reshape (retval, xr, xc); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/stft.m000644 000765 000120 00000007304 12402246361 025061 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Andreas Weingessel ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {@var{y} =} stft (@var{x}) ## @deftypefnx {Function File} {@var{y} =} stft (@var{x}, @var{win_size}) ## @deftypefnx {Function File} {@var{y} =} stft (@var{x}, @var{win_size}, @var{inc}) ## @deftypefnx {Function File} {@var{y} =} stft (@var{x}, @var{win_size}, @var{inc}, @var{num_coef}) ## @deftypefnx {Function File} {@var{y} =} stft (@var{x}, @var{win_size}, @var{inc}, @var{num_coef}, @var{win_type}) ## @deftypefnx {Function File} {[@var{y}, @var{c}] =} stft (@dots{}) ## Compute the short-time Fourier transform of the vector @var{x} with ## @var{num_coef} coefficients by applying a window of @var{win_size} data ## points and an increment of @var{inc} points. ## ## Before computing the Fourier transform, one of the following windows ## is applied: ## ## @table @asis ## @item @qcode{"hanning"} ## win_type = 1 ## ## @item @qcode{"hamming"} ## win_type = 2 ## ## @item @qcode{"rectangle"} ## win_type = 3 ## @end table ## ## The window names can be passed as strings or by the @var{win_type} number. ## ## The following defaults are used for unspecified arguments: ## @var{win_size} = 80, @var{inc} = 24, @var{num_coef} = 64, and ## @var{win_type} = 1. ## ## @code{@var{y} = stft (@var{x}, @dots{})} returns the absolute values ## of the Fourier coefficients according to the @var{num_coef} positive ## frequencies. ## ## @code{[@var{y}, @var{c}] = stft (@code{x}, @dots{})} returns the ## entire STFT-matrix @var{y} and a 3-element vector @var{c} containing ## the window size, increment, and window type, which is needed by the ## @code{synthesis} function. ## @seealso{synthesis} ## @end deftypefn ## Author: AW ## Description: Short-Time Fourier Transform function [y, c] = stft (x, win_size = 80, inc = 24, num_coef = 64, win_type = 1) if (nargin < 1 || nargin > 5) print_usage (); endif if (ischar (win_type)) switch (tolower (win_type)) case "hanning" win_type = 1; case "hamming" win_type = 2; case "rectangle" win_type = 3; otherwise error ("stft: unknown window type '%s'", win_type); endswitch endif ## Check whether X is a vector. if (! isvector (x)) error ("stft: X must be a vector"); endif x = x(:); ncoef = 2 * num_coef; if (win_size > ncoef) win_size = ncoef; printf ("stft: window size adjusted to %f\n", win_size); endif num_win = fix ((rows (x) - win_size) / inc); ## compute the window coefficients switch (win_type) case 1 win_coef = hanning (win_size); case 2 win_coef = hamming (win_size); case 3 win_coef = ones (win_size, 1); endswitch ## Create a matrix Z whose columns contain the windowed time-slices. z = zeros (ncoef, num_win + 1); start = 1; for i = 0:num_win z(1:win_size, i+1) = x(start:start+win_size-1) .* win_coef; start = start + inc; endfor y = fft (z); if (nargout == 1) y = abs (y(1:num_coef, :)); else c = [win_size, inc, win_type]; endif endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/synthesis.m000644 000765 000120 00000003654 12402246361 026136 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Andreas Weingessel ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {@var{x} =} synthesis (@var{y}, @var{c}) ## Compute a signal from its short-time Fourier transform @var{y} and a ## 3-element vector @var{c} specifying window size, increment, and ## window type. ## ## The values @var{y} and @var{c} can be derived by ## ## @example ## [@var{y}, @var{c}] = stft (@var{x} , @dots{}) ## @end example ## @seealso{stft} ## @end deftypefn ## Author: AW ## Description: Recover a signal from its short-term Fourier transform function x = synthesis (y, c) if (nargin != 2) print_usage (); endif if (numel (c) != 3) error ("synthesis: C must contain exactly 3 elements"); endif w_size = c(1); inc = c(2); w_type = c(3); if (w_type == 1) w_coeff = hanning (w_size); elseif (w_type == 2) w_coeff = hamming (w_size); elseif (w_type == 3) w_coeff = ones (w_size, 1); else error ("synthesis: window_type must be 1, 2, or 3"); endif z = real (ifft (y)); st = fix ((w_size-inc) / 2); z = z(st+1:st+inc, :); w_coeff = w_coeff(st+1:st+inc); nc = columns (z); for i = 1:nc z(:, i) ./= w_coeff; endfor x = reshape (z, inc * nc, 1); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/unwrap.m000644 000765 000120 00000012012 12402246361 025405 0ustar00jburnettadmin000000 000000 ## Copyright (C) 2000-2013 Bill Lash ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {@var{b} =} unwrap (@var{x}) ## @deftypefnx {Function File} {@var{b} =} unwrap (@var{x}, @var{tol}) ## @deftypefnx {Function File} {@var{b} =} unwrap (@var{x}, @var{tol}, @var{dim}) ## ## Unwrap radian phases by adding multiples of 2*pi as appropriate to ## remove jumps greater than @var{tol}. @var{tol} defaults to pi. ## ## Unwrap will work along the dimension @var{dim}. If @var{dim} ## is unspecified it defaults to the first non-singleton dimension. ## @end deftypefn ## Author: Bill Lash function retval = unwrap (x, tol, dim) if (nargin < 1 || nargin > 3) print_usage (); endif if (!isnumeric (x)) error ("unwrap: X must be a numeric matrix or vector"); endif if (nargin < 2 || isempty (tol)) tol = pi; endif ## Don't let anyone use a negative value for TOL. tol = abs (tol); nd = ndims (x); sz = size (x); if (nargin == 3) if (!(isscalar (dim) && dim == fix (dim)) || !(1 <= dim && dim <= nd)) error ("unwrap: DIM must be an integer and a valid dimension"); endif else ## Find the first non-singleton dimension. (dim = find (sz > 1, 1)) || (dim = 1); endif rng = 2*pi; m = sz(dim); ## Handle case where we are trying to unwrap a scalar, or only have ## one sample in the specified dimension. if (m == 1) retval = x; return; endif ## Take first order difference to see so that wraps will show up ## as large values, and the sign will show direction. idx = repmat ({':'}, nd, 1); idx{dim} = [1,1:m-1]; d = x(idx{:}) - x; ## Find only the peaks, and multiply them by the appropriate amount ## of ranges so that there are kronecker deltas at each wrap point ## multiplied by the appropriate amount of range values. p = ceil (abs (d)./rng) .* rng .* (((d > tol) > 0) - ((d < -tol) > 0)); ## Now need to "integrate" this so that the deltas become steps. r = cumsum (p, dim); ## Now add the "steps" to the original data and put output in the ## same shape as originally. retval = x + r; endfunction %!function t = __xassert (a,b,tol) %! if (nargin == 1) %! t = all (a(:)); %! else %! if (nargin == 2) %! tol = 0; %! endif %! if (any (size (a) != size (b))) %! t = 0; %! elseif (any (abs (a(:) - b(:)) > tol)) %! t = 0; %! else %! t = 1; %! endif %! endif %!endfunction %! %!test %! %! i = 0; %! t = []; %! %! r = [0:100]; # original vector %! w = r - 2*pi*floor ((r+pi)/(2*pi)); # wrapped into [-pi,pi] %! tol = 1e3*eps; # maximum expected deviation %! %! t(++i) = __xassert (r, unwrap (w), tol); #unwrap single row %! t(++i) = __xassert (r', unwrap (w'), tol); #unwrap single column %! t(++i) = __xassert ([r',r'], unwrap ([w',w']), tol); #unwrap 2 columns %! t(++i) = __xassert ([r;r], unwrap ([w;w],[],2), tol); #check that dim works %! t(++i) = __xassert (r+10, unwrap (10+w), tol); #check r(1)>pi works %! %! t(++i) = __xassert (w', unwrap (w',[],2)); #unwrap col by rows should not change it %! t(++i) = __xassert (w, unwrap (w,[],1)); #unwrap row by cols should not change it %! t(++i) = __xassert ([w;w], unwrap ([w;w])); #unwrap 2 rows by cols should not change them %! %! ## verify that setting tolerance too low will cause bad results. %! t(++i) = __xassert (any (abs (r - unwrap (w,0.8)) > 100)); %! %! assert (all (t)); %! %!test %! A = [pi*(-4), pi*(-2+1/6), pi/4, pi*(2+1/3), pi*(4+1/2), pi*(8+2/3), pi*(16+1), pi*(32+3/2), pi*64]; %! assert (unwrap (A), unwrap (A, pi)); %! assert (unwrap (A, pi), unwrap (A, pi, 2)); %! assert (unwrap (A', pi), unwrap (A', pi, 1)); %! %!test %! A = [pi*(-4); pi*(2+1/3); pi*(16+1)]; %! B = [pi*(-2+1/6); pi*(4+1/2); pi*(32+3/2)]; %! C = [pi/4; pi*(8+2/3); pi*64]; %! D = [pi*(-2+1/6); pi*(2+1/3); pi*(8+2/3)]; %! E(:, :, 1) = [A, B, C, D]; %! E(:, :, 2) = [A+B, B+C, C+D, D+A]; %! F(:, :, 1) = [unwrap(A), unwrap(B), unwrap(C), unwrap(D)]; %! F(:, :, 2) = [unwrap(A+B), unwrap(B+C), unwrap(C+D), unwrap(D+A)]; %! assert (unwrap (E), F); %! %!test %! A = [0, 2*pi, 4*pi, 8*pi, 16*pi, 65536*pi]; %! B = [pi*(-2+1/6), pi/4, pi*(2+1/3), pi*(4+1/2), pi*(8+2/3), pi*(16+1), pi*(32+3/2), pi*64]; %! assert (unwrap (A), zeros (1, length (A))); %! assert (diff (unwrap (B), 1) < 2*pi, true (1, length (B)-1)); %!error unwrap() usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/yulewalker.m000644 000765 000120 00000002741 12402246361 026265 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{a}, @var{v}] =} yulewalker (@var{c}) ## Fit an AR (p)-model with Yule-Walker estimates given a vector @var{c} ## of autocovariances @code{[gamma_0, @dots{}, gamma_p]}. ## ## Returns the AR coefficients, @var{a}, and the variance of white ## noise, @var{v}. ## @end deftypefn ## Author: FL ## Description: Fit AR model by Yule-Walker method function [a, v] = yulewalker (c) if (nargin != 1) print_usage (); endif p = length (c) - 1; if (columns (c) > 1) c = c'; endif cp = c(2 : p+1); CP = zeros (p, p); for i = 1:p for j = 1:p CP (i, j) = c (abs (i-j) + 1); endfor endfor a = inv (CP) * cp; v = c(1) - a' * cp; endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/private/rectangle_lw.m000644 000765 000120 00000002165 12402246361 030221 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} rectangle_lw (@var{n}, @var{b}) ## Rectangular lag window. Subfunction used for spectral density ## estimation. ## @end deftypefn ## Author: FL ## Description: Rectangular lag window function retval = rectangle_lw (n, b) retval = zeros (n, 1); t = floor (1 / b); retval(1:t, 1) = ones (t, 1); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/private/rectangle_sw.m000644 000765 000120 00000002304 12402246361 030223 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} rectangle_sw (@var{n}, @var{b}) ## Rectangular spectral window. Subfunction used for spectral density ## estimation. ## @end deftypefn ## Author: FL ## Description: Rectangular spectral window function retval = rectangle_sw (n, b) retval = zeros (n, 1); retval(1) = 2 / b + 1; l = (2:n)' - 1; l = 2 * pi * l / n; retval(2:n) = sin ((2/b + 1) * l / 2) ./ sin (l / 2); endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/private/triangle_lw.m000644 000765 000120 00000002156 12402246361 030062 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} triangle_lw (@var{n}, @var{b}) ## Triangular lag window. Subfunction used for spectral density ## estimation. ## @end deftypefn ## Author: FL ## Description: Triangular lag window function retval = triangle_lw (n, b) retval = 1 - (0 : n-1)' * b; retval = max ([retval'; (zeros (1, n))])'; endfunction usr/local/Cellar/octave//3.8.1/share/octave/3.8.1/m/signal/private/triangle_sw.m000644 000765 000120 00000002274 12402246361 030072 0ustar00jburnettadmin000000 000000 ## Copyright (C) 1995-2013 Friedrich Leisch ## ## This file is part of Octave. ## ## Octave 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. ## ## Octave 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 Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {} triangle_sw (@var{n}, @var{b}) ## Triangular spectral window. Subfunction used for spectral density ## estimation. ## @end deftypefn ## Author: FL ## Description: Triangular spectral window function retval = triangle_sw (n, b) retval = zeros (n,1); retval(1) = 1 / b; l = (2:n)' - 1; l = 2 * pi * l / n; retval(2:n) = b * (sin (l / (2*b)) ./ sin (l / 2)).^2; endfunction