view dsp/powspec.m @ 61:eff6bddf82e3 tip

Finally implemented perceptual brightness thing.

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