MAKEH - generates metric perturbation for GW source
hij = makeh(simID, ampl, startSamp, sampFreq);
simId string containing signal information
ampl fixed signal amplitude
startSamp start time of (in samples) of signal; mostly interested in
fractional part to properly line up signal samples with
data samples
sampFreq sampling frequency
hij metric perturbation
Currently supported signals are:
Gaussians 'G'
sine-Gaussians 'SG'
cosine-Gaussian 'CG'
BH ringdowns 'BH'
sine 'SIN'
cosine 'COS'
sinc 'SINC'
linear chirp 'CHIRP'
band limited white noise 'BWHITE'
Gaussian modulated white noise 'GWHITE'
Gaussian modulated linear chirp 'GCHIRP'
square wave 'SQUARE'
Gaussian modulated square wave 'SQUAREG'
sawtooth wave 'SAWTOOTH'
Gaussian modulated sawtooth wave 'SAWTOOTHG'
chirplet* 'CHIRPLET'
* 'CHIRPLET' is generated in graven.m, not makeh.m
ERRORS:
badSignalType - unrecognized simID
sampFreq - sampling frequency doesn't match IFO samp freq
durationSec - signal duration in seconds poorly formed
durationSamp - signal duration less than one sample
tiEmpty - sampling times vector is empty
badTau - tau value poorly formed
unknownType - unrecognized waveform type
$Id: makeh.m,v 1.33 2005/05/12 15:42:14 stuver Exp $
SEE ALSO: GRAVEN, SIMPARSE0001 function hij = makeh(simId, ampl, startSamp, sampFreq) 0002 % MAKEH - generates metric perturbation for GW source 0003 % 0004 % hij = makeh(simID, ampl, startSamp, sampFreq); 0005 % 0006 % simId string containing signal information 0007 % ampl fixed signal amplitude 0008 % startSamp start time of (in samples) of signal; mostly interested in 0009 % fractional part to properly line up signal samples with 0010 % data samples 0011 % sampFreq sampling frequency 0012 % 0013 % hij metric perturbation 0014 % 0015 % Currently supported signals are: 0016 % Gaussians 'G' 0017 % sine-Gaussians 'SG' 0018 % cosine-Gaussian 'CG' 0019 % BH ringdowns 'BH' 0020 % sine 'SIN' 0021 % cosine 'COS' 0022 % sinc 'SINC' 0023 % linear chirp 'CHIRP' 0024 % band limited white noise 'BWHITE' 0025 % Gaussian modulated white noise 'GWHITE' 0026 % Gaussian modulated linear chirp 'GCHIRP' 0027 % square wave 'SQUARE' 0028 % Gaussian modulated square wave 'SQUAREG' 0029 % sawtooth wave 'SAWTOOTH' 0030 % Gaussian modulated sawtooth wave 'SAWTOOTHG' 0031 % chirplet* 'CHIRPLET' 0032 % * 'CHIRPLET' is generated in graven.m, not makeh.m 0033 % 0034 % ERRORS: 0035 % badSignalType - unrecognized simID 0036 % sampFreq - sampling frequency doesn't match IFO samp freq 0037 % durationSec - signal duration in seconds poorly formed 0038 % durationSamp - signal duration less than one sample 0039 % tiEmpty - sampling times vector is empty 0040 % badTau - tau value poorly formed 0041 % unknownType - unrecognized waveform type 0042 % 0043 % $Id: makeh.m,v 1.33 2005/05/12 15:42:14 stuver Exp $ 0044 % 0045 % SEE ALSO: GRAVEN, SIMPARSE 0046 0047 hij = []; 0048 0049 % GET signal information from simId 0050 % CHECK for errors 0051 % GENERATE time vector 0052 0053 % Make simId all upper case for ease in using strcmp 0054 simId = upper(simId); 0055 % Get signal information from simId 0056 [id, f0, tau, durationSec, alpha] = simparse(simId); 0057 0058 % Check that simparse did not return an empty value 0059 if (isempty(id) || isempty(f0) || isempty(tau) || isempty(durationSec)) 0060 msgId = 'makeh:simparse'; 0061 fmtstr = [fmtstr, '%s\n\tnumel(id) = %d\n']; 0062 fmtstr = [fmtstr, '\tnumel(f0) = %d\n']; 0063 fmtstr = [fmtstr, '\tnumel(tau) = %d\n']; 0064 fmtstr = [fmtstr, '\tnumel(durationSec) = %d']; 0065 error(msgId,fmtstr,msgId,numel(id),numel(f0),... 0066 numel(tau),numel(durationSec)); 0067 end 0068 0069 % Check that sampFreq is > 0 0070 if (sampFreq <= 0 || ... 0071 ~isreal(sampFreq) || ... 0072 isnan(sampFreq) || ... 0073 isinf(sampFreq)) 0074 msgId = 'makeh:sampFreq'; 0075 error(msgId,'%s: %d',msgId,sampFreq); 0076 end 0077 0078 % Check that durationSec is > 0 0079 if (durationSec <= 0 || ... 0080 ~isreal(durationSec) || ... 0081 isnan(durationSec) || ... 0082 isinf(durationSec)) 0083 msgId = 'makeh:durationSec'; 0084 error(msgId,'%s: %d',msgId,durationSec'); 0085 end 0086 0087 % Deterimine length of time vector in samples (floor trucated) 0088 durationSamp = floor(durationSec*sampFreq); 0089 % Check that durationSamp is not < 1 0090 if (durationSamp < 1) 0091 msgId = 'makeh:durationSamp'; 0092 error(msgId,'%s: %d',msgId,durationSamp); 0093 end 0094 % Time vector = times at which IFO samples waveform 0095 ti = ((1:durationSamp)./sampFreq)+((1-mod(startSamp,1))/sampFreq); 0096 % N.B. This accounts for intersample start times 0097 0098 % Check to make sure ti is not empty 0099 if isempty(ti) 0100 msgId = 'makeh:tiEmpty'; 0101 error(msgId,msgId); 0102 end 0103 0104 % GENERATE needed base functions 0105 0106 % The logical(sum(strcmp(id, {...}))) makes it so that if any one the 0107 % strings match id, then the appropriate base is generated 0108 if sum(strcmp(id, {'G', 'SG', 'CG', 'SQUAREG', 'SAWTOOTHG', 'GWHITE', 'GCHIRP'})) > 0 0109 % A gaussian base needs to be generated 0110 0111 % Check that tau is > 0 0112 if (~(isreal(tau) && tau > 0 && ~isinf(tau))) 0113 msgId = 'makeh:badTau'; 0114 error(msgId,'%s: %d',msgId,tau); 0115 end 0116 0117 % Generate the Gaussian waveform that is also the base for SG and CG 0118 % Offset times for Gaussians 0119 tg = ti-(ti(end)/2); 0120 % Generate Gaussian 0121 gauss = exp(-(tg/tau).^2); 0122 end 0123 0124 if sum(strcmp(id, {'SIN', 'SG'})) > 0 0125 % A sine wave base needs to be generated 0126 0127 % Check that f >= 0 for periodic function 0128 checkperiodic(f0); 0129 0130 sine = ampl.*sin(2*pi*f0*ti); 0131 end 0132 0133 if sum(strcmp(id, {'COS', 'CG'})) > 0 0134 % Check that f >= 0 for periodic function 0135 checkperiodic(f0); 0136 0137 % A cosine wave base needs to be generated 0138 cosine = ampl.*cos(2*pi*f0*ti); 0139 end 0140 0141 if sum(strcmp(id, {'SQUARE', 'SQUAREG'})) > 0 0142 % Check that f >= 0 for periodic function 0143 checkperiodic(f0); 0144 0145 % A square wave base needs to be generated 0146 squ = ampl*square(2*pi*f0*ti); 0147 end 0148 0149 if sum(strcmp(id, {'SAWTOOTH', 'SAWTOOTHG'})) > 0 0150 % Check that f >= 0 for periodic function 0151 checkperiodic(f0); 0152 0153 % A sawtooth wave base needs to be generated 0154 saw = ampl*sawtooth(2*pi*f0*ti); 0155 end 0156 0157 % GENERATE desired metric perturbation based on simId 0158 % RETURN appropriate timeseries 0159 0160 switch id 0161 case 'G' 0162 hij = make33time(gauss); 0163 case 'SIN' 0164 hij = make33time(sine); 0165 case 'SG' 0166 hij = make33time(sine.*gauss); 0167 case 'COS' 0168 hij = make33time(cosine); 0169 case 'CG' 0170 hij = make33time(cosine.*gauss); 0171 case 'SINC' 0172 func = ampl*sinc(2*f0*ti); 0173 hij = make33time(func); 0174 case 'CHIRP' 0175 % Generate a chirp signal that goes from 0 to f0 frequency 0176 func = ampl*chirp(ti, 0, ti(end), f0); 0177 hij = make33time(func); 0178 case 'GCHIRP' 0179 func = ampl.*chirp(ti, 0, ti(end), f0).*gauss; 0180 hij = make33time(func); 0181 case 'SAWTOOTH' 0182 hij = make33time(saw); 0183 case 'SAWTOOTHG' 0184 hij = make33time(saw.*gauss); 0185 case 'SQUARE' 0186 hij = make33time(squ); 0187 case 'SQUAREG' 0188 hij = make33time(squ.*gauss); 0189 case 'BWHITE' 0190 if isempty(alpha) 0191 warning(['For BWHITE, alpha is used as the bandwidth, '... 0192 'therefore it cannot be empty.']) 0193 end 0194 func = makebwhite(f0, alpha, tau, sampFreq, durationSec); 0195 hij = make33time(func); 0196 case 'BH' 0197 % -- Sam did this (2 lines) without telling me and this is incorrect. 0198 %func = ampl*exp(-ti/tau).*sin(2*pi*f0*ti); 0199 %hij = make33time(func); 0200 0201 % Construct damped, periodic source for BH ringdown 0202 % hij = ampl*exp(-t/tau)*[sin, cos, 0; cos, sin, 0; 0, 0, 0] 0203 0204 % GENERATE periodic part of metric perturbation in 3x3xtime 0205 % array format 0206 % MAKE BH ringdown 0207 0208 % Initialize metric perturbation 0209 hij = zeros(3, 3, length(ti)); 0210 0211 % Generate vector of zeros in time dimension for use in 0212 % constructing the metric perturbation 0213 oh = zeros(1, 1, length(ti)); 0214 % Generate vector of the sine function in the time dimension for 0215 % use in constructing the metric perturbation 0216 xx = reshape(sin(2*pi*f0*ti), [1, 1, length(ti)]); 0217 % Generate vector of the cosine function in the time dimension 0218 % for use in constructing the metric perturbation 0219 xy = reshape(cos(2*pi*f0*ti), [1, 1, length(ti)]); 0220 0221 % Assemble the part of the metric perturbation that deals with 0222 % periodic functions 0223 periodic = [xx, xy, oh; xy, xx, oh; oh, oh, oh]; 0224 % This is in full 3x3xtime array format 0225 0226 % Multiply each time step of the periodic portion of the metric 0227 % of the coresponding time step of the exponential and the 0228 % amplitude 0229 for k = 1:length(ti) 0230 hij(:,:,k) = (ampl*exp(-ti(k)/tau)).*periodic(:,:,k); 0231 end 0232 case 'GWHITE' 0233 func = ampl.*gauss.*randn(1,length(gauss)); 0234 hij= make33time(func); 0235 otherwise 0236 msgId = 'makeh:unknownSignalType'; 0237 error(msgId,'%s: %s',msgId,id); 0238 end 0239 0240 return 0241 0242 % ~~~ END MAIN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0243 0244 function checkperiodic(f0) 0245 0246 % CHECKPERIODIC - check the frequency of a periodic function for errors 0247 % 0248 % checkperiodic(f0) 0249 % 0250 % f0 if f0 is < 0, imaginary, NaN or Inf, the program errors out 0251 0252 % Check that F0 is > 0, real, and finite 0253 if (~(isreal(f0) && f0 > 0 && ~isinf(f0))) 0254 msgId = 'makeh:checkPeriodic'; 0255 error(msgId,'%s: %d',msgId,f0); 0256 end 0257 0258 return 0259 0260 % ~~~ END CHECKPERIODIC ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0261 0262 function hij = make33time(func) 0263 % MAKE33TIME - put a timeseries vector into a 3x3xtime array 0264 % 0265 % hij = make33time(func) 0266 % 0267 % func function timeseries to be copied into the 3x3xtime array 0268 % 0269 % hij 3x3xtime array 0270 0271 % Using the notation of Thorne 1980 RMP: 0272 % Assume h_{jk} = h(t)\mathcal{Y}^{l=2,m=0}_{jk} 0273 % Recall \mathcal{Y}^{lm}_{A_l}N_{A_l} = Y^{lm} 0274 % h(t) for us is a damped sinusoid 0275 % 0276 % Recall Y^{l=2,m=0} = \sqrt{\frac{5}{16\pi}}(3\cos^2\theta-1) 0277 % = \sqrt{\frac{5}{16\pi}}(3n_z^2-1) 0278 % 0279 % Hence 0280 % 0281 % \mathcal{Y}_{jk}^{l=2,m=0} \propto [-1,0,0;0,-1,0,0;0,0,2] 0282 0283 % Choose normalization so that peak value of h+ (over all orientations) is 0284 % peak value of func 0285 func = func/1.5; 0286 0287 hij = zeros([3,3,length(func)]); 0288 hij(1,1,:) = -func; 0289 hij(1,1,:) = -func; 0290 hij(3,3,:) = 2*func; 0291 0292 return