Make some noise

I spent a long time looking at how to ‘create noise’ in order to make some synthetic seismograms, so I thought that I would put up my code in case anyone ends up in the same spot as me! I take several steps in order to model this:

  • Load in some typical seismic noise (I have taken mine from a quiet day near the Tunguruhua volcano in Equador), which has been detrended and demeaned.
  • Taking the Fast Fourier Transform (FFT) of this (this puts the data into the frequency domain).
  • Smooth the FFT data.
  • Multiply this by the FFT of white noise.
  • Take the Inverse Fast Fourier Transform (IFFT) of this (takes it back into the time domain).

The results of this are shown below, where the green is our white noise, the blue is our real seismic noise, and the pink is our synthetic seismic noise.

Createnoiseplot2_comp_of_white_and_T_and_created

Creating seismic noise

There are a few other intermediate steps to this code (such as looping through so that it is in segments), however it is quite a simple process! A few other libraries are loaded into this beforehand, such as Obspy and Numpy, however you will probably have loaded these in already if you are doing this.

Now go and make some noise!

— Roseanne


def noise_segmenting(poisson_times, st_event_2, st_t, noise_level, samp_rate, delta):
 """ Creates the noise array so that it is big enough to host all of the events.
 Creating the noise by multiplying white noise by the seismic noise, in the frequency domain.
 We then inverse FFT it and scale it to whatever SNR level is defined to output the full
 noise array.
poisson_times = array of times where we then put in the seismic events (boundary for the noise)
st_event_2 = size of events that we are putting in later (again, this is a boundary)
st_t = seismic noise array that you are basing your synthetic on
samp_rate, delta = trace properties of st_t
 """
 # end time for noise to cover all events
 noise_lim = (poisson_times[-1] + len(st_event_2)) *2 #gives some time after last event
# load in seismic noise to base the synthetic type on
 st_noise_start_t = UTCDateTime("2015-01-22T01:00:00")
 st_noise_end_t = UTCDateTime("2015-01-22T01:02:00")
 test_trace = st_t[0].slice(st_noise_start_t, st_noise_end_t)
 test_trace_length = int(len(test_trace) / test_trace.stats.sampling_rate)
# setting the boundary for how many loops etc
 minutes_long = (noise_lim)/st_event_2.stats.sampling_rate
 noise_loops = int(np.ceil(minutes_long/2.0)) #working out how many 2 minute loops we need
# zero array
 noise_array = np.zeros([noise_loops, len(test_trace)])
# loop for the amount of noise_loops needed (in segments)
 for j in range(noise_loops):
      # we average the seismic noise over twenty 2 minute demeaned samples
      tung_n_fft = np.zeros([20, int(np.ceil((len(test_trace)/2.0)))])
      for i in range(20):
         st_noise = st_t[0].slice(st_noise_start_t+(i*test_trace_length),st_noise_end_t+(i*test_trace_length))
         noise_detrended = st_noise.detrend()
         noise_demeaned = mlab.demean(noise_detrended)
         noise_averaging = Trace(noise_demeaned).normalize()
         tung_n_fft[i] = np.fft.rfft(noise_averaging.data)

      # work out the average fft
      ave = np.average(tung_n_fft, axis=0)
      # smooth the data
      aves = movingaverage(ave,20)
# create white noise
      whitenoise = np.random.normal(0, 1, len(noise_averaging))
      whitenoise_n = Trace(whitenoise).normalize()
# FFT the white noise
      wn_n_fft = np.fft.rfft(whitenoise_n.data)
# multiply the FFT of white noise and the FFT smoothed seismic noise
      newnoise_fft = wn_n_fft * aves
# IFFT the product
      newnoise = ifft(newnoise_fft, n = len(st_noise))

      noise_array[j] = np.real(newnoise)
# transform the noise into an Obspy trace
 full_noise_array = np.ravel(noise_array)
 full_noise_array_n = Trace(np.float32(full_noise_array)).normalize()
 full_noise_array_n_scaled = Trace(np.multiply(full_noise_array_n, noise_level))
 full_noise_array_n_scaled.stats.sampling_rate = samp_rate
 full_noise_array_n_scaled.stats.delta = delta

 return full_noise_array_n_scaled