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