Home > classes > @ao > fngen.m

fngen

PURPOSE ^

FNGEN creates an arbitrarily long time-series based on the input PSD.

SYNOPSIS ^

function varargout = fngen(varargin)

DESCRIPTION ^

 FNGEN creates an arbitrarily long time-series based on the input PSD.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: FNGEN creates an arbitrarily long time-series based on the input PSD.

 CALL:        b = fngen(axx, pl)

 PARAMETERS:
              'Nsecs'  - The number of seconds to produce
                         [default: inverse of PSD length]
              'Win'    - The spectral window to use for blending segments
                         [default: Kaiser -150dB]

 M-FILE INFO: Get information about this methods by calling
              >> ao.getInfo('fngen')

              Get information about a specified set-plist by calling:
              >> ao.getInfo('fngen', 'None')

 VERSION:     $Id: fngen.m,v 1.17 2008/09/05 11:05:29 ingo Exp $

 HISTORY:     12-03-07 M Hewitson
                 Creation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % FNGEN creates an arbitrarily long time-series based on the input PSD.
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: FNGEN creates an arbitrarily long time-series based on the input PSD.
0005 %
0006 % CALL:        b = fngen(axx, pl)
0007 %
0008 % PARAMETERS:
0009 %              'Nsecs'  - The number of seconds to produce
0010 %                         [default: inverse of PSD length]
0011 %              'Win'    - The spectral window to use for blending segments
0012 %                         [default: Kaiser -150dB]
0013 %
0014 % M-FILE INFO: Get information about this methods by calling
0015 %              >> ao.getInfo('fngen')
0016 %
0017 %              Get information about a specified set-plist by calling:
0018 %              >> ao.getInfo('fngen', 'None')
0019 %
0020 % VERSION:     $Id: fngen.m,v 1.17 2008/09/05 11:05:29 ingo Exp $
0021 %
0022 % HISTORY:     12-03-07 M Hewitson
0023 %                 Creation
0024 %
0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0026 
0027 function varargout = fngen(varargin)
0028 
0029   % Check if this is a call for parameters
0030   if utils.helper.isinfocall(varargin{:})
0031     varargout{1} = getInfo(varargin{3});
0032     return
0033   end
0034 
0035   import utils.const.*
0036   utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
0037   
0038   % Collect input variable names
0039   in_names = cell(size(varargin));
0040   for ii = 1:nargin,in_names{ii} = inputname(ii);end
0041 
0042   % Collect all AOs and plists
0043   [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0044   pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0045 
0046   if nargout == 0
0047     error('### fngen cannot be used as a modifier. Please give an output variable.');
0048   end
0049   
0050   % combine plists
0051   pl = combine(pl, getDefaultPlist());
0052 
0053   % Extract necessary parameters
0054   Nsecs = find(pl, 'Nsecs');
0055   swin  = find(pl, 'win');
0056 
0057   % Loop over input AOs
0058   bs = [];
0059   for j=1:length(as)
0060     if ~isa(as(j).data, 'fsdata')
0061       warning('!!! %s expects ao/fsdata objects. Skipping AO %s', mfilename, as(j).name);
0062     else
0063       % Properties of the input PSD
0064       N     = 2*(length(as(j).data.y)-1);
0065       fs    = as(j).data.x(end)*2;
0066       % Extract Fourier components
0067       Ak = sqrt(N*as(j).data.getY*fs);
0068       Ak = [Ak; Ak(end-1:-1:2)]; % make two-sided
0069       % redesign input window for this length
0070       switch swin.type
0071         case 'Kaiser'
0072           swin = specwin('Kaiser', N, swin.psll);
0073         otherwise
0074           swin = specwin(swin.type, N);
0075       end
0076       % Compute time-series segments
0077       Olap   = 1-swin.rov/100;
0078       win    = [swin.win].';
0079       segLen = N/fs;
0080       if segLen > Nsecs
0081         cNsecs = 2*segLen;
0082       else
0083         cNsecs = Nsecs;
0084       end
0085       Nsegs  = 1+floor(cNsecs/segLen/Olap);
0086 
0087       % Prepare for generation
0088       rphi = zeros(N,1);                   % Empty vector for random phases
0089       xs   = zeros(fs*(cNsecs+segLen), 1);  % Large empty vector for new time-series
0090       e1   = 1; e2 = segLen*fs;            % Indices into large vector
0091       step = round(segLen*fs*Olap);        % step size between each new segment
0092       lxs  = length(xs);
0093 
0094       % Loop over segments
0095       for s=1:Nsegs
0096         % Generate random phase vector
0097         rphi(2:N/2) = pi*rand(1,N/2-1);  % First half
0098         rphi(N/2+1) = pi*round(rand);    % mid point
0099         rphi(N/2+2:N) = -rphi(N/2:-1:2); % reflected half
0100         %---- Compute Fourier amplitudes
0101         % Use chi^2 distribution to randomize amplitudes.
0102         % - from Percival and Walden: S_est = S.*chi2rnd(2)/2
0103         %   so A_est = A.*sqrt(chi2rnd(2)/2)
0104         % Here we take the measured input data to be a good estimate of
0105         % the underlying power spectrum
0106         X = (Ak.*sqrt(chi2rnd(2)/2)) .*exp(1i.*rphi);
0107         % Inverse FFT
0108         x  = ifft(X, 'symmetric');
0109         % overlap the segments
0110         xs(e1:e2) = xs(e1:e2) + win.*x;
0111         % increase step
0112         e1 = e1 + step;
0113         e2 = e2 + step;
0114         if e2>lxs
0115           break
0116         end
0117       end
0118       % Make ao from the segment of data we want
0119       e1 = fs*segLen/2;
0120       e2 = fs*(Nsecs+segLen/2)-1;
0121       b  = ao(tsdata(xs(e1:e2).', fs));
0122       b.setName(sprintf('fngen(%s)', ao_invars{j}), 'internal');
0123       b.setXunits('s', 'internal');
0124       % Add history
0125       b.addHistory(getInfo, pl, ao_invars(j), as(j).hist);
0126       % Add to outputs
0127       bs = [bs b];
0128     end
0129   end
0130 
0131   % Set outputs
0132   varargout{1} = bs;
0133 end
0134 
0135 %--------------------------------------------------------------------------
0136 % Get Info Object
0137 %--------------------------------------------------------------------------
0138 function ii = getInfo(varargin)
0139   if nargin == 1 && strcmpi(varargin{1}, 'None')
0140     sets = {};
0141     pl   = [];
0142   else
0143     sets = {'Default'};
0144     pl   = getDefaultPlist;
0145   end
0146   % Build info object
0147   ii = minfo(mfilename, 'ao', '', utils.const.categories.sigproc, '$Id: fngen.m,v 1.17 2008/09/05 11:05:29 ingo Exp $', sets, pl);
0148 end
0149 
0150 %--------------------------------------------------------------------------
0151 % Get Default Plist
0152 %--------------------------------------------------------------------------
0153 function pl_default = getDefaultPlist()
0154   pl_default = plist('Nsecs', -1, 'Win', getappdata(0, 'ltpda_default_spectral_window'));
0155 end
0156

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