diff dsp/powspec.m @ 32:c3b0cd708782

Imported core dsp tools.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/powspec.m	Sun Jan 20 13:48:47 2013 +0000
@@ -0,0 +1,26 @@
+function [S,F]=powspec(X,N)
+% powspec - power spectra of columns of x
+%
+% powspec ::
+%    [[K,L]]    ~'sequence of L K-point signal frames',
+%    N:natural  ~'optional size of FFT (defaults to K)'
+% -> [[M,L]]    ~'sequence of L M-point power spectra, M=dftbins(N)',
+%    [[M]]      ~'normalised centre frequencies for each spectral bin'.
+%
+% This function is careful to account for the energy in the original
+% signal properly: we return the total energy in the signal in each 
+% band of the DFT such that sum(powspec(x))/4 = sum(x.^2).
+% Don't ask about the 4 - I don't know why it's in there. It just is
+% NB. this is not strictly the power spectral *density* because
+% when N is even, first and last bands, the singleton bands, are
+% half as wide as the others.
+
+if nargin<2, N=size(X,1); end
+
+select=1:(1+floor(N/2));         % non-redundant bands
+singles=1:ceil(N/2):select(end); % bands which don't come in pairs
+Y=fft(X,N);
+S=(8/N)*abs(Y(select,:)).^2;     % power in selected bands
+S(singles,:)=S(singles,:)/2;     % undo double counting of singleton bands
+if nargout>1, F=(select-1)/N; end % normalised centre frequencies
+