# FFT Spectrum Comparison

This is my latest coding project: a way to compare the spectra of two waveforms.

The plan was to be able to plot the PSD and/or the FFT from two waveforms on the one graph for clear comparison. This means both sample waveforms need to be the same length. I’ve also constrained the filter frequencies (if used) to be the same on both streams.

Most other things can be changed for the comparison. i.e. different stations, different channels, different outputs (displacement, velocity, acceleration, or Infrasound Pressure), raw counts or filtered units. Not all combinations of comparison units are necessarily meaningful, but may serve, say, an education purpose. For example, compare the same sample time and station and channel in both velocity and acceleration to show the bias towards higher frequencies when using acceleration output.

My main reason for doing this is to explore what frequencies are produced by some events, particularly with a Boom: lightning, wind, rain, etc.

Apart from displaying the PSD and FFT of two sample waveforms on the same graph, why not plot the differences? So I had a go at that as well for the FFT differences, but it turns out to not be as simple as it would appear at first thought. Somewhere in the FFT process complex numbers appear, so simple arithmetic subtraction of one FFT plot from another produces a warning that the complex component is being ignored.

I’ve come up with three ways to calculate the FFT difference (I suspect there are more - at least 1 for sure). At this stage, I don’t know which if any method is correct/valid but interestingly they are all consistent regardless of whether sample 1 is subtracted from sample 2 or vice versa. If anyone with some experience or understanding of FFT signal procressing can offer any constructive critcism on which method is valid, this would be greatly appreciated. For now, I’m plotting all three methods that I’ve coded.

Another query someone may be able to help with is how to subtract a PSD from another PSD. It’s not immediately obvious if this is possible the way the PSD’s are handled in Obspy.

Anyway, FWIW here is a first example of the output from the code. This is comparision of the arrival of the P wave from the M5.7 Kermedec Islands Quake with background noise just before the arrival:

As you can see, after comparing the raw counts, it is possible to apply a filter to both samples and hence “zoom in” to whatever areas of the frequency spectrum is of interest.

Here is the code so far. Probably not the tidiest, but happy for anyone to use and play with it and modify as they see fit.

``````
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 16 10:52:26 2023

@author: al72
"""

# Spectrum comparison for the Raspberry Shake and Boom

from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime
import matplotlib.pyplot as plt
import numpy as np

rs = Client('https://data.raspberryshake.org/')

def FFTsub(f1, f2):
n = len(f1)
diff = []
for i in range(0,n):
if f1[i] >= f2[i]:
diff.append(abs(f1[i]) - abs(f2[i]))
else:
diff.append(abs(f2[i]) - abs(f1[i]))
return diff

# Channel 1 data
startTime1 = UTCDateTime(2023, 11, 13, 14, 42, 54) # (YYYY, m, d, H, M, S) **** Enter data****

stn1 = 'R21C0'      # station name
nw1 = 'AM'          # network name
#ch1 = ['EHZ', 'DISP', 'm.']
ch1 = ['EHZ', 'VEL', 'm/s.']
#ch1 = ['EHZ', 'ACC', 'm/s².']
#ch1 = ['HDF', 'DEF', 'Pa.']
inv1 = rs.get_stations(network=nw1, station=stn1, level='RESP')  # get the instrument response
notes1 = 'Background Noise'

# Channel 2 data
startTime2 = UTCDateTime(2023, 11, 13, 14, 45, 45) # (YYYY, m, d, H, M, S) **** Enter data****

stn2 = 'R21C0'      # station name
nw2 = 'AM'          # network name
#ch2 = ['EHZ', 'DISP', 'm.']
ch2 = ['EHZ', 'VEL', 'm/s.']
#ch2 = ['EHZ', 'ACC', 'm/s².']
#ch2 = ['HDF', 'DEF', 'Pa.']
inv2 = rs.get_stations(network=nw2, station=stn2, level='RESP')  # get the instrument response
notes2 = 'M5.7 Quake Kermadec Islands P arrival'

#Setup the data plot
duration = 90               # duration of plots in seconds **** Enter data****
end1 = startTime1 + duration
end2 = startTime2 + duration
filtered = True
logspect = False
filt = [0.69, 0.7, 2, 2.1]

# calculate spectrum plot limits if filtered
if filtered:
if filt[2] >=4:
fxt = filt[2]+1
else:
fxt = filt[2]*1.25
if filt[1] >= 2:
fxb = filt[1]-1
else:
fxb = filt[1]*0.75
else:
fxt = 50
fxb = 0.05

# get and process the waveforms
trace1 = rs.get_waveforms(nw1, stn1, '00', ch1[0], startTime1, end1)
trace1.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace1.detrend(type='demean')                       #demean the data
raw1 = trace1.copy()
trace1.remove_response(inventory=inv1,pre_filt=filt,output=ch1[1],water_level=60, plot=False)

trace2 = rs.get_waveforms(nw2, stn2, '00', ch2[0], startTime2, end2)
trace2.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace2.detrend(type='demean')                       #demean the data
raw2 = trace2.copy()
trace2.remove_response(inventory=inv2,pre_filt=filt,output=ch2[1],water_level=60, plot=False)

# set up plot
fig = plt.figure(figsize=(20,14), dpi=150)       # set to page size in inches
ax1 = fig.add_subplot(5, 2, 1)      # waveform 1 plot
ax2 = fig.add_subplot(5, 2, 2)      # waveform 2 plot
ax3 = fig.add_subplot(5, 2, 3)      # waveform 1 spectrogram
ax4 = fig.add_subplot(5, 2, 4)      # waveform 2 spectrogram
ax5 = fig.add_subplot(5, 2, (5,6))      # PSD plots
ax6 = fig.add_subplot(5, 2, (7,8))      # FFT plots
ax7 = fig.add_subplot(5, 2, (9,10))     # FFT difference plots

#plot traces
if filtered:
ax1.plot(trace1[0].times(reftime=startTime1), trace1[0].data, lw=1, color='b')      # displacement waveform
ax2.plot(trace2[0].times(reftime=startTime2), trace2[0].data, lw=1, color='g')
ax1.set_ylabel(ch1[1]+', '+ch1[2],size='small')
ax2.set_ylabel(ch2[1]+', '+ch2[2],size='small')
else:
ax1.plot(raw1[0].times(reftime=startTime1), raw1[0].data, lw=1, color='b')      # displacement waveform
ax2.plot(raw2[0].times(reftime=startTime2), raw2[0].data, lw=1, color='g')
ax1.set_ylabel(ch1[1]+', counts.',size='small')
ax2.set_ylabel(ch2[1]+', counts.',size='small')
ax1.margins(x=0)
ax2.margins(x=0)

ax3.specgram(x=raw1[0], NFFT=256, noverlap=128, Fs=100, cmap='plasma')         # velocity spectrogram
ax3.set_ylabel(ch1[1]+' Spectrogram',size='small')
if logspect:
ax3.set_yscale('log')
ax3.set_ylim(0.05, 50)
ax4.specgram(x=raw2[0], NFFT=256, noverlap=128, Fs=100, cmap='plasma')         # velocity spectrogram
ax4.set_ylabel(ch2[1]+' Spectrogram',size='small')
if logspect:
ax4.set_yscale('log')
ax4.set_ylim(0.05, 50)

#plot PSD
#calculate NFFT for PSD
if duration >= 82:
nfft = 8192
else:
nfft = duration*100
if filtered:
ax5.psd(x=trace1[0], NFFT=nfft, noverlap=0, Fs=100, color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax5.psd(x=trace2[0], NFFT=nfft, noverlap=0, Fs=100, color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
#plot filter limits on PSD
ax5.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='r')
ax5.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='r')
else:
ax5.psd(x=raw1[0], NFFT=nfft, noverlap=0, Fs=100, color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax5.psd(x=raw2[0], NFFT=nfft, noverlap=0, Fs=100, color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
ax5.set_xlim(fxb, fxt)
#ax5.set_ylim(-290, -130)
ax5.legend(frameon=False, fontsize='x-small')
ax5.set_xscale('log')               #use logarithmic scale on PSD
#ax5.set_yscale('linear')
#ax5.set_yscale('log')
ax5.set_ylabel("PSD, dB",size='small')
ax5.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=-9)

# fourier analysis plot
if filtered:
fft1 = np.fft.rfft(trace1[0].data)
fft2 = np.fft.rfft(trace2[0].data)
xfft = np.fft.rfftfreq(trace1[0].data.size, d = 1/100)
# plot filter limits on FFT
ax6.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='r')
ax6.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='r')
else:
fft1 = np.fft.rfft(raw1[0].data)
fft2 = np.fft.rfft(raw2[0].data)
xfft = np.fft.rfftfreq(raw1[0].data.size, d = 1/100)
ax6.set_xlim(fxb, fxt)
ax6.plot(xfft, abs(fft1), color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax6.plot(xfft, abs(fft2), color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
ax6.legend(frameon=False, fontsize='x-small')
ax6.set_xscale('log')               #use logarithmic scale on FFT
#ax6.set_yscale('linear')
#ax6.set_yscale('log')
ax6.set_ylabel("FFT Spectrum",size='small')
ax6.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=0)

# fourier difference plot
dfft = FFTsub(fft1, fft2)
dfft1 = fft1 - fft2
dfft2 = abs(fft1) - abs(fft2)
ax7.plot(xfft, dfft, color='r', lw=1, label='FFT Magnitude Difference')
ax7.plot(xfft, abs(dfft1), color='purple', lw=1, label='FFT Abs Vector Difference')
ax7.plot(xfft, abs(dfft2), color='orange', lw=1, label='FFT Arithmetic Abs Difference')
ax7.legend(frameon=False, fontsize='x-small')
ax7.margins(x=0)
if filtered:
ax7.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='b')
ax7.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='b')
ax7.set_xlim(fxb, fxt)
ax7.set_xscale('log')               #use logarithmic scale on FFT Difference
ax7.set_ylabel("FFT Spectrum Difference",size='small')
ax7.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=0)

fig.suptitle("FFT Signal Comparison", weight='black', color='b', size='x-large')      #Title of the figure
fig.text(0.13, 0.94, nw1+'.'+stn1+'.00.'+ch1[0]+'.')
fig.text(0.13, 0.93, 'Start time = '+startTime1.strftime(' %d/%m/%Y %H:%M:%S UTC'))
fig.text(0.13, 0.92, 'Duration = '+str(duration)+' seconds.')
if filtered:
fig.text(0.13, 0.91, 'Filter = '+str(filt[1])+' to '+str(filt[2])+' Hz.')
else:
fig.text(0.13, 0.91, 'Raw Counts, Unfiltered.')
fig.text(0.13, 0.9, notes1)

fig.text(0.55, 0.94, nw2+'.'+stn2+'.00.'+ch2[0]+'.')
fig.text(0.55, 0.93, 'Start time = '+startTime2.strftime(' %d/%m/%Y %H:%M:%S UTC'))
fig.text(0.55, 0.92, 'Duration = '+str(duration)+' seconds.')
if filtered:
fig.text(0.55, 0.91, 'Filter = '+str(filt[1])+' to '+str(filt[2])+' Hz.')
else:
fig.text(0.55, 0.91, 'Raw Counts, Unfiltered.')
fig.text(0.55, 0.9, notes2)
``````
3 Likes

Fascinating!

I made a note to properly re-read this during my time off from work, so I can try with data recorded from my Shake. I am curious to see if differences between various cars/trucks that pass on a close-by road can be spotted with this method, alongside differences between events captured by the same instrument.

Thank you for sharing the code with us all in the community sheeny!

OK I solved the FFT difference dilemma.

dfft2 in the original code is the correct method. In the re-write I’ve renamed this to dfft.

I put a debug trap in after the dfft’s were calculated and printed just 10 instances of fft1, fft2, dfft, dfft1 and dfft2 (to avoid being swamped in data). I had to change the starting point for the trapped data by trial and error till I found different figures in dfft, dfft1 and dfft2. Then I copied this data into a spreadsheet, calculated the correct result from fft1 and fft2 and discarded the others.

Here’s the updated output and code:

``````
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 16 10:52:26 2023

@author: al72
"""

# Spectrum comparison for the Raspberry Shake and Boom

from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime
import matplotlib.pyplot as plt
import numpy as np

rs = Client('https://data.raspberryshake.org/')

# Channel 1 data
startTime1 = UTCDateTime(2023, 11, 13, 14, 42, 54) # (YYYY, m, d, H, M, S) **** Enter data****

stn1 = 'R21C0'      # station name
nw1 = 'AM'          # network name
#ch1 = ['EHZ', 'DISP', 'm.']
ch1 = ['EHZ', 'VEL', 'm/s.']
#ch1 = ['EHZ', 'ACC', 'm/s².']
#ch1 = ['HDF', 'DEF', 'Pa.']
inv1 = rs.get_stations(network=nw1, station=stn1, level='RESP')  # get the instrument response
notes1 = 'Background Noise'

# Channel 2 data
startTime2 = UTCDateTime(2023, 11, 13, 14, 45, 45) # (YYYY, m, d, H, M, S) **** Enter data****

stn2 = 'R21C0'      # station name
nw2 = 'AM'          # network name
#ch2 = ['EHZ', 'DISP', 'm.']
ch2 = ['EHZ', 'VEL', 'm/s.']
#ch2 = ['EHZ', 'ACC', 'm/s².']
#ch2 = ['HDF', 'DEF', 'Pa.']
inv2 = rs.get_stations(network=nw2, station=stn2, level='RESP')  # get the instrument response
notes2 = 'M5.7 Quake Kermadec Islands P arrival'

#Setup the data plot
duration = 90               # duration of plots in seconds **** Enter data****
end1 = startTime1 + duration
end2 = startTime2 + duration
filtered = False        # True to apply a frequency bandpass filter
logspect = False        # True to display logarithmic spectrogram Y axes
filt = [0.69, 0.7, 2, 2.1]      # Enter bandpass filter corners

# calculate spectrum plot limits if filtered
if filtered:
if filt[2] >=4:
fxt = filt[2]+1
else:
fxt = filt[2]*1.25
if filt[1] >= 2:
fxb = filt[1]-1
else:
fxb = filt[1]*0.75
else:
fxt = 50
fxb = 0.05

# get and process the waveforms
trace1 = rs.get_waveforms(nw1, stn1, '00', ch1[0], startTime1, end1)
trace1.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace1.detrend(type='demean')                       #demean the data
raw1 = trace1.copy()
trace1.remove_response(inventory=inv1,pre_filt=filt,output=ch1[1],water_level=60, plot=False)

trace2 = rs.get_waveforms(nw2, stn2, '00', ch2[0], startTime2, end2)
trace2.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace2.detrend(type='demean')                       #demean the data
raw2 = trace2.copy()
trace2.remove_response(inventory=inv2,pre_filt=filt,output=ch2[1],water_level=60, plot=False)

# set up plot
fig = plt.figure(figsize=(20,14), dpi=150)       # set to page size in inches
ax1 = fig.add_subplot(5, 2, 1)      # waveform 1 plot
ax2 = fig.add_subplot(5, 2, 2)      # waveform 2 plot
ax3 = fig.add_subplot(5, 2, 3)      # waveform 1 spectrogram
ax4 = fig.add_subplot(5, 2, 4)      # waveform 2 spectrogram
ax5 = fig.add_subplot(5, 2, (5,6))      # PSD plots
ax6 = fig.add_subplot(5, 2, (7,8))      # FFT plots
ax7 = fig.add_subplot(5, 2, (9,10))     # FFT difference plot

#plot traces
if filtered:
ax1.plot(trace1[0].times(reftime=startTime1), trace1[0].data, lw=1, color='b')      # displacement waveform
ax2.plot(trace2[0].times(reftime=startTime2), trace2[0].data, lw=1, color='g')
ax1.set_ylabel(ch1[1]+', '+ch1[2],size='small')
ax2.set_ylabel(ch2[1]+', '+ch2[2],size='small')
else:
ax1.plot(raw1[0].times(reftime=startTime1), raw1[0].data, lw=1, color='b')      # displacement waveform
ax2.plot(raw2[0].times(reftime=startTime2), raw2[0].data, lw=1, color='g')
ax1.set_ylabel(ch1[1]+', counts.',size='small')
ax2.set_ylabel(ch2[1]+', counts.',size='small')
ax1.margins(x=0)
ax2.margins(x=0)

ax3.specgram(x=raw1[0], NFFT=256, noverlap=128, Fs=100, cmap='plasma')         # velocity spectrogram
ax3.set_ylabel(ch1[1]+' Spectrogram',size='small')
if logspect:
ax3.set_yscale('log')
ax3.set_ylim(0.05, 50)
ax4.specgram(x=raw2[0], NFFT=256, noverlap=128, Fs=100, cmap='plasma')         # velocity spectrogram
ax4.set_ylabel(ch2[1]+' Spectrogram',size='small')
if logspect:
ax4.set_yscale('log')
ax4.set_ylim(0.05, 50)

#plot PSD
#calculate NFFT for PSD
if duration >= 82:
nfft = 8192
else:
nfft = duration*100
if filtered:
ax5.psd(x=trace1[0], NFFT=nfft, noverlap=0, Fs=100, color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax5.psd(x=trace2[0], NFFT=nfft, noverlap=0, Fs=100, color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
#plot filter limits on PSD
ax5.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='r')
ax5.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='r')
else:
ax5.psd(x=raw1[0], NFFT=nfft, noverlap=0, Fs=100, color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax5.psd(x=raw2[0], NFFT=nfft, noverlap=0, Fs=100, color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
ax5.set_xlim(fxb, fxt)
#ax5.set_ylim(-290, -130)
ax5.legend(frameon=False, fontsize='x-small')
ax5.set_xscale('log')               #use logarithmic scale on PSD
#ax5.set_yscale('linear')
#ax5.set_yscale('log')
ax5.set_ylabel("PSD, dB",size='small')
ax5.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=-9)

# fourier analysis plot
if filtered:
fft1 = np.fft.rfft(trace1[0].data)
fft2 = np.fft.rfft(trace2[0].data)
xfft = np.fft.rfftfreq(trace1[0].data.size, d = 1/100)
# plot filter limits on FFT
ax6.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='r')
ax6.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='r')
else:
fft1 = np.fft.rfft(raw1[0].data)
fft2 = np.fft.rfft(raw2[0].data)
xfft = np.fft.rfftfreq(raw1[0].data.size, d = 1/100)
ax6.set_xlim(fxb, fxt)
ax6.plot(xfft, abs(fft1), color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax6.plot(xfft, abs(fft2), color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
ax6.legend(frameon=False, fontsize='x-small')
ax6.set_xscale('log')               #use logarithmic scale on FFT
#ax6.set_yscale('linear')
#ax6.set_yscale('log')
ax6.set_ylabel("FFT Spectrum",size='small')
ax6.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=0)

# fourier difference plot
dfft = abs(fft1) - abs(fft2)
ax7.plot(xfft, abs(dfft), color='r', lw=1, label='FFT Difference')
ax7.legend(frameon=False, fontsize='x-small')
ax7.margins(x=0)
if filtered:
ax7.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='b')
ax7.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='b')
ax7.set_xlim(fxb, fxt)
ax7.set_xscale('log')               #use logarithmic scale on FFT Difference
ax7.set_ylabel("FFT Spectrum Difference",size='small')
ax7.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=0)

fig.suptitle("FFT Signal Comparison", weight='black', color='b', size='x-large')      #Title of the figure
fig.text(0.13, 0.94, nw1+'.'+stn1+'.00.'+ch1[0]+'.')
fig.text(0.13, 0.93, 'Start time = '+startTime1.strftime(' %d/%m/%Y %H:%M:%S UTC'))
fig.text(0.13, 0.92, 'Duration = '+str(duration)+' seconds.')
if filtered:
fig.text(0.13, 0.91, 'Filter = '+str(filt[1])+' to '+str(filt[2])+' Hz.')
else:
fig.text(0.13, 0.91, 'Raw Counts, Unfiltered.')
fig.text(0.13, 0.9, notes1)

fig.text(0.55, 0.94, nw2+'.'+stn2+'.00.'+ch2[0]+'.')
fig.text(0.55, 0.93, 'Start time = '+startTime2.strftime(' %d/%m/%Y %H:%M:%S UTC'))
fig.text(0.55, 0.92, 'Duration = '+str(duration)+' seconds.')
if filtered:
fig.text(0.55, 0.91, 'Filter = '+str(filt[1])+' to '+str(filt[2])+' Hz.')
else:
fig.text(0.55, 0.91, 'Raw Counts, Unfiltered.')
fig.text(0.55, 0.9, notes2)
``````
2 Likes

An example of a truck passing our place (going uphill).

Seismic minus background:

Infrasound minus background:

And seismic v infrasound signals:

lol! seismic v infrasound is not such a success to compare in linear FFT given the difference in signal scaling!

Seismic v infrasound with filter:

Clearly comparison across different sensors/units can be problematic.

Al.

2 Likes

Another couple of examples - Trucks passing uphill versus downhill:

The spectrogram confirms the direction of travel. The “starburst” pattern surrounds the corner at the bottom of our hill and is caused by the trucks slowing to take the corner and then accelerating again. If the truck is travelling uphill the starburst is on the left (before the strongest signal) and downhill it’s on the right (after the strongest signal).

Distinctly different shapes in the FFT signature (signal - minus background).

Al.

2 Likes

Your results are correct, FFT and PSD are similar for short duration records.
The difference is seen in long-time records in which PSD remove the transient content leaving only the permanent content, something that the FFT does not do.

3 Likes

``````
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 16 10:52:26 2023

@author: al72
"""

# Spectrum comparison for the Raspberry Shake and Boom

from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
import numpy as np

rs = Client('https://data.raspberryshake.org/')

# Channel 1 data
startTime1 = UTCDateTime(2023, 11, 22, 2, 20, 0) # (YYYY, m, d, H, M, S) **** Enter data****

stn1 = 'R21C0'      # station name
nw1 = 'AM'          # network name
#ch1 = ['*HZ', 'DISP', 'm.', 'Displacement']
#ch1 = ['*HZ', 'VEL', 'm/s.', 'Velocity']
#ch1 = ['*HZ', 'ACC', 'm/s².', 'Acceleration']
ch1 = ['HDF', 'DEF', 'Pa.', 'Pressure']
inv1 = rs.get_stations(network=nw1, station=stn1, level='RESP')  # get the instrument response
notes1 = '82 mm/h rain'

# Channel 2 data
startTime2 = UTCDateTime(2023, 11, 22, 4, 46, 0) # (YYYY, m, d, H, M, S) **** Enter data****

stn2 = 'R21C0'      # station name
nw2 = 'AM'          # network name
#ch2 = ['*HZ', 'DISP', 'm.', 'Displacement']
#ch2 = ['*HZ', 'VEL', 'm/s.', 'Velocity']
#ch2 = ['*HZ', 'ACC', 'm/s².', 'Acceleration']
ch2 = ['HDF', 'DEF', 'Pa.', 'Pressure']
inv2 = rs.get_stations(network=nw2, station=stn2, level='RESP')  # get the instrument response
notes2 = 'Background: 0 mm/h rain'

#Setup the data plot
duration = 120               # duration of plots in seconds **** Enter data****
end1 = startTime1 + duration
end2 = startTime2 + duration
filtered = True        # True to apply a frequency bandpass filter
logspect = False        # True to display logarithmic spectrogram Y axes

#Infrasund Filter bands
filt = [0.04, 0.05, 49, 50]   # For Booms with 20s mechanical filter fitted
#filt = [0.04, 0.05, 1, 1.1]
#filt = [0.09, 1, 49, 50]      # For all Booms...
#filt = [0.9, 1, 10, 10.1]
#filt = [9.9, 10, 20, 20.1]
#filt = [19.9, 20, 30, 30.1]
#filt = [29.9, 30, 40, 40.1]
#filt = [39.9, 40, 49, 50]

#seismic filter bands
#filt = [0.04, 0.05, 49, 50]
#filt = [0.69, 0.7, 2, 2.1]
#filt = [0.69, 0.7, 10, 10.1]
#filt = [0.69, 0.7, 20, 20.1]
#filt = [0.69, 0.7, 49, 50]

# calculate spectrum plot limits if filtered
if filtered:
if filt[2] >=4:
fxt = filt[2]+1
else:
fxt = filt[2]*1.25
if filt[1] >= 2:
fxb = filt[1]-1
else:
fxb = filt[1]*0.75
else:
fxt = 50
fxb = 0.05

# get and process the waveforms
trace1 = rs.get_waveforms(nw1, stn1, '00', ch1[0], startTime1, end1)
trace1.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace1.detrend(type='demean')                       #demean the data
raw1 = trace1.copy()
trace1.remove_response(inventory=inv1,pre_filt=filt,output=ch1[1],water_level=60, plot=False)

trace2 = rs.get_waveforms(nw2, stn2, '00', ch2[0], startTime2, end2)
trace2.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace2.detrend(type='demean')                       #demean the data
raw2 = trace2.copy()
trace2.remove_response(inventory=inv2,pre_filt=filt,output=ch2[1],water_level=60, plot=False)

# set up plot
fig = plt.figure(figsize=(20,14), dpi=150)       # set to page size in inches
ax1 = fig.add_subplot(5, 2, 1)      # waveform 1 plot
ax2 = fig.add_subplot(5, 2, 2)      # waveform 2 plot
ax3 = fig.add_subplot(5, 2, 3)      # waveform 1 spectrogram
ax4 = fig.add_subplot(5, 2, 4)      # waveform 2 spectrogram
ax5 = fig.add_subplot(5, 2, (5,6))      # PSD plots
ax6 = fig.add_subplot(5, 2, (7,8))      # FFT plots
ax7 = fig.add_subplot(5, 2, (9,10))     # FFT difference plot

#plot traces
if filtered:
ax1.plot(trace1[0].times(reftime=startTime1), trace1[0].data, lw=1, color='b')      # displacement waveform
ax2.plot(trace2[0].times(reftime=startTime2), trace2[0].data, lw=1, color='g')
ax1.set_ylabel(ch1[3]+', '+ch1[2],size='small')
ax2.set_ylabel(ch2[3]+', '+ch2[2],size='small')
else:
ax1.plot(raw1[0].times(reftime=startTime1), raw1[0].data, lw=1, color='b')      # displacement waveform
ax2.plot(raw2[0].times(reftime=startTime2), raw2[0].data, lw=1, color='g')
ax1.set_ylabel(ch1[3]+', counts.',size='small')
ax2.set_ylabel(ch2[3]+', counts.',size='small')
ax1.margins(x=0)
ax2.margins(x=0)

# document max amplitudes
if filtered:
max1 = trace1.max()
max2 = trace2.max()
units1 = ch1[2]
units2 = ch2[2]
else:
max1 = raw1.max()
max2 = raw2.max()
units1 = 'Counts'
units2 = 'Counts'
ax1l, ax1r = ax1.get_xlim()
ax2l, ax2r = ax2.get_xlim()
ax1b, ax1t = ax1.get_ylim()
ax2b, ax2t = ax2.get_ylim()
ax1.text(ax1r*0.02, ax1t*0.8, 'Max amplitude = '+f"{abs(max1[0]):0.3E}"+units1)
ax2.text(ax2r*0.02, ax2t*0.8, 'Max amplitude = '+f"{abs(max2[0]):0.3E}"+units2)

#plot spectrograms
ax3.specgram(x=raw1[0], NFFT=256, noverlap=128, Fs=100, cmap='plasma')         # velocity spectrogram
ax3.set_ylabel(ch1[3]+' Spectrogram',size='small')
ax4.specgram(x=raw2[0], NFFT=256, noverlap=128, Fs=100, cmap='plasma')         # velocity spectrogram
ax4.set_ylabel(ch2[3]+' Spectrogram',size='small')
if filtered:    # plot filter limits
ax3.axhline(y=filt[1], lw=1, color='w', linestyle='dotted')
ax3.axhline(y=filt[2], lw=1, color='w', linestyle='dotted')
ax4.axhline(y=filt[1], lw=1, color='w', linestyle='dotted')
ax4.axhline(y=filt[2], lw=1, color='w', linestyle='dotted')
if logspect:
ax3.set_yscale('log')
ax4.set_yscale('log')
ax4.set_ylim(0.05, 50)

#plot PSD
#calculate NFFT for PSD
if duration >= 82:
nfft = 8192
else:
nfft = duration*100
if filtered:
ax5.psd(x=trace1[0], NFFT=nfft, noverlap=0, Fs=100, color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax5.psd(x=trace2[0], NFFT=nfft, noverlap=0, Fs=100, color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
#plot filter limits on PSD
ax5.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='r')
ax5.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='r')
else:
ax5.psd(x=raw1[0], NFFT=nfft, noverlap=0, Fs=100, color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax5.psd(x=raw2[0], NFFT=nfft, noverlap=0, Fs=100, color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
ax5.set_xlim(fxb, fxt)
#ax5.set_ylim(-290, -130)
ax5.legend(frameon=False, fontsize='x-small')
#ax5.set_xscale('log')               #use logarithmic scale on PSD
#ax5.set_yscale('linear')
#ax5.set_yscale('log')
ax5.set_ylabel("PSD, dB",size='small')
ax5.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=-9)
ax5.xaxis.set_minor_locator(AutoMinorLocator(10))

# fourier analysis plot
if filtered:
fft1 = np.fft.rfft(trace1[0].data)
fft2 = np.fft.rfft(trace2[0].data)
xfft = np.fft.rfftfreq(trace1[0].data.size, d = 1/100)
# plot filter limits on FFT
ax6.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='r')
ax6.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='r')
else:
fft1 = np.fft.rfft(raw1[0].data)
fft2 = np.fft.rfft(raw2[0].data)
xfft = np.fft.rfftfreq(raw1[0].data.size, d = 1/100)
ax6.set_xlim(fxb, fxt)
ax6.plot(xfft, abs(fft1), color='b', lw=1, label=nw1+'.'+stn1+'.00.'+ch1[0]+' '+ch1[1])
ax6.plot(xfft, abs(fft2), color='g', lw=1, label=nw2+'.'+stn2+'.00.'+ch2[0]+' '+ch2[1])
#ax6.legend(frameon=False, fontsize='x-small')
#ax6.set_xscale('log')               #use logarithmic scale on FFT
#ax6.set_yscale('linear')
#ax6.set_yscale('log')
ax6.set_ylabel("FFT Spectrum",size='small')
ax6.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=0)
ax6.xaxis.set_minor_locator(AutoMinorLocator(10))

# fourier difference plot
dfft = abs(fft1) - abs(fft2)
ax7.plot(xfft, abs(dfft), color='r', lw=1, label='FFT Difference')
#ax7.legend(frameon=False, fontsize='x-small')
ax7.margins(x=0)
if filtered:
ax7.axvline(x=filt[1], linewidth=1, linestyle='dotted', color='b')
ax7.axvline(x=filt[2], linewidth=1, linestyle='dotted', color='b')
ax7.set_xlim(fxb, fxt)
#ax7.set_xscale('log')               #use logarithmic scale on FFT Difference
ax7.set_ylabel("FFT Spectrum Difference",size='small')
ax7.set_xlabel('F r e q u e n c y ,   H z', size='small', alpha=0.5, labelpad=0)
ax7.xaxis.set_minor_locator(AutoMinorLocator(10))

# Calculate Trend Lines
# Note high order polynomial fits may generate Rank Warnings for poor fit conditioning.
# Visually check the fit in these cases. If not acceptable reduce/adjust torder until
# fit is aceptable or Rank Warning disappears.
n = len(xfft)
trx = []
#tr1 = []
#tr2 = []
trd = []
torder = 15     #enter polynomial order for trend line

for i in range(0,n):
if xfft[i] >= filt[1] and xfft[i] <= filt[2]:   # only use data in the filter range!
trx.append(xfft[i])
#tr1.append(fft1[i])
#tr2.append(fft2[i])
trd.append(dfft[i])

#z1 = np.polyfit(trx, tr1, torder)
#z2 = np.polyfit(trx, tr2, torder)
zd = np.polyfit(trx, trd, torder)   # max polynomial trend line ... visually check fit especially if a rank warning is issued!
zav = np.polyfit(trx, trd, 0)       # order 0 = constant = average
zlin = np.polyfit(trx, trd, 1)      # order 1 = linear ... general trend.

#tl1 = np.poly1d(z1)
#tl2 = np.poly1d(z2)
tld = np.poly1d(zd)
tlav = np.poly1d(zav)
tllin = np.poly1d(zlin)

#ax6.plot(trx, abs(tl1(trx)), color = 'r', lw = 2, linestyle = 'dotted', label = 'FFT Blue Trend')
#ax6.plot(trx, abs(tl2(trx)), color = 'purple', lw = 2, linestyle = 'dotted', label = 'FFT Green Trend')
ax7.plot(trx, abs(tld(trx)), color = 'k', lw = 2, linestyle = 'dotted', label = 'FFT Difference Trend (order = '+str(torder)+')')
ax7.plot(trx, abs(tlav(trx)), color = 'b', lw = 2, linestyle = 'dotted', label = 'FFT Mean = '+f"{tlav.c[0]:0.3E}")
if tllin.c[1] < 0:
tlt = 'x - '
else:
tlt = 'x + '
ax7.plot(trx, abs(tllin(trx)), color = 'g', lw = 2, linestyle = 'dotted', label = 'FFT y = '+f"{tllin.c[0]:0.3E}"+tlt+f"{abs(tllin.c[1]):0.3E}")
ax6.legend(frameon=False, fontsize='x-small')
ax7.legend(frameon=False, fontsize='x-small')

fig.suptitle("FFT Signal Comparison", weight='black', color='b', size='x-large')      #Title of the figure
fig.text(0.13, 0.94, nw1+'.'+stn1+'.00.'+ch1[0]+'.')
fig.text(0.13, 0.93, 'Start time = '+startTime1.strftime(' %d/%m/%Y %H:%M:%S UTC'))
fig.text(0.13, 0.92, 'Duration = '+str(duration)+' seconds.')
if filtered:
fig.text(0.13, 0.91, 'Filter = '+str(filt[1])+' to '+str(filt[2])+' Hz.')
else:
fig.text(0.13, 0.91, 'Raw Counts, Unfiltered.')
fig.text(0.13, 0.9, notes1)

fig.text(0.55, 0.94, nw2+'.'+stn2+'.00.'+ch2[0]+'.')
fig.text(0.55, 0.93, 'Start time = '+startTime2.strftime(' %d/%m/%Y %H:%M:%S UTC'))
fig.text(0.55, 0.92, 'Duration = '+str(duration)+' seconds.')
if filtered:
fig.text(0.55, 0.91, 'Filter = '+str(filt[1])+' to '+str(filt[2])+' Hz.')
else:
fig.text(0.55, 0.91, 'Raw Counts, Unfiltered.')
fig.text(0.55, 0.9, notes2)
``````

The basics are the same, but I’ve tarted it up a little and added some functionality to compare trends and amplitudes, etc.

I intended it to do things like compare diurnal cultural noise variation, isolating “signatures” of different types of events, compare noise levels between stations, etc. From an educational point of view it could also be used to compare different outputs from the same shake and time (e.g. velocity v acceleration) to show the frequency bias of each output.

Al.

2 Likes

I’ve followed this topic during my days off, sheeny, and as usual, top quality!

Thank you for the polished code you’ve put at the end; I’m sure it will be handy as a ready-to-run program for whoever wants to take a look at this.

As I live close to a well-trafficked road, I wanted to see the differences between cars, trucks, buses and the like, as I wrote in my previous message. Alas, no time for that, but I’ll hopefully find it in the coming weeks.

1 Like