In [1]:
%load_ext autoreload
%autoreload 2
%reload_ext autoreload

In [2]:
from __future__ import division
import numpy as np 
from scipy import interpolate
from scipy.stats import pearsonr
import matplotlib.pyplot as plt
%matplotlib inline

In [11]:
basis = np.linspace(0,100,100)
data = 2 * np.sin(2*np.pi*basis/10) + 3 * np.sin(2*np.pi*basis/20) #+ 5* np.random.random(100)

In [12]:
signals = data[:, None]
base = len(basis)
min_env = np.zeros(base)
max_env = min_env.copy()
min_env = np.logical_and(
                    np.r_[True, signals[1:,0] < signals[:-1,0]],
                    np.r_[signals[:-1,0] < signals[1:,0], True])
max_env = np.logical_and(
                    np.r_[True, signals[1:,0] > signals[:-1,0]],
                    np.r_[signals[:-1,0] > signals[1:,0], True])
max_env[-1] = min_env[0] = False
min_env = min_env.nonzero()[0]
max_env = max_env.nonzero()[0]
signal_min = signals[min_env,0]
signal_max = signals[max_env,0]
order_min = 3
order_max = 3
inter_per = 0

In [15]:
n = 2
num = range(1, n+1)

min_arr = np.array([min_env, signal_min])
max_arr = np.array([max_env, signal_max])

left_min = np.zeros([2,n])
right_min = np.zeros([2, n])

left_max = np.zeros([2,n])
right_max = np.zeros([2, n])

for i in num:
    left_max[:, i-1] = [-1*min_env[n-i], signal_max[n-i]]
    left_min[:, i-1] = [-1*max_env[n-i], signal_min[n-i]]
    
    right_max[:, i-1] = [(basis[-1] - min_env[-i]) + basis[-1], signal_max[-i]]
    right_min[:, i-1] = [(basis[-1] - max_env[-i]) + basis[-1], signal_min[-i]]

    
ex_min_arr = np.concatenate([left_min, min_arr, right_min], axis=1)
ex_max_arr = np.concatenate([left_max, max_arr, right_max], axis=1)

In [16]:
plt.plot(basis,data)
plt.plot(ex_min_arr[0],ex_min_arr[1],'go')
plt.plot(ex_max_arr[0],ex_max_arr[1],'ro')


Out[16]:
[<matplotlib.lines.Line2D at 0x7f31c5cb1050>]

In [7]:
## NEW WAY
t = interpolate.splrep(*ex_min_arr, k=order_min, per=inter_per)
top = interpolate.splev(np.arange(len(signals[:,0])), t)
b = interpolate.splrep(*ex_max_arr, k=order_max, per=inter_per)
bot = interpolate.splev(np.arange(len(signals[:,0])), b)
env_new = (top+bot)/2
plt.plot(env_new)


Out[7]:
[<matplotlib.lines.Line2D at 0x7fcb6c7b6250>]

In [8]:
## OLD WAY
t = interpolate.splrep(min_env, signals[min_env,0],k=order_min, per=inter_per)
top = interpolate.splev(np.arange(len(signals[:,0])), t)
b = interpolate.splrep(max_env, signals[min_env,0], k=order_max, per=inter_per)
bot = interpolate.splev(np.arange(len(signals[:,0])), b)
env = (top+bot)/2
plt.plot(env)


Out[8]:
[<matplotlib.lines.Line2D at 0x7fcb6c6e4390>]

In [9]:
alpha = pearsonr(data,env_new)
alpha


Out[9]:
(0.82828781395462314, 2.1404776617442374e-26)

In [10]:
x = data - alpha[0]*env_new

In [11]:
plt.plot(data-env_new)
plt.plot(x,'r')


Out[11]:
[<matplotlib.lines.Line2D at 0x7fcb6c8804d0>]