|
samer@32
|
1 function [S,F]=powspecd(X,N)
|
|
samer@32
|
2 % powspecd - power spectral density of x
|
|
samer@32
|
3 %
|
|
samer@32
|
4 % powspecd ::
|
|
samer@32
|
5 % [[K,L]] ~'sequence of L K-point signal frames',
|
|
samer@32
|
6 % N:natural ~'optional size of FFT (defaults to K)'
|
|
samer@32
|
7 % -> [[M,L]] ~'sequence of L M-point power spectra, M=dftbins(N)',
|
|
samer@32
|
8 % [[M]] ~'normalised bandwidths for each spectral bin'.
|
|
samer@32
|
9 %
|
|
samer@32
|
10 % This function is careful to account for the power in the original
|
|
samer@32
|
11 % signal properly: we return the total power density in the signal in
|
|
samer@32
|
12 % each band of the DFT such that:
|
|
samer@32
|
13 % sum(powspecd(x,N).*dftbw(N)) = mean(x.^2)
|
|
samer@32
|
14 % dftbw(N) is the bandwidth of the N'th bin in the spectrum, in
|
|
samer@32
|
15 % units of fs/N, where fs is the original sampling frequency.
|
|
samer@32
|
16 % When N is odd, dftbw(N) =[0.5, 1, ..., 1] but when N is
|
|
samer@32
|
17 % even, dftbw(N) = [0.5, 1, ..., 1, 0.5].
|
|
samer@32
|
18 %
|
|
samer@32
|
19 % The second return, if requested, is dftbw(N), so if
|
|
samer@32
|
20 % [bw,y]=powspecd(x)
|
|
samer@32
|
21 % then
|
|
samer@32
|
22 % sum(y.*bw) = mean(x.^2);
|
|
samer@32
|
23
|
|
samer@32
|
24 T=size(X,1);
|
|
samer@32
|
25 if nargin<2, N=T; end
|
|
samer@32
|
26
|
|
samer@32
|
27 select=1:(1+floor(N/2)); % non-redundant bands
|
|
samer@32
|
28 Y=fft(X,N);
|
|
samer@32
|
29 S=(2/(N*T))*abs(Y(select,:)).^2; % power in selected bands
|
|
samer@32
|
30 if nargout>1, F=dftbw(N); end % return bandwidths
|
|
samer@32
|
31
|
