Home > m > sigproc > frequency_domain > ltpda_lpsd.m

ltpda_lpsd

PURPOSE ^

LTPDA_LPSD implement LPSD algorithm for analysis objects.

SYNOPSIS ^

function bs = ltpda_lpsd(varargin)

DESCRIPTION ^

 LTPDA_LPSD implement LPSD algorithm for analysis objects.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 LTPDA_LPSD implement LPSD algorithm for analysis objects.
 
   >> bs = ltpda_lpsd(as)
 or
   >> bs = ltpda_lpsd(as, pl)
 
   Inputs:
     as  - array of analysis objects
     pl  - parameter list (see below)
 
   Outputs:
     bs  - array of analysis objects, one for each input
 
   Parameter list:
     Kdes  - desired number of averages   (default 100)
     Kmin  - minimum number of averages   (default 1)
     Jdes  - number of spectral frequencies to compute (default fs/2)
     win   - a specwin window object
             Only the design parameters of the window object are used; the
             window is recomputed for each DFT length inside the lpsd_core
             algorithm.
     Nolap - desired overlap (default 0.77; later taken from window)
     Order - order of detrending
 
 M Hewitson 02-02-07
 
 $Id: ltpda_lpsd.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function bs = ltpda_lpsd(varargin)
0002 % LTPDA_LPSD implement LPSD algorithm for analysis objects.
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % LTPDA_LPSD implement LPSD algorithm for analysis objects.
0006 %
0007 %   >> bs = ltpda_lpsd(as)
0008 % or
0009 %   >> bs = ltpda_lpsd(as, pl)
0010 %
0011 %   Inputs:
0012 %     as  - array of analysis objects
0013 %     pl  - parameter list (see below)
0014 %
0015 %   Outputs:
0016 %     bs  - array of analysis objects, one for each input
0017 %
0018 %   Parameter list:
0019 %     Kdes  - desired number of averages   (default 100)
0020 %     Kmin  - minimum number of averages   (default 1)
0021 %     Jdes  - number of spectral frequencies to compute (default fs/2)
0022 %     win   - a specwin window object
0023 %             Only the design parameters of the window object are used; the
0024 %             window is recomputed for each DFT length inside the lpsd_core
0025 %             algorithm.
0026 %     Nolap - desired overlap (default 0.77; later taken from window)
0027 %     Order - order of detrending
0028 %
0029 % M Hewitson 02-02-07
0030 %
0031 % $Id: ltpda_lpsd.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $
0032 %
0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0034 
0035 ALGONAME = mfilename;
0036 VERSION  = '$Id: ltpda_lpsd.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $';
0037 
0038 % capture input variable names
0039 invars = {};
0040 as     = [];
0041 pl     = [];
0042 for j=1:nargin
0043   invars = [invars cellstr(inputname(j))];
0044   if isa(varargin{j}, 'ao')
0045     as = [as varargin{j}];
0046   end
0047   if isa(varargin{j}, 'plist')
0048     pl = varargin{j};
0049   end
0050 end
0051 
0052 % initialise output array
0053 bs = []; 
0054 
0055 % check plist
0056 if isempty(pl)
0057   pl = plist();
0058 end
0059 
0060 % Go through each input AO
0061 na = length(as);
0062 for i=1:na
0063   a = as(i);
0064   d = a.data;
0065   
0066   % check this is a time-series object
0067   if ~isa(d, 'tsdata')
0068     error('### lpsd requires tsdata (time-series) inputs.');
0069   end
0070   
0071   %--- check input parameters
0072 
0073   % Averaging
0074   Kdes = find(pl, 'Kdes');   % desired averages
0075   if isempty(Kdes)
0076     Kdes = 100;
0077     disp(sprintf('! Using default Kdes of %d', Kdes))
0078     pl = append(pl, param('Kdes', Kdes));
0079   end
0080   Kmin = find(pl, 'Kmin');    % minimum averages
0081   if isempty(Kmin)
0082     Kmin = 1;
0083     disp(sprintf('! Using default Kmin of %d', Kmin))
0084     pl = append(pl, param('Kmin', Kmin));
0085   end
0086 
0087   Jdes = find(pl, 'Jdes');    % num desired spectral frequencies
0088   if isempty(Jdes)
0089     Jdes = round(d.fs/4);
0090     disp(sprintf('! Using default Jdes of %d', Jdes))
0091     pl = append(pl, param('Jdes', Jdes));
0092   end
0093 
0094   Win = find(pl, 'Win');
0095   if isempty(Win)
0096     Win = specwin('Kaiser', 100, 250);
0097     pl = append(pl, param('Win', Win));
0098   end
0099   
0100   Nolap = find(pl, 'Nolap');    % desired overlap
0101   if isempty(Nolap)
0102     % use window function rov
0103     Nolap = Win.rov/100;
0104     disp(sprintf('! Using recommended Nolap of %d', Nolap))
0105     pl = append(pl, param('Nolap', Nolap));
0106   end
0107 
0108   Order = find(pl, 'Order');    % desired overlap
0109   if isempty(Order)
0110     Order = 0;
0111     disp(sprintf('! Using default detrending order of %d', Order))
0112     pl = append(pl, param('Order', Order));
0113   end
0114   
0115   % Get frequency vector
0116   [f,r,m,L,rr,rrr] = ltpda_compute_f(d.fs, length(d.x), Kdes, Kmin, Jdes, Nolap);
0117   
0118   % compute LPSD
0119   [Q, Qxx, ENBW] = mlpsd(d.x, f, r, m, L, d.fs, Win, Order);
0120 
0121   % create new output fsdata
0122   fs = fsdata(f, sqrt(Qxx), d.fs);
0123   fs = set(fs, 'name', sprintf('lpsd of %s', d.name));
0124   fs = set(fs, 'yunits', [d.yunits '/\surdHz']);
0125   fs = set(fs, 'enbw', ENBW);
0126   
0127   % create new output history
0128   h = history(ALGONAME, VERSION, pl, a.hist);
0129   h = set(h, 'invars', invars);
0130 
0131   % make output analysis object
0132   b = ao(fs, h);
0133   
0134   % set name
0135   if isempty(invars{1})
0136     n = a.name;
0137   else
0138     n = invars{1};
0139   end
0140   b = set(b, 'name', sprintf('lpsd(%s)', n));
0141   
0142   bs = [bs b];
0143 end
0144 
0145 %--------------------------------------------------------------------------
0146 %
0147 function varargout = mlpsd(varargin)
0148 
0149 % MLPSD
0150 %
0151 % [S,Sxx,ENBW] = mlpsd(x,f,r,m,L,fs,win,order)
0152 %
0153 % M Hewitson 19-05-07
0154 %
0155 %
0156 
0157 % Get inputs
0158 x     = varargin{1};
0159 f     = varargin{2};
0160 r     = varargin{3};
0161 m     = varargin{4};
0162 L     = varargin{5};
0163 fs    = varargin{6};
0164 win   = varargin{7};
0165 order = varargin{8};
0166 
0167 twopi    = 2.0*pi;
0168 
0169 nx   = length(x);
0170 nf   = length(f);
0171 Sxx  = zeros(nf,1);
0172 S    = zeros(nf,1);
0173 ENBW = zeros(nf,1);
0174 
0175 disp_each = round(nf/100)*10;
0176 
0177 for j=1:nf
0178   
0179   if mod(j,disp_each) == 0 || j == 1 || j == nf
0180     disp(sprintf('++ computing frequency %d of %d: %f Hz', j, nf, f(j)));
0181   end
0182 
0183   % compute DFT exponent and window
0184   l = L(j);
0185   switch win.name
0186     case 'Kaiser'
0187       win = specwin(win.name, l, win.psll);
0188     otherwise
0189       win = specwin(win.name, l);
0190   end
0191 
0192   p     = 1i * twopi * m(j)/l.*[0:l-1];
0193   C     = win.win .* exp(p);
0194   Xolap = (1-win.rov/100);
0195   % do segments
0196   A  = 0.0;
0197   start = 1;
0198   nsegs = 0;
0199   while start + l - 1 <= nx
0200 
0201     % get segment
0202     xs  = x(start:start+l-1);
0203     
0204     % detrend segment
0205     switch order
0206       case -1
0207         % do nothing
0208       case 0
0209         xs = xs - mean(xs);
0210       case 1
0211         xs = detrend(xs);
0212       otherwise
0213         xs = polydetrend(xs, order);
0214     end
0215 
0216     % make DFT
0217     a   = C*xs;
0218     A   = A + a*conj(a);
0219     
0220     % make next step
0221     start = floor(start + l*Xolap);
0222     nsegs = nsegs + 1;
0223   end
0224 
0225   if mod(j,disp_each) == 0 || j == 1 || j == nf
0226     disp(sprintf('   averaged %d segments', nsegs));
0227   end
0228 
0229   A2ns    = 2.0*A/nsegs;
0230   S1      = win.ws;
0231   S12     = S1*S1;
0232   S2      = win.ws2;
0233   ENBW(j) = fs*S2/S12;
0234   Sxx(j)  = A2ns/fs/S2;
0235   S(j)    = A2ns/S12;
0236   
0237 end
0238 
0239 varargout{1} = S;
0240 varargout{2} = Sxx;
0241 varargout{3} = ENBW;
0242 
0243 
0244 
0245 
0246 % END

Generated on Fri 08-Jun-2007 16:09:11 by m2html © 2003