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

transfer.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:
00006  *
00007  *    Revision history:
00008  *     $Log:
00009  *
00010  *
00011  *
00012  */
00013 
00014 /* transfer.c: Routines for dealing with instrument transfer functions */
00015 
00016 #include <stdio.h>
00017 #include <stdlib.h>
00018 #include <string.h>
00019 #include <math.h>
00020 #include <fft_prep.h>
00021 #include <fft99.h>
00022 #include <transfer.h>
00023 
00024 #define MAXLINE 80
00025 static int Debug = 0;
00026 
00027 /* Internal Function Prototypes */
00028 static void drop(PZNum *, int *, int);
00029 
00030 /*
00031  * response: compute frequency response from the pole-zero-gain information.
00032  *  arguments:  nfft: the number of points that will be used in the FFT
00033  *            deltat: the time interval between data points in the time-domain
00034  *               pRS: pointer to the Response Structure holding the poles,
00035  *                    zeros and gain information for the desired function
00036  *               tfr: pointer to the real part of the frequency response
00037  *               tfi: pointer to the imaginary part of the frequency 
00038  *                    response. Both tfr and tfi must be allocated
00039  *                    by the caller to contain at least nfft/2+1 values.
00040  */
00041 void response(long nfft, double deltat, ResponseStruct *pRS, 
00042               double *tfr, double *tfi)
00043 {
00044   double delomg, omega, mag2;
00045   double sr, si, srn, sin, srd, sid, sr0, si0;
00046   long i, j, ntr;
00047   
00048   ntr = nfft / 2 + 1;
00049   delomg = 2.0 * PI / (nfft * deltat);
00050   
00051   /* The (almost) zero frequency term */
00052   /* The zeros, in the numerator */
00053   srn = 1.0;
00054   sin = 0.0;
00055   omega = delomg * 0.001;
00056   for (j = 0; j < pRS->iNumZeros; j++)
00057   {
00058     sr = - pRS->Zeros[j].dReal;
00059     si = omega - pRS->Zeros[j].dImag;
00060     sr0 = srn * sr - sin * si;
00061     si0 = srn * si + sin * sr;
00062     srn = sr0;
00063     sin = si0;
00064   }
00065   
00066   /* The poles; in the denominator */
00067   srd = 1.0;
00068   sid = 0.0;
00069   
00070   for (j = 0; j < pRS->iNumPoles; j++)
00071   {
00072     sr = - pRS->Poles[j].dReal;
00073     si = omega - pRS->Poles[j].dImag;
00074     sr0 = srd * sr - sid * si;
00075     si0 = srd * si + sid * sr;
00076     srd = sr0;
00077     sid = si0;
00078   }
00079   
00080   /* Combine numerator, denominator and gain using complex arithemetic */
00081   mag2 = pRS->dGain / (srd * srd + sid * sid);
00082   tfr[0] = mag2 * (srn * srd + sin * sid);
00083   tfi[0] = 0.0; /* Actually the Nyqust part; we don't want it */
00084   
00085   /* The non-zero frequency parts */
00086   for (i = 1; i < ntr; i++)
00087   {
00088     /* The zeros, in the numerator */
00089     srn = 1.0;
00090     sin = 0.0;
00091     omega = delomg * i;
00092     for (j = 0; j < pRS->iNumZeros; j++)
00093     {
00094       sr = - pRS->Zeros[j].dReal;
00095       si = omega - pRS->Zeros[j].dImag;
00096       sr0 = srn * sr - sin * si;
00097       si0 = srn * si + sin * sr;
00098       srn = sr0;
00099       sin = si0;
00100     }
00101     
00102     /* The poles; in the denominator */
00103     srd = 1.0;
00104     sid = 0.0;
00105     
00106     for (j = 0; j < pRS->iNumPoles; j++)
00107     {
00108       sr = - pRS->Poles[j].dReal;
00109       si = omega - pRS->Poles[j].dImag;
00110       sr0 = srd * sr - sid * si;
00111       si0 = srd * si + sid * sr;
00112       srd = sr0;
00113       sid = si0;
00114     }
00115     
00116     /* Combine numerator, denominator and gain using complex arithemetic */
00117     mag2 = pRS->dGain / (srd * srd + sid * sid);
00118     tfr[i] = mag2 * (srn * srd + sin * sid);
00119     tfi[i] = mag2 * (sin * srd - srn * sid);
00120   }
00121   return;
00122 }
00123 
00124 
00125 /*
00126  * readPZ: read a SAC-format pole-zero file.
00127  * Arguments: pzfile: the name of the pole-zero file to read
00128  *               pRS: pointer to the response structure to be filled in
00129  *                    The calling program must allocate the ResponseStruct;
00130  *                    the individual pole and zero structures will be 
00131  *                    allocated here.
00132  *            
00133  *            Pole-zero-gain files must be for input displacement in 
00134  *            nanometers, output in digital counts, poles and zeros of
00135  *            the LaPlace transform, frequency in radians per second.
00136  * returns: 0 on success
00137  *         -1 on out-of-memory error
00138  *         -2 if unable to read or parse the file
00139  *         -3 for invalid arguments
00140  *         -4 error opeing file
00141  */
00142 int readPZ( char *pzfile, ResponseStruct *pRS )
00143 {
00144   FILE *pzFILE;
00145   int retval = 0, status = 0;
00146   int i, nz = 0, np = 0;
00147   char line[MAXLINE], word[21];
00148   enum states {Key, Pole, Zero};
00149   enum states state = Key;
00150   
00151   if ( pzfile == (char *)NULL || strlen(pzfile) == 0)
00152   {
00153     /* empty or missing file name */
00154     return -3;
00155   }
00156   
00157   if ( (pzFILE = fopen(pzfile, "r")) == (FILE *)NULL)
00158   {
00159     /* Error opening file */
00160     return -4;
00161   }
00162 
00163   while ( fgets( line, MAXLINE, pzFILE) != (char *)NULL)
00164   {
00165     /* if (line[strlen(line)-1] == '\n')
00166        line[strlen(line)-1] = '\0';*/
00167     
00168     switch (state)
00169     {
00170     case Key:  /* Looking for next keyword */
00171       if (sscanf(line, "%20s", word) == 0)
00172         continue;
00173       if ( word[0] == '*' || word[0] == '#') continue;  /* a comment line */
00174       
00175       if ( strcmp(word, "CONSTANT") == 0)
00176       {
00177         if (sscanf(line, "%*s %lf", &pRS->dGain) == 0 )
00178         {
00179           retval = -2;
00180           goto abort;
00181         }
00182         status |= 1;  /* Found the constant or gain */
00183       }
00184       else if ( strcmp(word, "ZEROS") == 0)
00185       {
00186         if (sscanf(line, "%*s %d", &pRS->iNumZeros) == 0 || pRS->iNumZeros < 0)
00187         {
00188           /* invalid or missing number after ZEROS keyword */
00189           retval = -2;
00190           goto abort;
00191         }
00192         if (pRS->iNumZeros > 0)
00193         {
00194           if ( (pRS->Zeros = (PZNum *)malloc(pRS->iNumZeros * sizeof(PZNum))) 
00195                == (PZNum *)0 )
00196           {
00197             retval = -1;
00198             goto abort;
00199           }
00200           for (i = 0; i < pRS->iNumZeros; i++)
00201           {
00202             pRS->Zeros[i].dReal = 0.0;
00203             pRS->Zeros[i].dImag = 0.0;
00204           }
00205                     
00206           state = Zero;  /* There are some zeros; go find them */
00207           continue;
00208         }
00209         status |= 2;  /* Got the number of zeros: none */
00210       }
00211       else if ( strcmp(word, "POLES") == 0)
00212       {
00213         if (sscanf(line, "%*s %d", &pRS->iNumPoles) == 0 || pRS->iNumPoles < 0)
00214         {
00215           /* invalid or missing number after POLES keyword */
00216           retval = -2;
00217           goto abort;
00218         }
00219         if (pRS->iNumPoles > 0)
00220         {
00221           if ( (pRS->Poles = (PZNum *)malloc(pRS->iNumPoles * sizeof(PZNum))) 
00222                == (PZNum *)0 )
00223           {
00224             retval = -1;
00225             goto abort;
00226           }
00227           for (i = 0; i < pRS->iNumPoles; i++)
00228           {
00229             pRS->Poles[i].dReal = 0.0;
00230             pRS->Poles[i].dImag = 0.0;
00231           }
00232           
00233           state = Pole;  /* There are some poles; go find them */
00234           continue;
00235         }
00236         status |= 4;  /* Got the number of poles: none */
00237       }
00238       else
00239       {
00240         /* Invalid keyword */
00241         retval = -2;
00242         goto abort;
00243       }
00244       break;
00245     case Zero:
00246       /* Looking for Zeros */
00247       if (nz >= pRS->iNumZeros)
00248       {
00249         /* Too many zeros! */
00250         retval = -2;
00251         goto abort;
00252       }
00253       if (sscanf(line, "%lf %lf", 
00254                  &pRS->Zeros[nz].dReal, &pRS->Zeros[nz].dImag) != 2)
00255       {
00256         /* Couldn't read a line of zeros */
00257         retval = -2;
00258         goto abort;
00259       }
00260       if (++nz == pRS->iNumZeros)
00261       {
00262         /* Found all the zeros we expected */
00263         status |= 2;
00264         state = Key;
00265         continue;
00266       }
00267       break;
00268     case Pole:
00269       /* Looking for poles */
00270       if (np >= pRS->iNumPoles)
00271       {
00272         /* Too many poles! */
00273         retval = -2;
00274         goto abort;
00275       }
00276       if (sscanf(line, "%lf %lf", 
00277                  &pRS->Poles[np].dReal, &pRS->Poles[np].dImag) != 2)
00278       {
00279         /* Couldn't read a line of poles */
00280         retval = -2;
00281         goto abort;
00282       }
00283       if (++np == pRS->iNumPoles)
00284       {
00285         /* Found all the poles we expected */
00286         status |= 4;
00287         state = Key;
00288         continue;
00289       }
00290       break;
00291     }
00292   }
00293 
00294   if (status != 7)
00295   {
00296     /* One of the keywords was missing */
00297     retval = -2;
00298   }
00299   
00300  abort:
00301   if (retval != 0)
00302   {
00303     /* Something went wrong; clean up the mess */
00304     cleanPZ(pRS);
00305   }
00306   fclose(pzFILE);
00307   
00308   return retval;
00309 }
00310 
00311 
00312 /*
00313  * ftaper: produce a cosine taper between unity (beyond fon) and zero
00314  *        (beyond foff). The cosine taper is between fon and foff.
00315  * Arguments: freq: the frequency at which the taper value is desired
00316  *             fon: the unity end of the taper
00317  *            foff: the zero end of the taper
00318  *    if fon and foff are equal, then taper returns 1.0, the all-pass filter.
00319  * returns: the value of the taper
00320  */
00321 double ftaper(double freq, double fon, double foff)
00322 {
00323   double t, pi = PI;
00324   
00325   if (fon > foff)
00326   {   /* high-pass taper */
00327     if (freq < foff)
00328       t = 0.0;
00329     else if (freq > fon)
00330       t = 1.0;
00331     else
00332       t = 0.5 * (1.0 - cos(pi * (freq - foff) / (fon - foff)));
00333   }
00334   else if (fon < foff)
00335   {   /* low-pass case */
00336     if (freq < fon)
00337       t = 1.0;
00338     else if (freq > foff)
00339       t = 0.0;
00340     else
00341       t = 0.5 * (1.0 + cos(pi * (freq - fon) / (foff - fon)));
00342   }
00343   else
00344     t = 1.0;
00345   
00346   return t;
00347 }
00348 
00349 /*
00350  * convertWave: converts a waveform (time series) from its original response
00351  *              function to a new response function. This conversion is done
00352  *              in the frequency domain. The frequency response of the 
00353  *              transfer function may be tapered. The input data will be
00354  *              padded in the time-domain. The amount of padding is determined
00355  *              automatically unless the user provides her own pad length.
00356  * Arguments: input: array of data for preocessing
00357  *             npts: number of data points to process
00358  *           deltat: time interval between samples, in seconds
00359  *           origRS: structure defining process that generated the input data
00360  *                   that is, the response function to be removed
00361  *          finalRS: structure defining desired response function
00362  *             freq: array of four frequencies (f0, f1, f2, f3) defining the
00363  *                   taper to be applied to the frequency response function
00364  *                   before it is convolved with the data. Below f0 and above
00365  *                   f3, the taper is 0; between f2 and f3 the taper is 1;
00366  *                   between f0-f1 and f2-f3 is a cosine taper.
00367  *            retFD: flag to return result in frequency-domain (if retFD == 1)
00368  *                   or in time-domain (if retFD == 0)
00369  *                   If the output is to stay in the frequency domain,
00370  *                   be sure you understand how the results are laid out. 
00371  *                   See the comments in the FFT package: currently sing.c
00372  *           padlen: The pad length to be applied to data before transforming
00373  *                   to frequency domain. If padlen < 0, pad length will be
00374  *                   estimated here and the value chosen will be returned
00375  *                   in this return-value parameter.
00376  *             nfft: The size of the FFT chosen, based on npts + *padlen
00377  *                   If the returned value of nfft + padlen is less than
00378  *                   npts, then convertWave had to effectively truncate the 
00379  *                   raw trace in order to fit the processed trace in
00380  *                   the limit of outBufLen.
00381  *           output: array of values output from the conversion
00382  *                   This array must be allocated by the caller. 
00383  *        outBufLen: size of `output' array.
00384  *             work: a work array that must be allocated by the caller.
00385  *                   Its size must be outBufLen+2
00386  *          workFFT: a work array needed by fft99. 
00387  *                   Its size must be outBufLen+1
00388  *
00389  * Returns: 0 on success
00390  *         -1 on out-of-memory errors                   
00391  *         -2 on too-small impulse response
00392  *         -3 on invalid arguments
00393  *         -4 on FFT error
00394  */
00395 
00396 int convertWave(double input[], long npts, double deltat, 
00397                 ResponseStruct *origRS, ResponseStruct *finalRS, 
00398                 double freq[4], int retFD, long *pPadLen, long *pnfft, 
00399                 double output[], long outBufLen, double *work, double *workFFT)
00400 {
00401   ResponseStruct rs;  /* the combined response function */
00402   FACT *pfact;
00403   long i, ii, retval = 0;
00404   int nz = 0, np = 0;
00405   long nfft, nfreq, trial_nfft;
00406   double delfreq, tpr, dre, dim, f, delomg, omega;
00407   double *fre, *fim;
00408   
00409   /* Validate arguments */
00410   if (origRS == (ResponseStruct *)NULL || finalRS == (ResponseStruct *)NULL ||
00411       npts < 2 || deltat <= 0.0 || freq == (double *)NULL || outBufLen < npts)
00412   {
00413     return -3;
00414   }
00415   if (freq[0] > freq[1] || freq[1] >= freq[2] || freq[2] > freq[3] ||
00416       origRS->dGain == 0.0 || finalRS->dGain == 0.0)
00417   {
00418     return -3;
00419   }
00420   
00421 
00422   /* Combine the response functions into one: finalRS / origRS */
00423   rs.dGain = finalRS->dGain / origRS->dGain;
00424   rs.iNumPoles = finalRS->iNumPoles + origRS->iNumZeros;
00425   rs.iNumZeros = finalRS->iNumZeros + origRS->iNumPoles;
00426   if ( (rs.Poles = (PZNum *)malloc(sizeof(PZNum) * rs.iNumPoles)) == 
00427        (PZNum *)NULL ||
00428        (rs.Zeros = (PZNum *)malloc(sizeof(PZNum) *rs.iNumZeros)) ==
00429        (PZNum *)NULL)
00430   {
00431     retval = -1;
00432     goto exit;
00433   }
00434   /* Copy the poles and zeros, using structure copy */
00435   for (i = 0; i < origRS->iNumPoles; i++)
00436     rs.Zeros[nz++] = origRS->Poles[i];
00437   for (i = 0; i < origRS->iNumZeros; i++)
00438     rs.Poles[np++] = origRS->Zeros[i];
00439   for (i = 0; i < finalRS->iNumPoles; i++)
00440     rs.Poles[np++] = finalRS->Poles[i];
00441   for (i = 0; i < finalRS->iNumZeros; i++)
00442     rs.Zeros[nz++] = finalRS->Zeros[i];
00443   
00444   if (Debug & TR_DBG_PZG)
00445   {
00446     printf("Input response function: gain %10.3e\n", origRS->dGain);
00447     printf("Poles: %d\n", origRS->iNumPoles);
00448     for (i = 0; i < origRS->iNumPoles; i++)
00449       printf("%10.3e   %10.3e\n", origRS->Poles[i].dReal, 
00450              origRS->Poles[i].dImag);
00451     printf("\nZeros: %d\n", origRS->iNumZeros);
00452     for (i = 0; i < origRS->iNumZeros; i++)
00453       printf("%10.3e   %10.3e\n", origRS->Zeros[i].dReal, 
00454              origRS->Zeros[i].dImag);
00455     printf("\nOutput response function: gain %10.3e\n", finalRS->dGain);
00456     printf("Poles: %d\n", finalRS->iNumPoles);
00457     for (i = 0; i < finalRS->iNumPoles; i++)
00458       printf("%10.3e   %10.3e\n", finalRS->Poles[i].dReal, 
00459              finalRS->Poles[i].dImag);
00460     printf("\nZeros: %d\n", finalRS->iNumZeros);
00461     for (i = 0; i < finalRS->iNumZeros; i++)
00462       printf("%10.3e   %10.3e\n", finalRS->Zeros[i].dReal, 
00463              finalRS->Zeros[i].dImag);
00464   }
00465 
00466   /* Determine how much padding we need, unless the caller told us. */
00467   if ( *pPadLen < 0)
00468   {
00469     if ( (*pPadLen = respLen( &rs, deltat, freq)) < 0)
00470     {
00471       /* Some error occured */
00472       retval = *pPadLen;
00473       if (Debug)
00474         printf("\nrespLen error: %ld\n", *pPadLen);
00475       goto exit;
00476     }
00477     if (Debug)
00478       printf("estimated pad length: %ld\n", *pPadLen);
00479   }
00480   
00481   /* Find a convenient FFT size for our data plus padding that will fit in
00482    *  outBuf */
00483   trial_nfft = *pPadLen + npts;
00484   while ( (nfft = prepFFT(trial_nfft, &pfact)) > outBufLen)
00485   {
00486     if (nfft < 0)
00487     {
00488       /* Out of memory */
00489       retval = nfft;
00490       goto exit;
00491     }
00492     trial_nfft -= 100;  /* Try a liitle bit smaller */
00493   }
00494   if (nfft - *pPadLen < npts)
00495     npts = nfft - *pPadLen;   /* Chop some off the end if it won't fit */
00496   /* We aren't passing this new value of npts back to our caller; she'll
00497    * have to figure it out from the values of pPadLen and nfft passed back */
00498 
00499   nfreq = nfft / 2 + 1;
00500   fre = work;
00501   fim = work + 1 + outBufLen/2;
00502 
00503   for (i = 0; i < npts; i++)
00504     output[i] = input[i];
00505   /* Fill the remainder of output buffer with zeros.                 *
00506    * For fft99, we must fill two slots past the normal end of output *
00507    * buffer; this space must be allocated in initBufs().             */
00508   for (i = outBufLen + FFT_EXTRA; i < nfft; i++)
00509     output[i] = 0.0;
00510   
00511   /* Transform the padded data into frequency domain */
00512   fft991(output, workFFT, pfact->trigs, pfact->ifax, 1, nfft, nfft, 1, -1);
00513 
00514   response(nfft, deltat, &rs, fre, fim);
00515   
00516   /* Convolve the tapered frequency response with the data. Since we  *
00517    * are in the frequency domain, convolution becomes `multiply'.     *
00518    * We skip the zero-frequency part; this only affects the mean      *
00519    * of the data, which should have been removed long ago.            */
00520   delfreq = 1.0 / (nfft * deltat);
00521   output[0] = 0.0;   /* Remove the mean, if there is any */
00522   for (i = 1; i < nfreq - 1; i++)
00523   {
00524     ii = i+i;
00525     f = i * delfreq;
00526     tpr = ftaper(f, freq[1], freq[0]) * ftaper(f, freq[2], freq[3]);
00527     dre = output[ii];   /* Real part of transformed data */
00528     dim = output[ii+1]; /* Imaginary part of transformed data */
00529     output[ii] = (dre * fre[i] - dim * fim[i]) * tpr;
00530     output[ii+1] = (dre * fim[i] + dim * fre[i]) * tpr;
00531   }
00532   f = i * delfreq;
00533   tpr = ftaper(f, freq[1], freq[0]) * ftaper(f, freq[2], freq[3]);
00534   dre = output[nfft];  /* Real part of transformed data; imaginary part is 0 */
00535   output[nfft] = dre * fre[i] * tpr;
00536 
00537   if (Debug & TR_DBG_ARS)
00538   {
00539     delomg = 2.0 * PI / (nfft * deltat);
00540     printf("    omega        tapered\n");
00541     printf("%4ld  %10.3e  %10.3e  %10.3e\n", 0L, 0.0, 0.0, 0.0);
00542     for (i = 1; i < nfreq - 1; i++ )
00543     {
00544       omega = i * delomg;
00545       printf("%4ld  %10.3e  %10.3e  %10.3e\n", i, omega, fre[i], fim[i]);
00546     }
00547   }
00548   
00549   /* If the user wants data in the time domain, transform it now */
00550   if (!retFD)
00551     fft991(output, workFFT, pfact->trigs, pfact->ifax, 1, nfft, nfft, 1, +1);
00552 
00553   /* We're done! Give the data back to the caller */
00554   *pnfft = nfft;
00555 
00556   exit:
00557     /* Clean up before we quit: free any memory that was allocated */
00558   cleanPZ( &rs );
00559   return retval;
00560 }
00561 
00562 /* NFFT_TEST is the size of the iFFT used for estimated the pad length. *
00563  * The value of NFFT_TEST MUST be a multiple of 2, 3, and/or 5.         */
00564 #define NFFT_TEST 1024
00565 /*
00566  * ALIAS_CUTOFF is the fraction of the maximum of the impulse response
00567  * function used for testing pad length. The length (in points) of the
00568  * impulse response that is above this threshold is the pad length
00569  */
00570 #define ALIAS_CUTOFF 0.01
00571 
00572 /*
00573  * respLen: estimate the length of the impulse response (the number of
00574  *          points where it is greater than some throshold) so we know
00575  *          how much padding we need for the convolution.      
00576  *          This is a trial algorithm that may or may not work.
00577  *          We assume the the impulse response looks something like a
00578  *          broadened and shifted impulse. We asssume that the width of
00579  *          its peak is independent of the number of points used in
00580  *          this trial FFT, as long as the peak isn't too broad.  
00581  *  Returns: the length of the peak (in data points) on success
00582  *           NFFT_TEST when impulse response never drops below threshold
00583  *          -1 when out of memory
00584  *          -2 when impulse response is too small to analyze
00585  *          -4 on FFT failure
00586  *          -NFFT_TEST when search tp left of peak finds drop-off
00587  *           but search to right doesn't find drop-off: logic error.
00588  *    Exits if NFFT_TEST is not a multiple of powers of 2, 3, 5
00589  *          That would be a coding error only.
00590  */
00591 int respLen( ResponseStruct *rs, double deltat, double freq[4])
00592 {
00593   FACT *pfact;
00594   double fre[NFFT_TEST/2+1], fim[NFFT_TEST/2+1], data[NFFT_TEST + 2];
00595   double work[NFFT_TEST+1];
00596   long i, ii, nf_test, nfreq;
00597   double imp_max, thresh, delfreq, f, tpr, delomg, omega;
00598   long imax, left_lim, right_lim;
00599   
00600   /*
00601    * Determine how much padding we need.
00602    * We do this by computing the frequency response for a small number of
00603    * points, getting the impulse response (iFFT of freq. response),
00604    * and measuring the length the non-zero part of the response.
00605    */
00606   if ( (nf_test = prepFFT( NFFT_TEST, &pfact)) < 0)
00607   {
00608     /* Out of memory */
00609     return -1;
00610   }
00611   else if (nf_test != NFFT_TEST)
00612   {
00613     fprintf(stderr, "respLen fatal error: NFFT_TEST (%d) not factorable by 2, 3, 5\n",
00614             NFFT_TEST);
00615     return -1;
00616   }
00617   nfreq = nf_test/2+1;
00618   
00619   response(nf_test, deltat, rs, fre, fim);
00620   
00621   delfreq = 1.0 / (nf_test * deltat);
00622 
00623   if (Debug & TR_DBG_TRS)
00624   {   /* Print the values while they are available */
00625     delomg = 2.0 * PI / (nf_test * deltat);
00626     printf("\nTest response function in Frequency Domain\n");
00627     printf("  i      omega             raw                   tapered\n");
00628     f = 0.0;
00629     tpr = ftaper(f, freq[1], freq[0]);
00630     data[0] = fre[0] * tpr;
00631     data[1] = 0.0;
00632     omega = 0.001 * delomg;
00633     printf("%4ld  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n", 0L, omega, fre[0],
00634            0.0, data[0], 0.0);
00635     for (i = 1; i < nfreq - 1; i++ )
00636     {
00637       ii = i+i;
00638       f = i * delfreq;
00639       tpr = ftaper(f, freq[1], freq[0]) * ftaper(f, freq[2], freq[3]);
00640       data[ii] = fre[i] * tpr;
00641       data[ii+1] = fim[i] * tpr;
00642       omega = i * delomg;
00643       printf("%4ld  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n", i, omega, 
00644              fre[i], fim[i], data[ii], data[ii+1]);
00645     }
00646     data[nf_test] = 0.0;
00647     data[nf_test+1] = 0.0;
00648   }
00649   else
00650   {  /* Don't print, just calculate quickly with no /if/ inside the loops */
00651     f = 0.0;
00652     tpr = ftaper(f, freq[1], freq[0]);
00653     data[0] = fre[0] * tpr;
00654     data[1] = 0.0;
00655     for (i = 1; i < nfreq - 1; i++ )
00656     {
00657       ii = i+i;
00658       f = i * delfreq;
00659       tpr = ftaper(f, freq[1], freq[0]) * ftaper(f, freq[2], freq[3]);
00660       data[ii] = fre[i] * tpr;
00661       data[ii+1] = fim[i] * tpr;
00662     }
00663     data[nf_test] = 0.0;
00664     data[nf_test+1] = 0.0;
00665   }
00666   
00667   /* Transform test response function into time domain */
00668   fft991(data, work, pfact->trigs, pfact->ifax, 1, nf_test, nf_test, 1, +1);
00669 
00670   
00671   if (Debug & TR_DBG_TRS)
00672   {
00673     printf("\nTest response function in TD\n");
00674     for (i = 0; i < nf_test; i++)
00675     {
00676       printf("%5ld  %10.3e\n", i, data[i]);
00677     }
00678   }
00679  
00680   /*
00681    * We assume that the impulse response has a peak somewhere and 
00682    * falls off on both sides of that peak, allowing for wrap-around.
00683    * So we locate the peak and set a threshold of ALIAS_CUTOFF (currently 1%)
00684    * of the peak. If the response does not drop off below this threshold
00685    * for any values of i, then we're screwed. Some other measure will have
00686    * to be dreamed for preventing time-domain aliasing in our convolution.
00687    */
00688   imp_max = 0.0;
00689   for (i = 0; i < nf_test; i++)
00690   {
00691     if (fabs(data[i]) > imp_max)
00692     {
00693       imp_max = fabs(data[i]);  /* the extremal value */
00694       imax = i;                 /* the location of the extremum */
00695     }
00696   }
00697   
00698   if (imp_max > 0.001 * rs->dGain)
00699   {   /* make sure the response isn't too small to measure */
00700     thresh = imp_max * ALIAS_CUTOFF;
00701     /* Start searching to the left from the peak to find where we drop   *
00702      * below the threshold. This point will be left_lim; we may have to  *
00703      * wrap the search around to the right end. If this earch continues  *
00704      * all the way back to the peak, then the peak is too broad and this *
00705      * algorithm fails.                                                  */
00706     right_lim = left_lim = imax;
00707     for (i = imax; i >= 0; i--)
00708     {
00709       if ( fabs(data[i]) < thresh)
00710       {
00711         left_lim = i;
00712         break;
00713       }
00714     }
00715     if (left_lim == imax) /* Didn't find it; wrap around */
00716     {
00717       for (i = nf_test - 1; i > imax; i--)
00718       {
00719         if ( fabs(data[i]) < thresh)
00720         {
00721           left_lim = i;
00722           break;
00723         }
00724       }
00725       if (left_lim == imax)
00726       {   /* Still didn't find it; peak doesn't drop of anywhere */
00727         return NFFT_TEST;
00728       }
00729     }
00730     /* Now search to the right for right_lim; may have to wrap around */
00731     for (i = imax; i < nf_test; i++)
00732     {
00733       if ( fabs(data[i]) < thresh)
00734       {
00735         right_lim = i;
00736         break;
00737       }
00738     }
00739     if (right_lim == imax)
00740     {  /* Didn't find it; wrap around to the left side */
00741       for (i = 0; i < imax; i++)
00742       {
00743         if ( fabs(data[i]) < thresh)
00744         {
00745           right_lim = i;
00746           break;
00747         }
00748       }
00749       if (right_lim == imax)
00750       {   /* How come we found left_lim but not right_lim: shouldn't happen! */
00751         return -NFFT_TEST;
00752       }
00753     }
00754     if (left_lim < right_lim)
00755       return (right_lim - left_lim + 1);
00756     else
00757       return (right_lim - left_lim + nf_test + 1);
00758   }
00759   /* Impulse response is too small to analyze */
00760   return -2;
00761 }
00762 
00763 /*
00764  * pzCancel: Remove cancelling pole-zero pairs from a response structure.
00765  *           Search for pairs of poles and zeros whose real and imaginary
00766  *           parts are within `tol' of each other. Remove any such pairs.
00767  *           This will remove useless computations from the calculation
00768  *           of the frequency response function in response().
00769  */
00770 void pzCancel(ResponseStruct *rs, double tol)
00771 {
00772   int ip, iz, mz;
00773   
00774   if (rs->iNumZeros == 0 || rs->iNumPoles == 0)
00775     return;
00776   
00777   for (ip = 0; ip < rs->iNumPoles; ip++)
00778   {
00779     mz = -1;
00780     for (iz = 0; iz < rs->iNumZeros; iz++)
00781     {
00782       if (fabs(rs->Poles[ip].dReal - rs->Zeros[iz].dReal) < tol &&
00783           fabs(rs->Poles[ip].dImag - rs->Zeros[iz].dImag) < tol)
00784       {
00785         mz = iz;
00786         break;
00787       }
00788     }
00789     if (mz != -1)
00790     {
00791       /* Found a match; remove the pole and the zero; move the other *
00792        * poles and zeros into those empty slots.                     */
00793       drop(rs->Poles, &rs->iNumPoles, ip);
00794       drop(rs->Zeros, &rs->iNumZeros, mz);
00795       ip--;  /* We have to look at this pole slot again */
00796       iz--;  /* and also this zero slot */
00797     }
00798   }
00799   return;
00800 }
00801 
00802 /*
00803  * drop: remove the `ipz' pole or zero from the PZNum array;
00804  *       update the pNumPZ counter.
00805  */
00806 static void drop(PZNum *pPZ, int *pNumPZ, int ipz)
00807 {
00808   int i;
00809   
00810   (*pNumPZ)--;
00811   for (i = ipz; i < *pNumPZ; i++)
00812   {
00813     pPZ[i] = pPZ[i+1];   /* structure copy */
00814   }
00815   return;
00816 }
00817 
00818 
00819 /* 
00820  * taper: Apply a cosine taper to a data series.
00821  * Arguments:  data: array of data to be tapered
00822  *             npts: number of points in data array
00823  *             tLen: width of taper (at each end) in number of points
00824  *      The end of the taper (where it has 0 value) is assumed to occur 
00825  *      at the first point before the start and after the end of the data;
00826  *      these data points are not modified here. There are `tLen' points
00827  *      that have taper values between 0 and 1 exclusive at each end of
00828  *      the data.
00829  *      If tLen is less then half of npts, this function returns
00830  *      silently without applying the taper.
00831  */
00832 void taper(double *data, long npts, long tLen)
00833 {
00834   long i, jb, je;
00835   double tap, omega;
00836   
00837   tLen++;
00838   if (tLen < 2 || tLen > npts / 2)
00839     return;
00840   
00841   omega = PI / tLen;
00842   
00843   for (i = 1; i < tLen; i++)
00844   {
00845     jb = i - 1;
00846     je = npts - i;
00847     tap = 0.5 * (1.0 - cos(omega * i));
00848     data[jb] *= tap;
00849     data[je] *= tap;
00850   }
00851   return;
00852 }
00853 
00854 
00855 void deMean( double *data, long npts, double *mean)
00856 {
00857   long i;
00858   double sum = 0.0;
00859   
00860   if (npts < 1)
00861   {
00862     *mean = 0.0;
00863     return;
00864   }
00865   
00866   for (i = 0; i < npts; i++)
00867     sum += data[i];
00868   
00869   sum /= (double)npts;
00870   for (i = 0; i < npts; i++)
00871     data[i] -= sum;
00872   
00873   *mean = sum;
00874   return;
00875 }
00876 
00877   
00878 void cleanPZ( ResponseStruct *pRS)
00879 {
00880   if (pRS->Zeros != (PZNum *)NULL)
00881   {
00882     free( pRS->Zeros );
00883     pRS->Zeros = (PZNum *)NULL;
00884   }
00885   if (pRS->Poles != (PZNum *)NULL)
00886   {
00887     free( pRS->Poles );
00888     pRS->Poles = (PZNum *)NULL;
00889   }
00890   pRS->iNumPoles = 0;
00891   pRS->iNumZeros = 0;
00892   return;
00893 }
00894 
00895 void transferDebug( int level )
00896 {
00897   Debug = level;
00898   return;
00899 }

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