Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add new different and less brutal whitening options #93

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion msnoise/api.py
Original file line number Diff line number Diff line change
@@ -226,7 +226,7 @@ def get_filters(session, all=False):


def update_filter(session, ref, low, mwcs_low, high, mwcs_high,
rms_threshold, mwcs_wlen, mwcs_step, used):
rms_threshold, mwcs_wlen, mwcs_step, used, whitening_method):
"""Updates or Insert a new Filter in the database.

.. seealso:: :class:`msnoise.msnoise_table_def.Filter`
@@ -255,6 +255,8 @@ def update_filter(session, ref, low, mwcs_low, high, mwcs_high,
:param mwcs_step: Step (in seconds) of the windowing procedure in MWCS
:type used: bool
:param used: Is the filter activated for the processing
:type whitening_method: str
:param whitening_method: Whitening method to use
"""
filter = session.query(Filter).filter(Filter.ref == ref).first()
if filter is None:
@@ -267,6 +269,7 @@ def update_filter(session, ref, low, mwcs_low, high, mwcs_high,
filter.mwcs_wlen = mwcs_wlen
filter.mwcs_step = mwcs_step
filter.used = used
filter.whitening_method = whitening_method
session.add(filter)
else:
filter.low = low
@@ -277,6 +280,7 @@ def update_filter(session, ref, low, mwcs_low, high, mwcs_high,
filter.mwcs_wlen = mwcs_wlen
filter.mwcs_step = mwcs_step
filter.used = used
filter.whitening_method = whitening_method
session.commit()
return

1 change: 1 addition & 0 deletions msnoise/default.py
Original file line number Diff line number Diff line change
@@ -27,6 +27,7 @@
default['overlap'] = ["Amount of overlap between data windows [0:1[ [0.]",'0.0']
default['windsorizing'] = ["Windsorizing at N time RMS , 0 disables windsorizing, -1 enables 1-bit normalization [3]",'3']
default['whitening'] = ["Whiten Traces before cross-correlation: All (except for autocorr), None, or only if components are different [A]/N/C",'A']
default['whitening_method'] = ["Whitening Method to use if Whitening is enabled: 'brutal', 'smoothing', 'smooth_whitening' or 'none'"]

default['stack_method'] = ["Stack Method: Linear Mean or Phase Weighted Stack: [linear]/pws ",'linear']
default['pws_timegate'] = ["If stack_method='pws', width of the smoothing in seconds : 10.0 ",'10.0']
116 changes: 80 additions & 36 deletions msnoise/move2obspy.py
Original file line number Diff line number Diff line change
@@ -6,6 +6,8 @@
import scipy.optimize
import scipy.signal
from obspy.signal.invsim import cosine_taper
from obspy.signal.konnoohmachismoothing \
import konno_ohmachi_smoothing_window as ko_window
from scipy.fftpack.helper import next_fast_len

try:
@@ -71,28 +73,32 @@ def myCorr(data, maxlag, plot=False, nfft=None):
return corr


def whiten(data, Nfft, delta, freqmin, freqmax, plot=False):
"""This function takes 1-dimensional *data* timeseries array,
goes to frequency domain using fft, whitens the amplitude of the spectrum
in frequency domain between *freqmin* and *freqmax*
and returns the whitened fft.
def whiten(data, Nfft, delta, freqmin, freqmax, plot=False, wtype='brutal'):
"""This function whitens the amplitude spectrum of the input trace. It is
modified from MSNoise to include different types of whitening."""

:type data: :class:`numpy.ndarray`
:param data: Contains the 1D time series to whiten
:type Nfft: int
:param Nfft: The number of points to compute the FFT
:type delta: float
:param delta: The sampling frequency of the `data`
:type freqmin: float
:param freqmin: The lower frequency bound
:type freqmax: float
:param freqmax: The upper frequency bound
:type plot: bool
:param plot: Whether to show a raw plot of the action (default: False)
def smoothen(spectrum, bandwidth, low, high):
"""Smoothens the given spectrum between low and high using a
Konno Ohmachi window. Returns an array with the smoothened part only"""

:rtype: :class:`numpy.ndarray`
:returns: The FFT of the input trace, whitened between the frequency bounds
"""
# half-width of the sliding window
wl = 90
# create the filter
ko_filter = ko_window(freqVec[:2*wl+1], freqVec[wl],
bandwidth, normalize=False)

# add zeros to the low end of the spectrum if filter doesn't fit
if low < wl:
offset = wl - low
spectrum = np.append(np.zeros(offset), spectrum)
low += offset
high += offset

return np.array([np.dot(spectrum[x - wl:x + wl + 1], ko_filter) \
for x in range(low, high)])

if wtype not in ['brutal', 'smoothing', 'smooth_whitening', 'none']:
raise ValueError('Unknown whitening type %s' % wtype)

if plot:
plt.subplot(411)
@@ -109,8 +115,8 @@ def whiten(data, Nfft, delta, freqmin, freqmax, plot=False):
if low <= 0:
low = 1

porte1 = J[0]
porte2 = J[-1]
p1 = J[0]
p2 = J[-1]
high = J[-1] + Napod
if high > Nfft / 2:
high = int(Nfft // 2)
@@ -124,33 +130,71 @@ def whiten(data, Nfft, delta, freqmin, freqmax, plot=False):
plt.xlim(0, max(axis))
plt.title('FFTRawSign')

# Left tapering:
if wtype == 'brutal':
# Left tapering: amplitude spectrum is the cosine taper itself
FFTRawSign[low:p1] = np.cos(
np.linspace(np.pi / 2., np.pi, p1 - low)) ** 2 * np.exp(
1j * np.angle(FFTRawSign[low:p1]))

# Pass band: amplitude spectrum is 1
FFTRawSign[p1:p2] = np.exp(1j * np.angle(FFTRawSign[p1:p2]))

# Right tapering: amplitude spectrum is the cosine taper itself
FFTRawSign[p2:high] = np.cos(
np.linspace(0., np.pi / 2., high - p2)) ** 2 * np.exp(
1j * np.angle(FFTRawSign[p2:high]))

else:
if wtype == 'smoothing':
# Pass band: a smoothened spectrum for the pass band
FFTRawSign[low:high] = (smoothen(np.abs(FFTRawSign), 10, low, high) *
np.exp(1j * np.angle(FFTRawSign[low:high])))

elif wtype == 'smooth_whitening':
# compute a smooth amplitude spectrum
smooth = smoothen(np.abs(FFTRawSign), 10, low, high)
waterlevel = 0.001 * np.max(smooth)
# add waterlevel
smooth[smooth < waterlevel] = waterlevel

# Pass band: the amplitudes are whitened using a smooth spectrum
FFTRawSign[low:high] = ((np.abs(FFTRawSign[low:high]) / smooth) *
np.exp(1j * np.angle(FFTRawSign[low:high])))

# Left tapering: amplitude spectrum is multiplied by cosine taper
FFTRawSign[low:p1] = (np.abs(FFTRawSign[low:p1]) *
np.cos(np.linspace(np.pi / 2., np.pi, p1 - low)) ** 2 *
np.exp(1j * np.angle(FFTRawSign[low:p1])))

# Right tapering: amplitude spectrum is multiplied by cosine taper
FFTRawSign[p2:high] = (np.abs(FFTRawSign[p2:high]) *
np.cos(np.linspace(0., np.pi / 2., high - p2)) ** 2 *
np.exp(1j * np.angle(FFTRawSign[p2:high])))

# set lower part of spectrum to zero
FFTRawSign[0:low] *= 0
FFTRawSign[low:porte1] = np.cos(
np.linspace(np.pi / 2., np.pi, porte1 - low)) ** 2 * np.exp(
1j * np.angle(FFTRawSign[low:porte1]))
# Pass band:
FFTRawSign[porte1:porte2] = np.exp(1j * np.angle(FFTRawSign[porte1:porte2]))
# Right tapering:
FFTRawSign[porte2:high] = np.cos(
np.linspace(0., np.pi / 2., high - porte2)) ** 2 * np.exp(
1j * np.angle(FFTRawSign[porte2:high]))

# set upper part of spectrum to zero
FFTRawSign[high:Nfft + 1] *= 0

# Hermitian symmetry (because the input is real)
FFTRawSign[-(Nfft // 2) + 1:] = FFTRawSign[1:(Nfft // 2)].conjugate()[::-1]

# Normalize the amplitude spectrum
FFTRawSign = (np.abs(FFTRawSign) / np.amax(FFTRawSign) *
np.exp(1j * np.angle(FFTRawSign)))

if plot:
plt.subplot(413)
axis = np.arange(len(FFTRawSign))
plt.axvline(low, c='g')
plt.axvline(porte1, c='g')
plt.axvline(porte2, c='r')
plt.axvline(p1, c='g')
plt.axvline(p2, c='r')
plt.axvline(high, c='r')

plt.axvline(Nfft - high, c='r')
plt.axvline(Nfft - porte2, c='r')
plt.axvline(Nfft - porte1, c='g')
plt.axvline(Nfft - p2, c='r')
plt.axvline(Nfft - p1, c='g')
plt.axvline(Nfft - low, c='g')

plt.plot(axis, np.abs(FFTRawSign))
13 changes: 11 additions & 2 deletions msnoise/msnoise_admin.py
Original file line number Diff line number Diff line change
@@ -182,18 +182,27 @@ def mwcs_step(form, field):
if field.data > form.data['mwcs_wlen']:
raise ValidationError("'mwcs_step' should be smaller or equal to"
" 'mwcs_wlen'")

def whitening_method(form, field):
whitening_methods = ('brutal', 'smoothing', 'smooth_whitening', 'none')
if field.data not in whitening_methods:
raise ValidationError(
"'whitening_method' should be one of %s" % whitening_methods)

form_args = dict(
mwcs_low=dict(validators=[mwcs_low]),
mwcs_high=dict(validators=[mwcs_high]),
high=dict(validators=[high]),
mwcs_step=dict(validators=[mwcs_step]),
whitening_method=dict(validators=[whitening_method]),
)

column_list = ('ref', 'low', 'mwcs_low', 'mwcs_high', 'high',
'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used')
'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used',
'whitening_method')
form_columns = ('low', 'mwcs_low', 'mwcs_high', 'high',
'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used')
'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used',
'whitening_method')

def __init__(self, session, **kwargs):
# You can pass name and other parameters if you want to
4 changes: 4 additions & 0 deletions msnoise/msnoise_table_def.py
Original file line number Diff line number Diff line change
@@ -31,6 +31,8 @@ class Filter(Base):
:param mwcs_step: Step (in seconds) of the windowing procedure in MWCS
:type used: bool
:param used: Is the filter activated for the processing
:type whitening_method: str
:param whitening_method: Whitening method to use
"""

__tablename__ = "filters"
@@ -44,6 +46,7 @@ class Filter(Base):
mwcs_wlen = Column(Float())
mwcs_step = Column(Float())
used = Column(Boolean(True))
whitening_method = Column(String(50))

def __init__(self, **kwargs):
""""""
@@ -55,6 +58,7 @@ def __init__(self, **kwargs):
# self.mwcs_wlen = mwcs_wlen
# self.mwcs_step = mwcs_step
# self.used = used
# self.whitening_method = whitening_method


class Job(Base):
8 changes: 6 additions & 2 deletions msnoise/s001configurator.py
Original file line number Diff line number Diff line change
@@ -56,6 +56,8 @@ class TFilter(HasTraits):
mwcs_wlen = Float
mwcs_step = Float
used = CBool
whitening_method = Str
# XXX do we have to add whitening_method to View??
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no, those views were for the TraitsUI, I think

traits_view = View('low', 'high', 'rms_threshold', 'used',
buttons=['OK', 'Cancel'])

@@ -105,6 +107,7 @@ class TStation(HasTraits):
StationColumn(name='rms_threshold'),
StationColumn(name='mwcs_wlen'),
StationColumn(name='mwcs_step'),
StationColumn(name='whitening_method'),
])

config_editor = TableEditor(
@@ -225,7 +228,8 @@ def main():
TFilter(
ref=f.ref, low=f.low, high=f.high, mwcs_low=f.mwcs_low,
mwcs_high=f.mwcs_high, rms_threshold=f.rms_threshold,
mwcs_wlen=f.mwcs_wlen, mwcs_step=f.mwcs_step, used=f.used))
mwcs_wlen=f.mwcs_wlen, mwcs_step=f.mwcs_step, used=f.used,
whitening_method=f.whitening_method))

configs = []
for name in default.keys():
@@ -263,7 +267,7 @@ def main():
for f in demo.company.filters:
update_filter(db, f.ref, f.low, f.mwcs_low, f.high,
f.mwcs_high, f.rms_threshold, f.mwcs_wlen,
f.mwcs_step, f.used)
f.mwcs_step, f.used, f.whitening_method)

db.close()
print("Done !")
8 changes: 7 additions & 1 deletion msnoise/s03compute_cc.py
Original file line number Diff line number Diff line change
@@ -84,6 +84,10 @@
example, computing ZZ components for very close by stations (much closer than
the wavelength sampled), leading to spatial autocorrelation issues.

The whitening method to use can be set in each filter's configuration as key
``whitening_method``. Currently supported options are "brutal" (MSNoise
default), "smoothing", "smooth_whitening" and "none".

When both traces are ready, the cross-correlation function is computed
(see :ref:`mycorr`). The function returned contains data for time lags
corresponding to ``maxlag`` in the acausal (negative lags) and causal
@@ -390,6 +394,7 @@ def main():
low = float(filterdb.low)
high = float(filterdb.high)
rms_threshold = filterdb.rms_threshold
whitening_method = filterdb.whitening_method

trames2hWb = np.zeros((2, int(nfft)), dtype=np.complex)
skip = False
@@ -399,7 +404,8 @@ def main():
#logging.debug("Whitening %s" % components)
trames2hWb[i] = whiten(tmp[i].data, nfft,
dt, low, high,
plot=False)
plot=False,
method=whitening_method)
else:
#logging.debug("Autocorr %s"%components)
tmp[i].filter("bandpass", freqmin=low,
7 changes: 4 additions & 3 deletions msnoise/scripts/msnoise.py
Original file line number Diff line number Diff line change
@@ -173,7 +173,7 @@ def d(path):
click.echo('')
click.echo('Filters:')
print('ID: [low:high] [mwcs_low:mwcs_high] mwcs_wlen mwcs_step'
' used')
' used whitening_method')
for f in get_filters(db, all=True):
data = (f.ref,
f.low,
@@ -182,8 +182,9 @@ def d(path):
f.mwcs_high,
f.mwcs_wlen,
f.mwcs_step,
['N', 'Y'][f.used])
print('%02i: [%.03f:%.03f] [%.03f:%.03f] %.03i %.03i %s' % data)
['N', 'Y'][f.used],
f.whitening_method)
print('%02i: [%.03f:%.03f] [%.03f:%.03f] %.03i %.03i %s %s' % data)

click.echo('')
click.echo('Stations:')
8 changes: 6 additions & 2 deletions msnoise/test/tests.py
Original file line number Diff line number Diff line change
@@ -63,6 +63,7 @@ def test_004_set_and_get_filters(self):
f.mwcs_wlen = 10
f.mwcs_step = 5
f.used = True
f.whitening_method = 'brutal'
filters.append(f)
f = Filter()
f.low = 0.1
@@ -73,16 +74,19 @@ def test_004_set_and_get_filters(self):
f.mwcs_wlen = 10
f.mwcs_step = 5
f.used = True
f.whitening_method = 'smooth_whitening'
filters.append(f)

for f in filters:
update_filter(db, f.ref, f.low, f.mwcs_low, f.high, f.mwcs_high,
f.rms_threshold, f.mwcs_wlen, f.mwcs_step, f.used)
f.rms_threshold, f.mwcs_wlen, f.mwcs_step, f.used,
f.whitening_method)

dbfilters = get_filters(db)
for i, filter in enumerate(dbfilters):
for param in ['low', 'mwcs_low', 'high', 'mwcs_high',
'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used']:
'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used',
'whitening_method']:
self.failUnlessEqual(eval("filter.%s" % param),
eval("filters[i].%s" % param))