view m-toolbox/classes/@ao/mlpsd_m.m @ 0:f0afece42f48
Import.
author |
Daniele Nicolodi <nicolodi@science.unitn.it> |
date |
Wed, 23 Nov 2011 19:22:13 +0100 (2011-11-23) |
parents |
|
children |
|
line source
% MLPSD_M m-file only version of the LPSD algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: MLPSD_M m-file only version of the LPSD algorithm
%
% CALL: [S,Sxx,ENBW] = mlpsd_m(x,f,r,m,L,fs,win,order,olap)
%
% VERSION: $Id: mlpsd_m.m,v 1.9 2011/05/22 22:22:50 mauro Exp $
%
% HISTORY: 19-05-2007 M Hewitson
% Creation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = mlpsd_m(varargin)
import utils.const.*
utils.helper.msg(msg.PROC1, 'using MATLAB to compute core DFT');
% Get inputs
x = varargin{1};
f = varargin{2};
r = varargin{3};
m = varargin{4};
L = varargin{5};
fs = varargin{6};
win = varargin{7};
order = varargin{8};
olap = varargin{9};
twopi = 2.0*pi;
nx = length(x);
nf = length(f);
Sxx = zeros(nf,1);
S = zeros(nf,1);
ENBW = zeros(nf,1);
disp_each = round(nf/100)*10;
for jj = 1:nf
if mod(jj, disp_each) == 0 || jj == 1 || jj == nf
utils.helper.msg(msg.PROC1, 'computing frequency %04d of %d: %f Hz', jj, nf, f(jj));
end
% compute DFT exponent and window
l = L(jj);
switch lower(win.type)
case 'kaiser'
win = specwin(win.type, l, win.psll);
otherwise
win = specwin(win.type, l);
end
p = 1i * twopi * m(jj)/l.*[0:l-1];
C = win.win .* exp(p);
Xolap = (1-olap);
% do segments
A = 0.0;
% Compute the number of averages we want here
segLen = l;
nData = length(x);
ovfact = 1 / (1 - olap);
davg = (((nData - segLen)) * ovfact) / segLen + 1;
navg = round(davg);
% Compute steps between segments
if navg == 1
shift = 1;
else
shift = (nData - segLen) / (navg - 1);
end
if shift < 1 || isnan(shift)
shift = 1;
end
% disp(sprintf('Seglen: %d\t | Shift: %f\t | navs: %d', segLen, shift, navg))
start = 1.0;
for ii = 1:navg
% compute start index
istart = round (start);
start = start + shift;
% get segment
xs = [x(istart:istart+l-1)].';
% detrend segment
switch order
case -1
% do nothing
case 0
xs = xs - mean(xs);
case 1
xs = detrend(xs);
otherwise
xs = polydetrend(xs, order);
end
% make DFT
a = C*xs;
A = A + a*conj(a);
end
if mod(jj, disp_each) == 0 || jj == 1 || jj == nf
utils.helper.msg(msg.PROC2, 'averaged %d segments', navg);
end
A2ns = 2.0*A/navg;
S1 = win.ws;
S12 = S1*S1;
S2 = win.ws2;
ENBW(jj) = fs*S2/S12;
Sxx(jj) = A2ns/fs/S2;
S(jj) = A2ns/S12;
end % for j=1:nf
varargout{1} = S;
varargout{2} = Sxx;
varargout{3} = ENBW;
end