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

Generated on Tue May 6 09:16:03 2003 for Earthworm Libs by doxygen1.3-rc3