# -*- coding: utf-8 -*-
# 2009-02-03 Moritz

import obspy
import numpy as np
from obspy.signal import cosTaper
import matplotlib.pyplot as plt
from matplotlib.mlab import detrend
plt.close('all')

# Load data, China earthquake
tr = obspy.read("data/GR.FUR..LHZ.D.2008.133")[0]
# Uncomment to zoom to a certain region
#tr.trim(tr.stats.starttime + 3000, tr.stats.starttime + 15000) #zoom in
df, npts = (tr.stats.sampling_rate, tr.stats.npts)

# Preprocessing
tr.data = tr.data.astype('float64') #convert to double
tr.data = detrend(tr.data, 'linear')
# Use cosinus Taper, percent of taper is currently 20%
win = cosTaper(npts, 0.2)
tr.data *= win
# Do the Fast Fourier Transformation for Real Valued Data
fdat = np.fft.rfft(tr.data, n=npts)
fdat /= abs(fdat).max() #normalize to 1

#
# The plotting part
#
freq = np.linspace(0,df/2,len(fdat))
time = np.arange(npts,dtype='f')/df
plt.figure(1)
plt.subplot(121)
plt.semilogx(freq, np.abs(fdat))
plt.ylabel('Norm. Amplitude')
plt.xlabel('Frequency [Hz]')
plt.subplot(122)
plt.plot(freq, np.angle(fdat))
plt.ylabel('Phase')
plt.xlabel('Freq [Hz]')

plt.figure(2)
plt.xlabel('Time [s]')
plt.ylabel('Normalized Counts')
plt.plot(time, tr.data/tr.data.max() , label='Raw Data')
plt.plot(time, win, label='Cosine Taper', linewidth=3 )
plt.title( "%s %s" % (tr.id, tr.stats.starttime))
plt.ylim(-1.2,1.2)
plt.legend()
plt.show()