samer@32: function [S,F]=powspecd(X,N)
samer@32: % powspecd - power spectral density of x
samer@32: %
samer@32: % powspecd ::
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 bandwidths for each spectral bin'.
samer@32: %
samer@32: % This function is careful to account for the power in the original
samer@32: % signal properly: we return the total power density in the signal in
samer@32: % each band of the DFT such that:
samer@32: %     sum(powspecd(x,N).*dftbw(N)) = mean(x.^2)
samer@32: % dftbw(N) is the bandwidth of the N'th bin in the spectrum, in
samer@32: % units of fs/N, where fs is the original sampling frequency.
samer@32: % When N is odd, dftbw(N) =[0.5, 1, ..., 1] but when N is
samer@32: % even, dftbw(N) = [0.5, 1, ..., 1, 0.5].
samer@32: %
samer@32: % The second return, if requested, is dftbw(N), so if
samer@32: %    [bw,y]=powspecd(x)
samer@32: % then
samer@32: %    sum(y.*bw) = mean(x.^2);
samer@32: 
samer@32: T=size(X,1);
samer@32: if nargin<2, N=T; end
samer@32: 
samer@32: select=1:(1+floor(N/2));         % non-redundant bands
samer@32: Y=fft(X,N);
samer@32: S=(2/(N*T))*abs(Y(select,:)).^2; % power in selected bands
samer@32: if nargout>1, F=dftbw(N); end % return bandwidths
samer@32: