


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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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