help-octave
[Top][All Lists]

## 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

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