view dsp/powspecd.m @ 61:eff6bddf82e3 tip

Finally implemented perceptual brightness thing.

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