Main Page Class Hierarchy Compound List File List Compound Members File Members

# gma.c

Go to the documentation of this file.

    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,v 1.1 2006/05/19 16:55:02 paulf Exp $
    00006  *
    00007  *    Revision history:
    00008  *     $Log: gma_8c-source.html,v $
    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 }