Home > classes > @ao > mlpsd_m.m

mlpsd_m

PURPOSE ^

MLPSD_M m-file only version of the LPSD algorithm

SYNOPSIS ^

function varargout = mlpsd_m(varargin)

DESCRIPTION ^

 MLPSD_M m-file only version of the LPSD algorithm

 [S,Sxx,ENBW] = mlpsd_m(x,f,r,m,L,fs,win,order,olap)

 M Hewitson 19-05-07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % MLPSD_M m-file only version of the LPSD algorithm
0002 %
0003 % [S,Sxx,ENBW] = mlpsd_m(x,f,r,m,L,fs,win,order,olap)
0004 %
0005 % M Hewitson 19-05-07
0006 %
0007 %
0008 
0009 function varargout = mlpsd_m(varargin)
0010 
0011   import utils.const.*
0012   
0013   utils.helper.msg(msg.PROC1, 'using MATLAB to compute core DFT');
0014   
0015   % Get inputs
0016   x     = varargin{1};
0017   f     = varargin{2};
0018   r     = varargin{3};
0019   m     = varargin{4};
0020   L     = varargin{5};
0021   fs    = varargin{6};
0022   win   = varargin{7};
0023   order = varargin{8};
0024   olap  = varargin{9};
0025 
0026   twopi    = 2.0*pi;
0027 
0028   nx   = length(x);
0029   nf   = length(f);
0030   Sxx  = zeros(nf,1);
0031   S    = zeros(nf,1);
0032   ENBW = zeros(nf,1);
0033 
0034   disp_each = round(nf/100)*10;
0035 
0036   for j=1:nf
0037 
0038     if mod(j,disp_each) == 0 || j == 1 || j == nf
0039       utils.helper.msg(msg.PROC1, 'computing frequency %04d of %d: %f Hz', j, nf, f(j));
0040     end
0041 
0042     % compute DFT exponent and window
0043     l = L(j);
0044     switch win.type
0045       case 'Kaiser'
0046         win = specwin(win.type, l, win.psll);
0047       otherwise
0048         win = specwin(win.type, l);
0049     end
0050 
0051     p     = 1i * twopi * m(j)/l.*[0:l-1];
0052     C     = win.win .* exp(p);
0053     Xolap = (1-olap);
0054     % do segments
0055     A  = 0.0;
0056 
0057     % Compute the number of averages we want here
0058     segLen = l;
0059     nData  = length(x);
0060     ovfact = 1 / (1 - olap);
0061     davg   = (((nData - segLen)) * ovfact) / segLen + 1;
0062     navg   = round(davg);
0063 
0064     % Compute steps between segments
0065     if navg == 1
0066       shift = 1;
0067     else
0068       shift = (nData - segLen) / (navg - 1);
0069     end
0070     if shift < 1 || isnan(shift)
0071       shift = 1;
0072     end
0073 
0074     %   disp(sprintf('Seglen: %d\t | Shift: %f\t | navs: %d', segLen, shift, navg))
0075 
0076     start = 1.0;
0077     for ii = 1:navg
0078       % compute start index
0079       istart = round (start);
0080       start  = start + shift;
0081 
0082       % get segment
0083       xs  = [x(istart:istart+l-1)].';
0084 
0085       % detrend segment
0086       switch order
0087         case -1
0088           % do nothing
0089         case 0
0090           xs = xs - mean(xs);
0091         case 1
0092           xs = detrend(xs);
0093         otherwise
0094           xs = polydetrend(xs, order);
0095       end
0096 
0097       % make DFT
0098       a   = C*xs;
0099       A   = A + a*conj(a);
0100 
0101     end
0102 
0103     if mod(j,disp_each) == 0 || j == 1 || j == nf
0104       utils.helper.msg(msg.PROC2, 'averaged %d segments', navg);
0105     end
0106     
0107     A2ns    = 2.0*A/navg;
0108     S1      = win.ws;
0109     S12     = S1*S1;
0110     S2      = win.ws2;
0111     ENBW(j) = fs*S2/S12;
0112     Sxx(j)  = A2ns/fs/S2;
0113     S(j)    = A2ns/S12;
0114 
0115   end % for j=1:nf
0116 
0117   varargout{1} = S;
0118   varargout{2} = Sxx;
0119   varargout{3} = ENBW;
0120 
0121 end
0122

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003