# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Lomb-Scargle performance on Evenly Spaced Signal
#
# This notebook tests the performance of the Lomb-Scargle algorithm as implemented in Pyleoclim and therefore Scipy. In using this algorithm, we found spurious peaks at the Nyquist frequency on "perfect" sinusoidal signal. 
#
# An example of such spurious results can be found in [this notebook](https://github.com/KnowledgeCaptureAndDiscovery/autoTS/blob/master/notebooks/Methods/Lomb-Scargle%20Performance_RvsPython.ipynb). 
#
#

import pyleoclim as pyleo
import matplotlib.pyplot as plt
import scipy.signal as signal
import numpy as np

#
# ## Example from a unevenly-spaced signal
#
# ### Using the scipy implementation directly
#
# The code is taken from the [Scipy example](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lombscargle.html).

# +
A = 2.
w = 1.
phi = 0.5 * np.pi
nin = 1000
nout = 100000
frac_points = 0.9 # Fraction of points to select

r = np.random.rand(nin)
x = np.linspace(0.01, 10*np.pi, nin)
x = x[r >= frac_points]
y = A * np.sin(w*x+phi)

# Use the scipy.signal implementation on this signal
f = np.linspace(0.01, 10, nout)
pgram = signal.lombscargle(x, y, f, normalize=True)

#plot
plt.plot(f, pgram)
plt.show()
# -

# ### Using the Pyleoclim implementation
#
# The implementation in Pyleoclim uses windows (hanning by default) and WOSA (number of segments with 50% ovelap =3). The first step is to create a pyleoclim.Series object and apply the spectral method with parameter 'lomb_scargle'.
#
# First, let's use the frequency vector defined in the example from scipy, which is passed through the settings parameter.

ts=pyleo.Series(time=x,value=y)
psd=ts.spectral(method='lomb_scargle',settings={'freq':f})
psd.plot()

# Next, let's use the default frequency vector for this method which is based on the REDFIT implementation of Lomb-Scargle. 

psd2=ts.spectral(method='lomb_scargle', freq_method='lomb_scargle')
psd2.plot()

# There are no spurious peaks at the Nyquist frequency.
#
# ## Using an evenly-spaced signal
#
#

# +
t = np.arange(1,2001)
sig = np.cos(2*np.pi*1/20*t)

ts2=pyleo.Series(time=t,value=sig)
psd3=ts2.spectral(method='lomb_scargle', freq_method='lomb_scargle')
ts2.plot()
psd3.plot()
# -

# We find a spurious peak at the Nyquist frequency. Let's cut off earlier.

freq = np.linspace(1/5, 1/2000, nout)
psd4=ts2.spectral(method='lomb_scargle', settings={'freq':freq})
psd4.plot()

# There is still a spurious value at the last frequency, independent of Nyquist.
#
# Let's remove data points (40%) to our original signal.

# +
del_percent = 0.4
n_del = int(del_percent*np.size(t))
deleted_idx = np.random.choice(range(np.size(t)), n_del, replace=False)
time_unevenly =  np.delete(t, deleted_idx)
signal_unevenly =  np.delete(sig, deleted_idx)

ts3 = pyleo.Series(time=time_unevenly, value=signal_unevenly)
psd5=ts3.spectral(method='lomb_scargle', freq_method='lomb_scargle')
ts3.plot()
psd5.plot()
# -

# No spurious peak at the end but a less clearly defined peak and decay, which is to be expected from this degraded signal. The method is still capable of clearly identifyin the 20-year periodicity. 

# ### Scipy Implementation of the perfect harmonic

# +
t = np.arange(1,2001)
sig = np.cos(2*np.pi*1/20*t)

f=np.linspace(1/2000,1/2,5000)

pgram = signal.lombscargle(t,sig, f, normalize=True)

#plot
fig, ax = plt.subplots()
ax.plot(f, pgram_window)
ax.set_xscale('log', nonpositive='clip')
ax.set_yscale('log', nonpositive='clip')
plt.show()
# -

# ### Windowing effect
#

# +
pgram_window = signal.lombscargle(t,sig*signal.get_window('hann',len(t)), f)

#plot
fig, ax = plt.subplots()
ax.plot(f, pgram_window)
ax.set_xscale('log', nonpositive='clip')
ax.set_yscale('log', nonpositive='clip')
plt.show()
# -



StarcoderdataJupyterScriptsDedupFiltered