Using SAC

Using the SAC Library

Overview

In addition to being able to read and write SAC data files in one's own C or FORTRAN programs (see SAC Reading and Writing Routines ), one can use many of SAC's data-processing routines in stand-alone codes if one uses specific flags in the compiling stage and the SAC libraries in the linking stage. These libraries libsac.a and libsacio.a are in ${SACHOME}/lib. The .f and .c files, as well as the sample input files, are in ${SACHOME}/doc/examples.

Compiling

To ease the requirements for compilation and linking, a helper script is provided, ${SACHOME}/bin/sac-config, which should output the necessary flags and libraries. Try the following:

gcc -o program source.c `sac-config --cflags --libs sac`

f77 -o program source.f `sac-config --cflags --libs sac`

Filtering

Example use of filtering with xapiir() to Apply an Infinite Impulse Recursive filter

Fortran

      program filter
      implicit none

      include "sacf.h"

!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=1000)

!     Define the Data Array of size MAX
      real yarray, xarray
      dimension yarray(MAX)

!     Declare Variables used in the rsac1() subroutine
      real beg, delta
      integer nlen
      character*20 KNAME
      integer nerr

!     Define variables used in the filtering routine
      real *8 low, high
      real *8 transition_bandwidth, attenuation
      real *8 delta_d
      integer order, passes

      kname = 'filterf_in.sac'

      call rsac1(kname, yarray, nlen, beg, delta, MAX, nerr)
      delta_d = delta

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

      low    = 0.10
      high   = 1.00
      passes = 2
      order  = 4
      transition_bandwidth = 0.0
      attenuation = 0.0
!     Call xapiir ( Apply a IIR Filter )
!        - yarray - Original Data
!        - nlen   - Number of points in yarray
!        - proto  - Prototype of Filter
!                 - SAC_FILTER_BUTTERWORK        - Butterworth
!                 - SAC_FILTER_BESSEL            - Bessel
!                 - SAC_FILTER_CHEBYSHEV_TYPE_I  - Chebyshev Type I
!                 - SAC_FILTER_CHEBYSHEV_TYPE_II - Chebyshev Type II
!        - transition_bandwidth (Only for Chebyshev Filter)
!                 - Bandwidth as a fraction of the lowpass prototype
!                   cutoff frequency
!        - attenuation (Only for Chebyshev Filter)
!                 - Attenuation factor, equals amplitude reached at
!                   stopband egde
!        - order  - Number of poles or order of the analog prototype
!                   4 - 5 should be ample
!                   Cannot exceed 10
!        - type   - Type of Filter
!                 - SAC_FILTER_BANDPASS
!                 - SAC_FILTER_BANDREJECT
!                 - SAC_FILTER_LOWPASS
!                 - SAC_FILTER_HIGHPASS
!        - low    - Low Frequency Cutoff [ Hertz ]
!                   Ignored on SAC_FILTER_LOWPASS
!        - high   - High Frequency Cutoff [ Hertz ]
!                   Ignored on SAC_FILTER_HIGHPASS
!        - delta  - Sampling Interval [ seconds ]
!        - passes - Number of passes
!                 - 1 Forward filter only
!                 - 2 Forward and reverse (i.e. zero-phase) filtering
!
      call xapiir(yarray, nlen,
     +            SAC_BUTTERWORTH,
     +            transition_bandwidth, attenuation,
     +            order,
     +            SAC_BANDPASS,
     +            low, high, delta_d, passes)

!     Do more processing ....

      xarray = 0
      kname='filterf_out.sac'
!     Write the SAC file kname
!       - kname holds the name of the file to be written
!       - xdata Input Time Data
!       - yfunc Input Amplitude Data
!       - nerr Error return Flag
      call wsac0(kname,xarray,yarray,nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',kname
          call exit(-1)
      endif

      call exit(0)

      end program filter

C

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include "sac.h"
#include "sacio.h"

#define MAX 1000

int
main(int argc, char *argv[]) {

    /* Local variables */
    double low, high, attenuation, transition_bandwidth;;
    int nlen, nerr, max;

    float beg, delta;
    double delta_d;
    char *kname;
    int order;

    int passes;
    float yarray[1000];
    float xarray[1];

    max = MAX;

    /*     Define the Maximum size of the data Array
     * Filter Prototypes
     * Filter Types
     *     Define the Data Array of size MAX
     *     Declare Variables used in the rsac1() subroutine
     *     Define variables used in the filtering routine
     */

    kname = strdup("filterc_in.sac");
    rsac1(kname, yarray, &nlen, &beg, &delta, &max, &nerr, SAC_STRING_LENGTH);
    delta_d = delta;
    if (nerr != 0) {
      fprintf(stderr, "Error reading in file: %s\n", kname);
      exit(-1);
    }

    low    = 0.10;
    high   = 1.00;
    passes = 2;
    order  = 4;
    transition_bandwidth = 0.0;
    attenuation = 0.0;

    /*     Call xapiir ( Apply a IIR Filter )
     *        - yarray - Original Data
     *        - nlen   - Number of points in yarray
     *        - proto  - Prototype of Filter
     *                 - SAC_FILTER_BUTTERWORK        - Butterworth
     *                 - SAC_FILTER_BESSEL            - Bessel
     *                 - SAC_FILTER_CHEBYSHEV_TYPE_I  - Chebyshev Type I
     *                 - SAC_FILTER_CHEBYSHEV_TYPE_II - Chebyshev Type II
     *        - transition_bandwidth (Only for Chebyshev Filter)
     *                 - Bandwidth as a fraction of the lowpass prototype
     *                   cutoff frequency
     *        - attenuation (Only for Chebyshev Filter)
     *                 - Attenuation factor, equals amplitude reached at
     *                   stopband egde
     *        - order  - Number of poles or order of the analog prototype
     *                   4 - 5 should be ample
     *                   Cannot exceed 10
     *        - type   - Type of Filter
     *                 - SAC_FILTER_BANDPASS
     *                 - SAC_FILTER_BANDREJECT
     *                 - SAC_FILTER_LOWPASS
     *                 - SAC_FILTER_HIGHPASS
     *        - low    - Low Frequency Cutoff [ Hertz ]
     *                   Ignored on SAC_FILTER_LOWPASS
     *        - high   - High Frequency Cutoff [ Hertz ]
     *                   Ignored on SAC_FILTER_HIGHPASS
     *        - delta  - Sampling Interval [ seconds ]
     *        - passes - Number of passes
     *                 - 1 Forward filter only
     *                 - 2 Forward and reverse (i.e. zero-phase) filtering
     */
    xapiir(yarray, nlen, SAC_BUTTERWORTH,
           transition_bandwidth, attenuation,
           order,
           SAC_BANDPASS,
           low, high,
           delta_d, passes);

    /*     Do more processing ....  */
    xarray[0] = 0;
    kname = strdup("filterc_out.sac");
    wsac0(kname, xarray, yarray, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
      fprintf(stderr, "Error writing out file: %s\n", kname);
      exit(-1);
    }


    return 0;
}

Cut Data

Fortran

! Demonstrate applying a cut to a seismogram using the SAC library
program cut_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer :: npts, npts_cut, nstart, nstop
    integer :: nfillb, nfille, nerr, cuterr
    real    :: data(nmax), cut_data(nmax)
    real    :: beg, dt, begin_cut, end_cut

    ! Read in the data file
    call rsac1('raw.sac', data, npts, beg, dt, nmax, nerr)

    ! Set up cut parameters
    begin_cut = 10.0     ! Begin time of cut
    end_cut = 15.0       ! End time of cut
    cuterr = 3         ! fill with zeros if the window is too large

!   Call cut_define
!    - begin_cut    - Begin time for cut
!    - dt           - Sample rate of data
!    - end_cut      - End time for cut
!    - npts_cut     - Number of points in data after cutting
    call cut_define(begin_cut, dt, end_cut, npts_cut)

!   Call cut_define_check
!    - begin_cut - Begin time for cut
!    - end_cut   - End time for cut
!    - npts      - Number of points in data
!    - cuterr    - How to handle cuts outside the length of the trace. Three possible values:
!                - CUT_FATAL = 1 throws an error if the cut window is too large
!                - CUT_USEBE = 2 use the b and e values of the trace if the cut window is too large
!                - CUT_FILLZ = 3 fills with zeros if the cut windows is too large
!    - nstart    - Number of points corresponding to begin_cut
!    - nstop     - Number of points corresponding to end_cut
!    - nfillb    - Number of points filled before begin_time
!    - nfille    - Number of pionts filled after end_time
!    - nerr      - Error value returned
    call cut_define_check(begin_cut, end_cut, npts, cuterr, nstart, nstop, nfillb, nfille, nerr)

!   Call cut
!   - data       - Original data to cut
!   - nstart     - Number of points corresponding to begin_cut
!   - nstop      - Number of points corresponding to end_cut
!   - nfillb     - Number of points filled with zeros if begin_time is before data
!   - nfille     - Number of points filled with zeros if end_time is after data
!   - cut_data   - Cut data
    call cut(data, nstart, nstop, nfillb, nfille, cut_data)

    ! write the cut seismogram back to disk
    call wsac1('cut.sac', cut_data(nstart:nstop), npts_cut, begin_cut, dt, nerr)

end program cut_example

Interpolate Data

Fortran

! Demonstrate interpolating a seismogram using the SAC library
program interpolate_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer :: npts, newlen, nerr
    real*4 :: data(nmax)
    real*4 :: beg, eval, dt, tstart, dtnew
    real*4 :: eps = 0.01

    ! Read in the data file
    call rsac1('raw.sac', data, npts, beg, dt, nmax, nerr)

    ! Set up interpolation parameters
    eval = beg + float(npts - 1)*dt ! Ending time of original data
    newlen = npts/2 ! Number of point in interpolated data
    dtnew  = dt/2 ! New sample rate (half of current sample rate)

!   Call interp ( Interpolates the seismogram to a new sample rate )
!    - data   - Original Data
!    - npts   - Number of points in data
!    - interpolated_data   - Interpolated Data
!    - newlen  - Number of points in interpolated data
!    - beg     - Beginning time of original data
!    - eval    - Ending time of original data
!    - dt      - Sample rate of original data
!    - tstart  - Start time of interpolated data
!    - dtnew   - Sample rate of interpolated data
!    - eps     - Machine epsilon precision
    call interp(data, npts, data, newlen, beg, eval, dt, beg, dtnew, eps)

    ! write the interpolated seismogram back to disk
    call wsac1('interp.sac', data(1:newlen), newlen, beg, dtnew, nerr)

end program interpolate_example

Remove Mean

Fortran

! Demonstrate removing the mean from a seismogram using the SAC library
program rmean_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer :: nlen, npts, nerr
    real*4 :: data(nmax)
    real*4 :: beg, dt, mean

    ! Read in the data file
    call rsac1('raw.sac', data, npts, beg, dt, nmax, nerr)

!   Call rmean ( Removes the mean )
!    - data   - Original Data
!    - npts   - Number of points in data
!    - mean   - Mean value of the Original Data
    call rmean(data, npts, mean)

    ! write the seismogram with mean removed back to disk
    call wsac0('rmean.sac', data, data, nerr)

end program rmean_example

Remove Trend

Fortran

! Demonstrate removing the trend from a seismogram using the SAC library
program rtrend_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer :: nlen, npts, nerr
    real*4 :: data(nmax)
    real*4 :: beg, dt, yint, slope

    ! Read in the data file
    call rsac1('raw.sac', data, npts, beg, dt, nmax, nerr)

!   Call rtrend ( Removes the trend )
!    - data   - Original Data
!    - npts   - Number of points in data
!    - yint   - y-intercept of best-fitting line
!    - slope  - Slope of best-fitting line
!    - beg    - Beginning time of original data
!    - dt     - Sample rate of original data
    call rtrend(data, npts, yint, slope, beg, dt)

    ! write the seismogram with trend removed back to disk
    call wsac0('rtrend.sac', data, data, nerr)

end program rtrend_example

Taper Data

Fortran

! Demonstrate applying a taper to a seismogram using the SAC library
program taper_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer :: nlen, npts, nerr, ipts, taper_type
    real*4 :: data(nmax)
    real*4 :: beg, dt, width

    ! Read in the data file
    call rsac1('raw.sac', data, npts, beg, dt, nmax, nerr)

    ! Set up taper parameters
    width = .05     ! Width to taper original data
    taper_type = 2  ! HANNING taper

!   Call taper_width_to_points
!    - width    - Width to taper original data
!    - npts     - Number of points in data
!    - ipts     - Number of points that will be tapered in data
    call taper_width_to_points(width, npts, ipts)

!   Call taper ( Applies a Taper )
!    - data   - Original Data
!    - npts   - Number of points in data
!    - taper_type - Typer of taper: 1 = cosine; 2 = hanning; 3 = hamming
!    - ipts   - Number of points tapered
    call taper(data, npts, taper_type, ipts)

    ! write the seismogram with taper applied back to disk
    call wsac0('taper.sac', data, data, nerr)

end program taper_example

Rotate Data

Fortran

! Demonstrate performing a rotation on a pair of signals using the SAC library
program rotate_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer :: npts, newlen, nerr
    real :: signal1(nmax), signal2(nmax)
    real :: rotated_signal1(nmax), rotated_signal2(nmax)
    real :: beg, dt, angle
    real :: eps = 0.01
    logical :: lnpi, lnpo

    ! Read in the two signals to be rotated
    call rsac1('signal1.sac', signal1, npts, beg, dt, nmax, nerr)
    call rsac1('signal2.sac', signal2, npts, beg, dt, nmax, nerr)

    ! Set up parameters for rotation
    angle = 45.0      ! rotate components 45 degrees clockwise
    lnpi = .true.     ! input signals have "normal" polarity
    lnpo = .true.     ! output signals have "normal polarity

!   Call rotate ( Interpolates the seismogram to a new sample rate )
!   Assumes "normal" polariy is such that the second component leads the
!   first component by 90 degrees in a clockwise rotation
!    - signal1 - First input signal
!    - signal2 - Second input signal
!    - npts    - Number of points in input signal
!    - angle   - Angle of rotation, clockwise from direction of signal1
!    - lnpi    - True if input signals have "normal" polarity
!    - lnpo    - True if output signals have "normal" polarity
!    - rotated_signal1 - First input signal rotated
!    - rotated_signal2 - Second input signal rotated
    call rotate(signal1, signal2, npts, angle, lnpi, lnpo, rotated_signal1, rotated_signal2)

    ! write the seismogram with trend removed back to disk
    call wsac0('rotated1.sac', rotated_signal1, rotated_signal1, nerr)
    call wsac0('rotated2.sac', rotated_signal2, rotated_signal2, nerr)

end program rotate_example

Instrument Removal

Fortran

! Demonstrate removing the instrument response from a seismogram using evresp and the SAC library
program transfer_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer,parameter :: MIN_PERIOD = 60
    integer,parameter :: MAX_PERIOD = 120
    integer,parameter :: FFT_NPTS = 262144
    integer :: npts, newlen, nerr, nfreq, iflag, i, j
    integer :: start_stage, stop_stage, stdio_flag, evresp
    real :: data(nmax)
    real :: f1, f2, f3, f4
    real :: beg, eval, dt, tstart, dtnew
    real :: eps = 0.01
    double precision :: freq(FFT_NPTS), resp(FFT_NPTS*2)
    double precision :: delfrq, x_re(FFT_NPTS), x_im(FFT_NPTS)
    double precision :: resp_re(FFT_NPTS), resp_im(FFT_NPTS)
    real(kind=8), dimension(4) :: F
    character(len=5)  :: unts, sta, cha, net, locid, rtyp
    character(len=18) :: resp_file
    character*200        :: datime
    character*10         :: vbs


    ! Read in the data file
    call rsac1('raw.sac', data, npts, beg, dt, nmax, nerr)

    ! Set up parameters for evresp
    start_stage=-1
    stop_stage = 0
    stdio_flag=0
    unts='DIS'
    vbs='-v'
    rtyp = 'CS'
    do i = 1, FFT_NPTS
      freq(i) = float(i-1)/FFT_NPTS/dt
    enddo
    resp_file = 'RESP.AFI.II.00.BHZ'

    ! Use evresp to generate the instrument respone from RESP file
    iflag = evresp(sta,cha,net,locid,datime,unts,resp_file,freq,FFT_NPTS,resp, &
                    rtyp,vbs,start_stage,stop_stage,stdio_flag)

    ! Set up transfer parameters for deconvolution
    nfreq = 4
    f3 = 1.0/MIN_PERIOD
    f2 = 1.0/MAX_PERIOD
    f1 = f2 * 0.8
    f4 = f3 * 1.2
    F = (/ f1, f2, f3, f4 /)
    delfrq = 1.0/(FFT_NPTS*dt)
    j = 1
    do i = 1, FFT_NPTS
      resp_re(i) = resp(j)
      resp_im(i) = resp(j+1)
      j = j + 2
    enddo
    x_re(1:npts) = data(1:npts)

!   Call ztransfer ( Transfers the response from evalresp to seismogram )
!    - data      - Original Data
!    - npts      - Number of points in data
!    - dt        - Sample rate of original data
!    - x_re      - Real part of TO transfer function
!    - x_im      - Imaginary part of TO transfer function
!    - resp_re   - Real part of FROM transfer function
!    - resp_im   - Imaginary part of FROM transfer function
!    - nfreq     - Number of frequncies in F
!    - FFT_NPTS  - Number of points in frequency domain
!    - detfreq   - Sample interval of number of points in frequency domain
!    - F         - Array containing corner frequncies for deconvolution
    call ztransfer(data, npts, dt, x_re, x_im, resp_re, resp_im, nfreq, FFT_NPTS, delfrq, F)

    ! write the deconvolved seismogram back to disk
    call wsac0('deconvolved.sac', x_re(1:npts), x_re(1:npts), nerr)

end program transfer_example

Envelope Function

Example of taking the envelope of a data series

Fortran

      program envelopef
      implicit none

      include "sacf.h"

!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=1000)

!     Define the Data Array of size MAX
      real yarray, xarray, yenv
      dimension yarray(MAX), yenv(MAX)

!     Declare Variables used in the rsac1() subroutine
      real beg, delta
      integer nlen
      character*20 KNAME
      integer nerr

      kname = 'envelopef_in.sac'

      call rsac1(kname, yarray, nlen, beg, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

!     Call envelope ( Envelope Routine )
!        - nlen   - Number of points in yarray
!        - yarray - Original Data
!        - yenv   - Envelope of the Original Data
!
      call envelope(nlen, yarray, yenv)

!     Do more processing ....

      xarray = 0
      kname='envelopef_out.sac'
!     Write the SAC file kname
!       - kname holds the name of the file to be written
!       - xdata Input Time Data
!       - yfunc Input Amplitude Data
!       - nerr Error return Flag
      call wsac0(kname,xarray,yenv,nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',kname
          call exit(-1)
      endif

      call exit(0)

      end program envelopef

C

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include "sac.h"
#include "sacio.h"

#define MAX 1000

int
main(int argc, char *argv[]) {

    /* Local variables */
    int nlen, nerr, max;

    float beg, delta;
    char *kname;

    float yarray[1000];
    float xarray[1];
    float *yenv;

    max = MAX;

    /*     Define the Maximum size of the data Array
     *     Define the Data Array of size MAX
     *     Declare Variables used in the rsac1() subroutine
     *     Define variables used in the envelope routine
     */

    kname = strdup("envelopec_in.sac");
    rsac1(kname, yarray, &nlen, &beg, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file: %s\n", kname);
      exit(-1);
    }

    /* Allocate space for the envelope of yarray */
    yenv = (float *) malloc(sizeof(float) * nlen);
    if(yenv == NULL) {
      fprintf(stderr, "Error allocating memory for envelope\n");
      exit(-1);
    }

    /*     Call envelope ( Envelope Function )
     *        - nlen   - Number of points in yarray
     *        - yarray - Input array to take the envelope of
     *        - yenv   - Envelope of yarray
     */
    envelope(nlen, yarray, yenv);

    /*     Do more processing ....  */
    xarray[0] = 0;
    kname = strdup("envelopec_out.sac");
    wsac0(kname, xarray, yenv, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
      fprintf(stderr, "Error writing out file: %s\n", kname);
      exit(-1);
    }

    free(yenv);

    return 0;
}

Correlation

Correlation of two time series

Fortran

      program envelopef
      implicit none

      include "sacf.h"

      integer i
!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=4000)

!     Define the Data Array of size MAX
      real yarray1, yarray2, xarray, out, ytmp
      dimension yarray1(MAX), yarray2(MAX), out(MAX*4), ytmp(MAX*4)

!     Declare Variables used in the rsac1() subroutine
      real beg1, beg2, delta, endv
      integer nlen, nlen1, nlen2
      character*30 KNAME
      integer nerr
      real fval

!     Cross Correlation Variables
      integer nwin, wlen, nfft
      character *256 error
      real max_value, max_time

!     Read in the first data file
      kname = 'correlatef_in1.sac'
      call rsac1(kname, yarray1, nlen1, beg1, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

!     Read in the second data file
      kname = 'correlatef_in2.sac'
      call rsac1(kname, yarray2, nlen2, beg2, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

      nlen = nlen1
!
!     If the signals are not the same length, then find the longest
!     signal, make both signals that length by filling the remainder
!     with zeros (pad at the end) and then run them through crscor
!     This should be fixed in upcoming releases and the introduction
!     of a "correlate" function so you do not need to handle the
!     signal length, padding, and window lengths, ...
!

      nwin = 1
      wlen = nlen
      nfft = 0
!     Call crscor ( Cross Correlation )
!        - yarray1 - First  Input array to correlate
!        - yarray2 - Second Input array to correlate
!        - nlen    - Number of points in yarray and yarray2
!        - nwin    - Windows to use in the correlation
!        - wlen    - Length of the windows
!        - type    - Type of Window (SAC_RECTANGLE)
!        - out     - output sequence
!        - nfft    - Length of the output sequence
!        - error   - Error Message
!

      call crscor(yarray1, yarray2, nlen,
     &            nwin, wlen, SAC_RECTANGLE,
     &            ytmp, nfft, error)

      do i = 1,MAX*4
         out(i) = 0.0
      enddo
!
!     out[1 : nlen1 - 1 ] <-- ytmp[ nfft - nlen1 + 2 : nfft  ]
!     out[nlen1 : nlen1 + nlen2 - 1 ] <-- ytmp[ 1 : nlen2 ]
!
      do i = 1,nlen1-1
         out(i) = ytmp(nfft - nlen1 + i + 1)
      enddo
      do i = 1,nlen2
         out(nlen1 + i - 1) = ytmp(i)
      enddo


      nfft = nlen1 + nlen2 - 1
      xarray = 0
      beg1 = -delta * (nlen1 - 1.0) + (beg2 - beg1)
      endv = beg1 + delta * (nfft - 1)

      call setnhv('npts',   nfft,    nerr)
      call setfhv('delta',  delta,   nerr)
      call setlhv('leven',  .true.,  nerr)
      call setfhv('b',      beg1,    nerr)
      call setfhv('e',      endv,    nerr)
      call setihv('iftype', "itime", nerr)
      call setnhv('nzyear', SAC_NUMBER_UNDEFINED, nerr)
      call setnhv('nzhour', SAC_NUMBER_UNDEFINED, nerr)
!     Find the maximum value and time of the correlation function
      max_value = out(1)
      max_time  = 0
      do i = 2,nfft
         if(out(i) > max_value) then
            max_value = out(i)
!           Negative shifts are at the end of the correlation sequence
            if(i > nfft/2) then
               max_time = beg1 + (i - nfft) * delta
            else
               max_time = beg1 + i * delta
            endif
         endif
      enddo

!      call setfhv( 'user0', max_time,  nerr)
!      call setfhv( 'user1', max_value, nerr)


!     Write the SAC file
      kname='correlatef_out1.sac'
      call wsac0(kname, xarray, out, nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',kname
          call exit(-1)
      endif

      call exit(0)

      end program envelopef

C

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <sac.h>
#include <sacio.h>

#define MAX        4000
#define ERROR_MAX  256

int
main(int argc, char *argv[]) {

    /* Local variables */
    int i;
    int nlen, nlen1, nlen2, nerr, max;

    float beg1, beg2, delta, end;
    char *kname;

    float yarray1[MAX], yarray2[MAX], ytmp[MAX*4], xarray[1];
    float *out;
    float max_value, max_time;

    int nwin, wlen, nfft, leven, when;

    char error[ERROR_MAX];

    max = MAX;

    /* Read in the first file  */
    kname = strdup("correlatec_in1.sac");
    rsac1(kname, yarray1, &nlen1, &beg1, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file(%d): %s\n", nerr, kname);
      exit(-1);
    }

    /* Read in the second file  */
    kname = strdup("correlatec_in2.sac");
    rsac1(kname, yarray2, &nlen2, &beg2, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file: %s\n", kname);
      exit(-1);
    }

    nlen = nlen1;
    /**
     *  If the signals are not the same length, then find the longest
     *  signal, make both signals that length by filling the remainder
     *  with zeros (pad at the end) and then run them through crscor
     *  This should be fixed in upcoming releases and the introduction
     *  of a "correlate" function so you do not need to handle the
     *  signal length, padding, and window lengths, ...
     *
     */

    /* Allocate space for the correlation of yarray1 and yarray2 */
    max = next2((2 * nlen) - 1) * 2;
    out = (float *) malloc(sizeof(float) * max);
    if(out == NULL) {
      fprintf(stderr, "Error allocating memory for correlation\n");
      exit(-1);
    }

    /* Set up values for the cross correlation */
    nwin = 1;
    wlen = nlen;
    nfft = 0;

    /*     Call crscor ( Cross Correlation )
     *        - yarray1 - First  Input array to correlate
     *        - yarray2 - Second Input array to correlate
     *        - nlen    - Number of points in yarray and yarray2
     *        - nwin    - Windows to use in the correlation
     *        - wlen    - Length of the windows
     *        - type    - Type of Window (SAC_RECTANGLE)
     *        - out     - output sequence
     *        - nfft    - Length of the output sequence
     *        - error   - Error Message
     *        - err_len - Length of Error Message (on input)
     */
    crscor(yarray1, yarray2, nlen,
           nwin, wlen, SAC_RECTANGLE,
           ytmp, &nfft, error, ERROR_MAX);

    for(i = 0; i < max; i++) {
      out[i] = 0;
    }

    /*
     *     out[0 : nlen1 - 2 ] <-- ytmp[ nfft - nlen1 + 1 : nfft -1 ]
     *     out[nlen1 - 1 : nlen1 + nlen2 - 2 ] <-- ytmp[ 0 : nlen2-1 ]
     */
    for(i = 0; i <= nlen1 - 2; i++) {
      out[i] = ytmp[nfft - nlen1 + i + 1];
    }
    for(i = 0; i <= nlen2 - 1; i++) {
      out[nlen1 + i - 1] = ytmp[i];
    }


    nfft = nlen1 + nlen2 - 1;
    xarray[0] = 0;
    leven = TRUE;
    beg1 = -delta * (nlen1 - 1) + (beg2 - beg1);
    end = beg1 + delta * (nfft - 1);

    setnhv ( "npts",   &nfft,   &nerr, SAC_STRING_LENGTH);
    setfhv ( "delta",  &delta,  &nerr, SAC_STRING_LENGTH);
    setlhv ( "leven",  &leven,  &nerr, SAC_STRING_LENGTH);
    setfhv ( "b",      &beg1,   &nerr, SAC_STRING_LENGTH);
    setfhv ( "e",      &end,    &nerr, SAC_STRING_LENGTH);
    setihv ( "iftype", "itime", &nerr, SAC_STRING_LENGTH, SAC_STRING_LENGTH);
    when = SAC_NUMBER_UNDEFINED;
    setnhv ( "nzyear", &when,   &nerr, SAC_STRING_LENGTH);
    setnhv ( "nzhour", &when,   &nerr, SAC_STRING_LENGTH);

    /* Find the maximum value and time of the correlation function */
    max_value = out[0];
    max_time  = 0;
    for(i = 1; i < nfft; i++) {
      if(out[i] > max_value) {
        max_value = out[i];
        /* Negative shifts are at the end of the correlation sequence */
        if(i > nfft/2) {
          max_time  = (i - nfft) * delta;
        } else {
          max_time  = i * delta;
        }
      }
    }

    /*
      setfhv( "user0", &max_time,  &nerr, SAC_STRING_LENGTH);
      setfhv( "user1", &max_value, &nerr, SAC_STRING_LENGTH);
    */

    /*   Write out the correlation function   */
    kname = strdup("correlatec_out1.sac");
    wsac0(kname, xarray, out, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
      fprintf(stderr, "Error writing out file: %s\n", kname);
      exit(-1);
    }

    free(out);

    return 0;
}

Convolution

Convolution of two time series

Fortran

      program convolvef

!       Reads in a (short) pulse that is convolved with (longer)
!         waveform.  Lengths are n_p and n_w.  Output (conv)
!         has length n_w + n_p + 1.  delta must be same for both.
!       Easiily expanded to read in multiple long time series.
!       If discconv is "y", the output is not premultilied by delta
!         and the begin time for the pulse is treated as zero
!
!       gfortran -o convolvef convolvef.f -L/usr/local/sac/lib -lsacio

      implicit none

      integer i,j
!     Define the Maximum length of waveform
      integer MAX
      parameter (MAX=10000)

      real waveform, pulse, ytmp, conv
      dimension waveform(MAX),pulse(MAX),ytmp(MAX),conv(MAX)
      character*16 kevnm

!     Declare Variables used in the rsac1() calls
      real b_w, delta, b_p, delta_in, b_p_in
      integer n_w, n_p, nmarg, iargc
      character*80 wf_name, p_name, c_name, kname
      character*1 discconv
      integer nerr

      nmarg = iargc()
      if (nmarg .eq. 0) then
      write(*,*)
     1   'Usage: convolvef wf_name p_name c_name discconv'
      write(*,*) 'where the first three arguments are'
      write(*,*) '  filenames for pulse, waveform, and output conv'
        write(*,*) 'discconv is y if uses a discrete convolution'
        write(*,*) '  otherwise time-series convolution'
        stop
      end if
      call getarg(1,p_name)
      call getarg(2,wf_name)
      call getarg(3,c_name)
      call getarg(4,discconv)

!   Read in pulse time series
      call rsac1(p_name, pulse, n_p, b_p, delta, MAX, nerr)

      if(nerr .NE. 0) then
         write(*,*)'Error reading in file: ',p_name
         call exit(-1)
       endif

!  Test if want to do a discrete convolution

       delta_in = delta
       b_p_in = b_p
       if (discconv .eq. 'y') then
         delta_in = 1.0
           b_p_in = 0.0
       endif

! If wanted to do more than one waveform, do loop starts here

!    Read in waveform time series

      call rsac1(wf_name, waveform, n_w, b_w, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',wf_name
          call exit(-1)
       endif

!     Do the convolution

      call td_conv(waveform,n_w,pulse,n_p,conv,delta_in,b_p_in)

      call setnhv('npts',n_w+n_p-1,nerr)
      kevnm = 'Convolution'
      call setkhv ('kevnm', kevnm, nerr)
!     Write the SAC file
      call wsac0(c_name, ytmp, conv, nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',c_name,nerr
          call exit(-1)
      endif

! If doing more than one waveform, do loop stops here

      call exit(0)

      end program convolvef

c+
      subroutine td_conv(waveform,n_w,pulse,n_p,conv,delta,b_p)
C
C       waveform of length n_w is the time series against which pulse
C         of length n_p is convolved.  Output: conv of length n_w+n_w+1.
C         waveform and pulse are unchanged.
C       The convoluton is done as an inner product in the time
C         domain.
C       Stops if n_w < n_p.
C
C       Arthur Snoke 2015
C-
        implicit none
        real*4 waveform(*), pulse(*), conv(*)
        real*4 delta,b_p,temp
        integer n_p, n_w, j_1, i, j
C
        if (n_p .ge. n_w) then
                write(*,*) 'Need more points in waveform than pulse'
                write(*,*) 'n_w and n_p:',n_w, n_p
                stop
        end if
C
        j_1 = -nint(b_p/delta)+1
        do i=1,n_w+n_p-1
          temp = 0.0
          do j=1,n_w
            if (i.ge.(j-j_1) .and. n_p.ge.(i-j+j_1))
     1          temp = temp + waveform(j)*pulse(i-j+j_1)
          end do
          conv(i) = delta*temp
        end do
        return
        end

C

/**
 * @file   td_conv.c
 *
 * @brief  Compute convolution of a long series (waveform) with pulse
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/**
 * Compute the Cross-Convolution of a pulse with a time series.
   If discconv is y, does a discrete convolution.
 *
 * @param waveform
 *    Array containing input time series
 * @param n_w
 *    Number of samples in waveform
 * @param pulse
 *    Array containing the input pulse time series
 * @param n_p
 *    Number of samples in pulse
 * @param conv
 *    Array containing the output time series conv
 * @param n_conv
 *    Number of samples in conv
 * @param delta
 *    Time interval for waveform, pulse, conv
 * @param b_p
 *    begin time of pulse time series
 *
 * @return Nothing
 *
 * \author   Arthur Snoke
 *           VT
 *
 * \date: 151105  Created
 *
 */
static void td_conv(
             float     *waveform,
             int        n_w,
             float     *pulse,
             int        n_p,
             float     *conv,
             float      delta,
             float      b_p)
{
  int i, j, j_1;
  float temp;

  if (n_p >= n_w) {
    fprintf(stderr, "Error: waveform and pulse lengths %d %d\n", n_w, n_p);
    exit(-1);
  }

  j_1 = -lrint(b_p/delta);

  for(i=0; i < n_w+n_p-1; i++){
    temp = 0.0;
    for(j=0; j < n_w; j++) {
        if (i >= (j-j_1) && n_p > (i-j+j_1)) {
            temp = temp + waveform[j]*pulse[i-j+j_1];
        }
    }
    conv[i] = delta*temp;
  }
  return;
} /* end of function td_conv*/


/*  convolvec.c
       Reads in a (short) pulse that is convolved with a (longer)
         waveform.  Lengths are n_p and n_w.  Output (conv)
         has length n_w + n_p + 1.  delta must be same for both.
       Easiily expanded to read in multiple long time series.
 gcc -o convolvec convolvec.c -I/usr/local/sac/include  -L/usr/local/sac/lib  -lsacio
 */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <sacio.h>

#define MAX        10000
#define ERROR_MAX  256


int main(int argc, char *argv[]) {

    /* Local variables */
    int i, j;
    int n_w, n_p, n_conv, nerr, max;

    float b_w, b_p, delta, b_p_in, delta_in;
    char *wf_name, *p_name, *c_name, *discconv;

    float waveform[MAX], pulse[MAX], conv[MAX], dummy[MAX];

    char error[ERROR_MAX];

    if (argc == 1)
        {
          fprintf(stderr, "usage: convolvef o_name wf_name c_name discconv\n");
          fprintf(stderr, "where the first three arguments are\n");
          fprintf(stderr, "  filenames for pulse, waveform, and output conv\n");
          fprintf(stderr, "discconv is y if uses a discrete convolution\n");
          fprintf(stderr, "  otherwise time-series convolution\n");
          exit(-1);
        }
    p_name = argv[1];
    wf_name = argv[2];
    c_name = argv[3];
    discconv = argv[4];

    max = MAX;

    for(i = 0; i < MAX; i++) {
      pulse[i] = 0.0;
    }

    /* Read in the pulse time series  */

    rsac1(p_name, pulse, &n_p, &b_p, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
        fprintf(stderr, "Error reading in file(%d): %s\n", nerr, p_name);
        exit(-1);
    }

  /*  Test if want to do a discrete convolution*/

       delta_in = delta;
       b_p_in = b_p;
       if (discconv[0] == 'y') {
         delta_in = 1.0;
           b_p_in = 0.0;
       }

    /* If wanted to do more than one waveform, do loop starts here */

    for(i = 0; i < MAX; i++) {
      waveform[i] = 0.0;
      conv[i] = 0.0;
      dummy[i] = 0.0;
    }

    /* Read in the waveform time series */

    rsac1(wf_name, waveform, &n_w, &b_w, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
        fprintf(stderr, "Error reading in file: %s\n", wf_name);
        exit(-1);
    }


  td_conv(waveform,n_w,pulse,n_p,conv,delta_in,b_p_in);

  n_conv = n_w+n_p-1;
  setnhv ( "npts",   &n_conv,    &nerr, SAC_STRING_LENGTH);
  setkhv ( "kevnm",  "Convolution", &nerr, SAC_STRING_LENGTH, SAC_STRING_LENGTH);

  /* Write output SAC file */

    wsac0(c_name, dummy, conv, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
        fprintf(stderr, "Error writing out file: %s\n", c_name);
        exit(-1);
    }

    /* If doing more than one waveform, do loop stops here */

    return 0;

} /* end of program convolvec*/

Sample Runs for Convolution

The SAC convolve command uses two digital-convolution features that may be undesirable: no scaling of the output by the digitizing interval, and no check on the start time for the pulse. To get the default SAC output, the last argument in the programs is y. After compiling the programs, do the following:

SAC> fg triangle npts 8 delta 1 begin -4
SAC> write triangle_n8_d1.sac
SAC> fg impulse npts 12 delta 1 begin 0
SAC> write impulse_n12_d1.sac

Then:

convolvef
    Usage: convolvef wf_name p_name c_name discconv
    where the first three arguments are filenames for pulse, waveform, and output conv
    discconv is y if uses a discrete convolution otherwise time-series convolution
convolvef triangle_n8_d1.sac impulse_n12_d1.sac conv_f.sac n
or similar entries for convolvec.