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