00001 /* 00002 * THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE 00003 * CHECKED IT OUT USING THE COMMAND CHECKOUT. 00004 * 00005 * $Id: gma_8c-source.html 2161 2006-05-19 16:55:03Z paulf $ 00006 * 00007 * Revision history: 00008 * $Log$ 00008 * Revision 1.1 2006/05/19 16:55:02 paulf 00008 * first inclusion 00008 * 00009 * Revision 1.1 2002/02/28 17:00:39 lucky 00010 * Initial revision 00011 * 00012 * 00013 */ 00014 00015 /* gma.c: Routines for the ground motion analyzer: compute PGA, PGV, PGD * 00016 * and spectral response, using spectral techniques */ 00017 00018 #include <stdio.h> 00019 #include <stdlib.h> 00020 #include <string.h> 00021 #include <math.h> 00022 #include <fft_prep.h> 00023 #include <fft99.h> 00024 #include <transfer.h> 00025 #include <gma.h> 00026 00027 #define MAXLINE 80 00028 00029 static int Debug = 0; 00030 00031 /* Internal Function Prototypes */ 00032 static void makeInt(long nfft, double deltat, double *fIntI); 00033 static int makePsa(double, double, ResponseStruct *); 00034 00035 void cmpmax(long kug, double *ug, double w, double damp, double delta, 00036 double *zd, double *zv, double *za); 00037 00038 00039 /* 00040 * gma: converts a waveform (time series) from its original response 00041 * function to a new response function. This conversion is done 00042 * in the frequency domain. The frequency response of the 00043 * transfer function may be tapered. The input data will be 00044 * padded in the time-domain. The amount of padding is determined 00045 * automatically unless the user provides her own pad length. 00046 * Arguments: input: array of data for preocessing 00047 * npts: number of data points to process 00048 * deltat: time interval between samples, in seconds 00049 * origRS: structure defining process that generated the input data 00050 * that is, the response function to be removed 00051 * fTaper: array of four frequencies (f0, f1, f2, f3) defining the 00052 * taper to be applied to the frequency response function 00053 * before it is convolved with the data. Below f0 and above 00054 * f3, the taper is 0; between f2 and f3 the taper is 1; 00055 * between f0-f1 and f2-f3 is a cosine taper. 00056 * tTaper: length of cosine taper in seconds to be applied at 00057 * both ends of the input data to provide a smooth 00058 * transition into the zero-padded area. 00059 * padlen: The pad length to be applied to data before transforming 00060 * to frequency domain. If padlen < 0, pad length will be 00061 * estimated here and the value chosen will be returned 00062 * in this return-value parameter. 00063 * nfft: The size of the FFT chosen, based on npts + *padlen 00064 * If the returned value of nfft + padlen is less than 00065 * npts, then convertWave had to effectively truncate the 00066 * raw trace in order to fit the processed trace in 00067 * the limit of outBufLen. 00068 * spectPer: array of periods at which the Spectral Response is 00069 * to be calculated 00070 * nsp; number of periods in spectPer 00071 * spectDamp: array of damping values at which the Spectral Response 00072 * is to be calculated 00073 * nsd; number of damping values in spectDamp 00074 * output: array of the six sets of values output from the analysis. 00075 * This array must be allocated by the caller, of length 00076 * outBufLen * (3 + nsd * nsf). 00077 * outBufLen: size of one of the `output' arrays. 00078 * work: a work array that must be allocated by the caller. 00079 * It's size must be outBufLen * 3. 00080 * workFFT: a work array for the FFT routines. It's size must be 00081 * outBufLen * (2 + nsf * nsd). 00082 * 00083 * Returns: 0 on success 00084 * -1 on out-of-memory errors 00085 * -2 on too-small impulse response 00086 * -3 on invalid arguments 00087 */ 00088 00089 int gma(double input[], long npts, double deltat, ResponseStruct *origRS, 00090 double fTaper[4], double tTaper, long *pPadLen, long *pnfft, 00091 double spectPer[], int nsp, double spectDamp[], int nsd, 00092 double output[], long outBufLen, double *work, double *workFFT) 00093 { 00094 static ResponseStruct accRS, psaRS; 00095 FACT *pfact; 00096 double *acc, *vel, *disp, *psa; 00097 double *fAccR, *fAccI, *fIntR, *fIntI; 00098 double *fPsaR, *fPsaI; 00099 double dre, dim; /* temporary values */ 00100 double delfreq, f, rnfft, tpr, accr, acci, minDamp, maxPeriod; 00101 long flen, i, ii, retval = 0; 00102 long nfft, nfreq, trial_nfft, psaPadLen, ntaper; 00103 int nz, np, isd, isp, isdp; 00104 00105 /* Some handy pointers for all our arrays */ 00106 acc = output; 00107 vel = &output[outBufLen]; 00108 disp = &output[outBufLen * 2]; 00109 flen = outBufLen/2; 00110 fAccR = work; 00111 fAccI = &work[flen]; 00112 fIntR = &work[flen * 2]; /* We never use this; it is treated as zero */ 00113 fIntI = &work[flen * 3]; 00114 fPsaR = &work[flen * 4]; 00115 fPsaI = &work[flen * 5]; 00116 00117 /* Validate arguments */ 00118 if (origRS == (ResponseStruct *)NULL || npts < 2 || deltat <= 0.0 || 00119 fTaper == (double *)NULL || outBufLen < npts) 00120 { 00121 return -3; 00122 } 00123 if (fTaper[0] > fTaper[1] || fTaper[1] >= fTaper[2] || 00124 fTaper[2] > fTaper[3] || origRS->dGain == 0.0 ) 00125 { 00126 return -3; 00127 } 00128 00129 /* Adjust the input response function to get ground acceleration * 00130 * Input units: nanometers of displacement * 00131 * Output units: acceleration in cm/sec/sec */ 00132 00133 accRS.dGain = 1.0 / (1.0e+7 * origRS->dGain); 00134 00135 accRS.iNumPoles = origRS->iNumZeros; 00136 accRS.iNumZeros = origRS->iNumPoles + 2; /* Differentiate to acceleration */ 00137 00138 if ( (accRS.Poles = (PZNum *)malloc(sizeof(PZNum) * accRS.iNumPoles)) == 00139 (PZNum *)NULL || 00140 (accRS.Zeros = (PZNum *)malloc(sizeof(PZNum) *accRS.iNumZeros)) == 00141 (PZNum *)NULL) 00142 { 00143 retval = -1; 00144 goto exit; 00145 } 00146 00147 00148 /* Copy the poles and zeros, using structure copy */ 00149 np = nz = 0; 00150 00151 for (i = 0; i < origRS->iNumPoles; i++) 00152 accRS.Zeros[nz++] = origRS->Poles[i]; 00153 00154 for (i = 0; i < origRS->iNumZeros; i++) 00155 accRS.Poles[np++] = origRS->Zeros[i]; 00156 00157 for (i = 0; i < 2; i++) 00158 { 00159 accRS.Zeros[nz].dReal = 0.0; 00160 accRS.Zeros[nz].dImag = 0.0; 00161 nz++; 00162 } 00163 00164 /* Remove any cancelling poles and zeros */ 00165 pzCancel(&accRS, 1.0e-5); 00166 00167 /* Determine how much padding we need, unless the caller told us. */ 00168 if ( *pPadLen < 0) 00169 { 00170 if ( (*pPadLen = respLen( &accRS, deltat, fTaper)) < 0) 00171 { 00172 /* Some error occured */ 00173 retval = *pPadLen; 00174 if (Debug) 00175 logit ("", "\nrespLen error: %ld\n", *pPadLen); 00176 goto exit; 00177 } 00178 00179 /* 00180 * Find the length in samples for the lightly damped PSA oscillator 00181 * to decay to 1% of it's maximum: a measure of the length of its 00182 * convolution filter. Use the longest period for this. 00183 */ 00184 /* Find the minumum damping and maximum period, to get the longest decay */ 00185 minDamp = spectDamp[0]; 00186 for (isd = 1; isd < nsd; isd++) 00187 if (minDamp > spectDamp[isd]) minDamp = spectDamp[isd]; 00188 maxPeriod = spectPer[0]; 00189 for (isp = 1; isp < nsp; isp++) 00190 if (maxPeriod < spectPer[isp]) maxPeriod = spectPer[isp]; 00191 00192 psaPadLen = (long)( -1.0 * maxPeriod * log(0.01) / 00193 (minDamp * 2.0 * PI * deltat)); 00194 if (*pPadLen < psaPadLen) 00195 *pPadLen = psaPadLen; 00196 if (Debug) 00197 logit ("", "estimated pad length: %ld\n", *pPadLen); 00198 } 00199 00200 /* Find a convenient FFT size for our data plus padding that will fit in 00201 * outBuf */ 00202 trial_nfft = *pPadLen + npts; 00203 while ( (nfft = prepFFT(trial_nfft, &pfact)) > outBufLen - FFT_EXTRA) 00204 { 00205 if (nfft < 0) 00206 { 00207 /* Out of memory */ 00208 retval = nfft; 00209 goto exit; 00210 } 00211 trial_nfft -= 100; /* Try a liitle bit smaller */ 00212 } 00213 00214 if (nfft - *pPadLen < npts) 00215 npts = nfft - *pPadLen; /* Chop some off the end if it won't fit */ 00216 /* We aren't passing this new value of npts back to our caller; she'll 00217 * have to figure it out from the values of pPadLen and nfft passed back */ 00218 00219 nfreq = nfft / 2; 00220 00221 /* Copy the input data into the acc array in preparation for FFT */ 00222 for (i = 0; i < npts; i++) 00223 acc[i] = input[i]; 00224 for (i = npts; i < nfft + FFT_EXTRA; i++) /* Do the null-padding */ 00225 acc[i] = 0.0; 00226 00227 00228 /* Taper the data into the zero-padded area */ 00229 ntaper = (long)(tTaper / deltat); 00230 taper(acc, npts, ntaper); 00231 00232 /* Transform the padded data into frequency domain */ 00233 fft991(acc, workFFT, pfact->trigs, pfact->ifax, 1, nfft, nfft, 1, -1); 00234 00235 /* Fill in all the transfer functions */ 00236 response(nfft, deltat, &accRS, fAccR, fAccI); 00237 makeInt(nfft, deltat, fIntI); 00238 00239 /* Compute acceleration, velocity and displacement output data in the 00240 frequency domain */ 00241 rnfft = 1.0 / (2*nfft); 00242 00243 acc[0] = 0.0; /* Remove the mean, if there is any */ 00244 vel[0] = 0.0; 00245 disp[0] = 0.0; 00246 acc[1] = 0.0; 00247 vel[1] = 0.0; 00248 disp[1] = 0.0; 00249 00250 delfreq = 1.0 / (nfft * deltat); 00251 for (i = 1; i < nfreq; i++) 00252 { 00253 ii = i+i; 00254 f = i * delfreq; 00255 tpr = ftaper(f, fTaper[1], fTaper[0]) * ftaper(f, fTaper[2], fTaper[3]); 00256 dre = acc[ii]; /* Real part of transformed data */ 00257 dim = acc[ii+1]; /* Imaginary part of transformed data */ 00258 /* The acceleration */ 00259 /* WAS: scaled by 1 * rnfft */ 00260 accr = (dre * fAccR[i] - dim * fAccI[i]) * tpr; 00261 acci = (dre * fAccI[i] + dim * fAccR[i]) * tpr; 00262 acc[ii] = accr; 00263 acc[ii+1] = acci; 00264 /* Integrate once to get velocity */ 00265 vel[ii] = fIntI[i] * acci; 00266 vel[ii+1] = -fIntI[i] * accr; 00267 /* Integrate again to get displacement */ 00268 disp[ii] = -fIntI[i] * vel[ii+1]; 00269 disp[ii+1] = fIntI[i] * vel[ii]; 00270 } 00271 /* Remove the Nyquist frequency component */ 00272 acc[nfft] = 0.0; 00273 vel[nfft] = 0.0; 00274 disp[nfft] = 0.0; 00275 00276 /* Compute the Spectral Responses for all their periods and dampings */ 00277 isdp = -1; 00278 for (isp = 0; isp < nsp; isp++) 00279 { 00280 for (isd = 0; isd < nsd; isd++) 00281 { 00282 isdp++; 00283 psa = &output[outBufLen * (3 + isdp)]; 00284 if ( (retval = makePsa(spectPer[isp], spectDamp[isd], &psaRS)) < 0) 00285 goto exit; 00286 response(nfft, deltat, &psaRS, fPsaR, fPsaI); 00287 00288 psa[0] = 0.0; 00289 psa[1] = 0.0; 00290 00291 for (i = 1; i < nfreq; i++) 00292 { 00293 ii = i+i; 00294 accr = acc[ii]; 00295 acci = acc[ii+1]; 00296 /* The spectral response accelerations */ 00297 psa[ii] = accr * fPsaR[i] - acci * fPsaI[i] + accr; 00298 psa[ii+1] = accr * fPsaI[i] + acci * fPsaR[i] + acci; 00299 } 00300 /* Remove the Nyquist frequency component */ 00301 psa[nfft] = 0.0; 00302 psa[nfft+1] = 0.0; 00303 00304 cleanPZ( &psaRS ); 00305 } 00306 } 00307 00308 /* Transform the whole output array back to time domain */ 00309 fft991(output, workFFT, pfact->trigs, pfact->ifax, 1, outBufLen, nfft, 00310 3 + isd * isp, +1); 00311 00312 /* We're done! Give the data back to the caller */ 00313 *pnfft = nfft; 00314 00315 exit: 00316 /* Clean up before we go home */ 00317 cleanPZ( &accRS ); 00318 cleanPZ( &psaRS ); 00319 00320 return retval; 00321 } 00322 00323 /* 00324 * makeInt: make a `response' function for integrating once with respect 00325 * to time. This function is simply `i' times the frequency; 00326 * or the response function for a single zero at the origin. 00327 */ 00328 static void makeInt(long nfft, double deltat, double *fIntI) 00329 { 00330 double delomg, omega; 00331 long i, ntr; 00332 00333 ntr = nfft / 2; 00334 delomg = 2.0 * PI / (nfft * deltat); 00335 00336 fIntI[0] = 0.0; 00337 for (i = 1; i < ntr; i++) 00338 { 00339 omega = delomg * i; 00340 fIntI[i] = 1.0 / omega; 00341 } 00342 return; 00343 } 00344 00345 /* 00346 * makePsa: Generate the acceleration response function in pole-zero-gain 00347 * form for a simple oscillator with given free period and damping. 00348 */ 00349 static int makePsa(double period, double damp, ResponseStruct *rs) 00350 { 00351 double mu, omega; 00352 00353 omega = 2.0 * PI / period; 00354 mu = sqrt(1.0 - damp * damp); 00355 00356 rs->iNumZeros = 2; 00357 rs->iNumPoles = 2; 00358 if ( (rs->Zeros = (PZNum *)malloc(2 * sizeof(PZNum))) == (PZNum *)0 || 00359 (rs->Poles = (PZNum *)malloc(2 * sizeof(PZNum))) == (PZNum *)0 ) 00360 return -1; 00361 00362 rs->Zeros[0].dReal = rs->Zeros[0].dImag = 0.0; 00363 rs->Zeros[1].dReal = rs->Zeros[1].dImag = 0.0; 00364 00365 rs->dGain = -1.0; 00366 00367 rs->Poles[0].dReal = - omega * damp; 00368 rs->Poles[0].dImag = omega * mu; 00369 rs->Poles[1].dReal = - omega * damp; 00370 rs->Poles[1].dImag = - omega * mu; 00371 00372 return 0; 00373 } 00374 00375 void gmaDebug( int level ) 00376 { 00377 Debug = level; 00378 return; 00379 }