diff dsp/powspecd.m @ 32:c3b0cd708782

Imported core dsp tools.

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