Hello,

Im builing a program that analyses the ambient noise to infer the location of water. For this im using the H/V method proposed by Nakumura.

The technique originally proposed by Nogoshi and Igarashi (1971), and wide-spread by Nakamura (1989), consists in estimating the ratio between the Fourier amplitude spectra of the horizontal (H) to vertical (V) components of the ambient noise vibrations recorded at one single station.

The computation of the H/V ratio follows different steps (see Figure 1):

- record a 3-component ambient noise signal
- select of the most stationary time windows (e.g., using an anti-triggering algorithm) in order to avoid transient noise
- compute and smoothing of the Fourier amplitude spectra for each time windows
- Average the two horizontal component (using a quadractic mean)
- compute the H/V ratio for each window
- compute the average H/V ratio

Source: http://www.geopsy.org/documentation/geopsy/hv.html

I leave here the code that I been developing to test this method:

import numpy as np

from obspy.core import *

from obspy.arclink import Clientfrom obspy.signal import seisSim

from obspy.xseed import Parserimport matplotlib.pyplot as plt

st = read('http://examples.obspy.org/COP.BHE.DK.2009.050’)

st += read('http://examples.obspy.org/COP.BHN.DK.2009.050’)

st += read('http://examples.obspy.org/COP.BHZ.DK.2009.050’)

e_trace= st[0].copy()

n_trace = st[1].copy()

z_trace= st[2].copy()

client = Client(“webdc.eu”, 18001, user='test@obspy.org’)

for tr in [e_trace, n_trace, z_trace]:

start=tr.stats.starttime

end=start+3600

tr.trim(start,end) #only want 1hour of data for testing purposes

paz=client.getPAZ(tr.stats.network, tr.stats.station, tr.stats.location, tr.stats.channel,start)

df = tr.stats.sampling_rate

tr.data = seisSim(tr.data, df, paz_remove=paz) #desconvolution, converting the data from counts to m/s

#Get the amplitude values from the resulting traces

e_array=e_trace.data

n_array=n_trace.data

z_array=z_trace.data

#Removing the mean from the trace

e=e-np.mean(e)

n=n-np.mean(n)

z=z-np.mean(z)

#Apply the fourier transformation to the arrays

e_fft=np.fft.fft(e)

n_fft=np.fft.fft(n)

z_fft=np.fft.fft(z)

#Smoothing

e=smooth(e_fft)

n=smooth(n_fft)

V=smooth(z_fft)

ndouble=abs(np.power(n,2))

edouble=abs(np.power(e,2))

#Apply the quadratic mean

H=np.sqrt(np.add(edouble,ndouble))

#Devide the two components

HV=np.divide(H,V)

n=e.size

freq = np.fft.fftfreq(n, d=(1/df))

plt.figure()

plt.plot(freq,H)

plt.xlim(0,2)

plt.figure()

plot(freq,HV)

plt.xlim(0,2)

plt.show()

Where smooth is the following function (source):

def smooth(x,window_len=11,window=‘hanning’):

if x.ndim != 1:

raise ValueError, “smooth only accepts 1 dimension arrays.”

if x.size < window_len:

raise ValueError, “Input vector needs to be bigger than window size.”

if window_len<3:

return x

if not window in [‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’]:

raise ValueError, “Window is on of ‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’”

s=np.r_[2

x[0]-x[window_len-1::-1],x,2x[-1]-x[-1:-window_len:-1]]if window == ‘flat’: #moving average

w=np.ones(window_len,‘d’)

else:

w=eval(‘np.’+window+’(window_len)’)

y=np.convolve(w/w.sum(),s,mode=‘same’)

return y[window_len:-window_len+1]

My problem right now is finding the best way to smooth the waveform data. I researched quit a lot and tested many different methods and functions from the scipy. Im a newbie in obspy and python for that matter. Very few people in my college can help me.Does anybody have a good smoothing algorithm that can give me?

Thanks so much in advance.

Best regards,

Guilherme