Seisplotjs Tutorial

Filtering, Deconvolution and FFT

See it live in tutorial5.html.

Maybe we would like to see the body waves more clearly and filter out the surface waves. And instead of the raw velocity data, lets apply the full instrument response to correct to displacement. First, lets switch from the LH? channels, recorded at one sample per second, to the HH? channels, recorded at 100 sps to give us a wider frequency band to play with. Since we will need the response information as well, we will need to change from using queryChannels() to queryResponses() in the stationQuery. This will return stationXML with the full response filled in for each channel.


            let stationQuery = new seisplotjs.fdsnstation.StationQuery()
              .networkCode('CO')
              .stationCode('HODGE')
              .locationCode('00')
              .channelCode('HH?')
              .timeWindow(queryTimeWindow);
            Promise.all( [ eventQuery.query(), stationQuery.queryResponses() ] )
          

Of course 100 sps means that we will have 100 times as many samples, so lets reduce the time window to just the region where the P arrival is, say 300 seconds instead of 1800.


          let timeWindow = new seisplotjs.util.StartEndDuration(startTime, null, 300);
          

Now we insert the filtering code after the seismograms have arrived but before we plot them. We will insert another then() call even though we could just add to the existing one. We create a butterworth filter using the sampling rate of the seismogram, with a passband of .5 to 10 Hz. Removing the mean is usually a good idea, then we apply the filter. Tapering is important before a deconvolution and then we transfer the instrument response for the channel.


          }).then( ( [ quakeList, networks, ttimes, seismogramDataList ] ) => {
            seismogramDataList.forEach(sdd => {
              let butterworth = seisplotjs.filter.createButterworth(
                                     2, // poles
                                     seisplotjs.filter.BAND_PASS,
                                     .5, // low corner
                                     10, // high corner
                                     1/sdd.seismogram.sampleRate // delta (period)
                            );
              let rmeanSeis = seisplotjs.filter.rMean(sdd.seismogram);
              let filteredSeis = seisplotjs.filter.applyFilter(butterworth, rmeanSeis);
              let taperSeis = seisplotjs.taper.taper(filteredSeis);
              let correctedSeis = seisplotjs.transfer.transfer(taperSeis,
                  sdd.channel.response, .001, .02, 250, 500);
              sdd.seismogram = correctedSeis;
            });
            return Promise.all( [ quakeList, networks, ttimes, seismogramDataList ] );
          

Then, just to make sure we don't correct the data twice, we disable the gain correction.


            seisConfig.doGain = false;
          

We can also plot the amplitude spectra for the three seismograms. We need to add an additional div to hold them.


            <div id="fftplot">
            </div>
          

And an additional style to size it.


            div#fftplot {
              height: 600px;
            }
          

Then we calculate the fft and plot it. Because the resulting plot is just an svg element created with d3, we can add a plot label by modifying it with normal d3 techniques.


          }).then( ( [ quakeList, networks, ttimes, seismogramDataList ] ) => {
            let div = seisplotjs.d3.select('div#fftplot');
            let fftList = seismogramDataList.map(sdd => seisplotjs.fft.fftForward(sdd.seismogram));
            let seisConfig = new seisplotjs.seismographconfig.SeismographConfig();
            let fftPlot = new seisplotjs.fftplot.FFTPlot('div#fftplot', seisConfig, fftList, true);
            fftPlot.draw();
            fftPlot.svg.append("g").classed("title", true)
              .attr("transform", "translate(600, 100)")
              .append("text").classed("title label", true)
              .selectAll("tspan")
              .data(seismogramDataList)
              .enter()
              .append("tspan")
              .text(sdd => " "+sdd.seismogram.channelCode+" ");

            return Promise.all( [ quakeList, networks, ttimes, seismogramDataList ] );
          

See it live in tutorial5.html.

Previous: Predicted phase arrival times

Next: Helicorder