help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Spectrum Estimate functions


From: Francesco Potorti`
Subject: Re: Spectrum Estimate functions
Date: Fri, 9 May 2003 07:49:14 -0500

I have used this for some time.  I don't know if it does anything more
or better than the current spectrum evaluation functions.


===File ~/math/octavelib.old/signal/spectrum.m==============
function A = spectrum(x, n)

  ## usage: A = spectrum (X, N)
  ##
  ## Given a vector X in the time domain, compute its one-sided power
  ## spectrum in N frequency points plus the null frequency by repeatedly
  ## applying the fft to overlapping chunks of 2*N points. Faster if N is a
  ## power of 2.  Gives a vector of length N+1.  Changing N changes the
  ## frequency spacing (that is the frequency resolution), that is, it
  ## changes the minimum frequency computed, not the maximum one.
  ##
  ## The mean value of each element of the resulting vector is equal to the
  ## average power computed on a frequency window centered on the element
  ## itself. Its variance is equal to (11/9)*(N/length(X)) * mean^2.
  ##
  ## Algorithm inspired by Numerical recipes in C, II edition, pg.550ff.
  ## Adapted to Octave by Francesco Potortì <address@hidden>, 1997
  ## License: GPL <http://www.gnu.org/licenses/gpl.html>

  if (nargin != 2)
    usage ("A = spectrum (x, n)");
  endif

  if (!is_scalar(n))
    error ("spectrum: n must be a scalar");
  endif

  if (! is_vector (x))
    error ("lin2mu: x must be a vector");
  endif

  column_vector = (columns(x) == 1);

  if (length(x) < 2*n)
    x = postpad (x, 2*n);
  endif
  chunks = floor(length(x)/n)-1;
  x = reshape (x(1:n*(chunks+1)), n, chunks+1);
  a = [ x(:,1:chunks); x(:,2:chunks+1) ];
  A = fft(a).';
  A = sumsq(abs(A));
  A = A(1:n+1) + [0, fliplr(A(n+1:2*n))];
  A = A / (2*n)^2;

  if (column_vector)
    A = A';
  endif
endfunction
============================================================



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------



reply via email to

[Prev in Thread] Current Thread [Next in Thread]