kopia lustrzana https://github.com/gabrielegilardi/SignalFilters
complete comments
rodzic
753245f285
commit
87ef343228
|
@ -4,68 +4,61 @@ Signal Filtering/Smoothing and Generation of Synthetic Time-Series.
|
||||||
Copyright (c) 2020 Gabriele Gilardi
|
Copyright (c) 2020 Gabriele Gilardi
|
||||||
|
|
||||||
|
|
||||||
X (n_samples, n_series) Dataset to filter
|
X (n_samples, ) Dataset to filter (input)
|
||||||
b (n_b, ) Numerator coefficients
|
b (n_b, ) Transfer response coefficients (numerator)
|
||||||
a (n_a, ) Denominator coefficients
|
a (n_a, ) Transfer response coefficients (denominator)
|
||||||
Y (n_samples, n_series) Filtered dataset
|
Y (n_samples, ) Filtered dataset (output)
|
||||||
idx scalar First filtered element in Y
|
idx scalar First filtered element in Y
|
||||||
|
|
||||||
n_samples Number of data to filter
|
n_samples Number of samples in the input dataset
|
||||||
n_series Number of series to filter
|
nb Number of coefficients in array <b>
|
||||||
nb Number of coefficients (numerator)
|
na Number of coefficients in array <a>
|
||||||
na Number of coefficients (denominator)
|
|
||||||
|
|
||||||
Notes:
|
Notes:
|
||||||
- the filter is applied starting from index.
|
- the filter is applied starting from index idx = MAX(0, nb-1, na-1).
|
||||||
- non filtered data are set equal to the original input, i.e.
|
- non filtered data are set equal to the input, i.e. Y[0:idx-1] = X[0:idx-1]
|
||||||
Y[0:idx-1,:] = X[0:idx-1,:]
|
- X needs to be a 1D array.
|
||||||
- if n_series = 1 then must be ( ..., 1)
|
|
||||||
|
|
||||||
Filters:
|
|
||||||
|
|
||||||
Generic b,a Generic case
|
Filter list:
|
||||||
SMA N Simple Moving Average
|
-----------
|
||||||
EMA N/alpha Exponential Moving Average
|
Generic b, Generic
|
||||||
|
SMA N Simple moving average
|
||||||
|
EMA N/alpha Exponential moving average
|
||||||
WMA N Weighted moving average
|
WMA N Weighted moving average
|
||||||
MSMA N Modified Simple Moving Average
|
MSMA N Modified simple moving average
|
||||||
MLSQ N Modified Least-Squares Quadratic (N=5,7,9,11)
|
MLSQ N Modified least-squares quadratic (N = 5, 7, 9, 11)
|
||||||
ButterOrig P,N Butterworth original (N=2,3)
|
ButterOrig P, N Butterworth original filter (N = 2, 3)
|
||||||
ButterMod P,N Butterworth modified (N=2,3)
|
ButterMod P, N Butterworth modified filter (N = 2, 3)
|
||||||
SuperSmooth P,N Super smoother (N=2,3)
|
SuperSmooth P, N Supersmoother filter (N = 2, 3)
|
||||||
GaussLow P,N Gauss low pass (P>=2)
|
GaussLow P, N Gauss low pass filter (P > 1)
|
||||||
GaussHigh P,N Gauss high pass (P>=5)
|
GaussHigh P, N Gauss high pass filter (P > 4)
|
||||||
BandPass P,delta Band-pass filter
|
BandPass P, delta Band-pass filter
|
||||||
BandStop P,delta Band-stop filter
|
BandStop P, delta Band-stop filter
|
||||||
ZEMA1 N/alpha,K,Vn Zero-lag EMA (type 1)
|
ZEMA1 N/alpha, K, Vn Zero-lag EMA (type 1)
|
||||||
ZEMA2 N/alpha,K Zero-lag EMA (type 2)
|
ZEMA2 N/alpha, K Zero-lag EMA (type 2)
|
||||||
InstTrend N/alpha Instantaneous trendline
|
InstTrend N/alpha Instantaneous trendline
|
||||||
SincFunction N Sinc function
|
SincFilter P, nel Sinc function filter (N > 1)
|
||||||
Decycler P Decycler, 1-GaussHigh (P>=5)
|
Decycler P De-cycler filter (P >= 4)
|
||||||
DecyclerOsc P1,P2 Decycle oscillator, GH(P1) - GH(P2), (P1>=5)
|
DecyclerOsc P1, P2 De-cycle oscillator (P >= 4)
|
||||||
ABG
|
ABG alpha, beta, Alpha-beta-gamma filter (0 < alpha, beta < 1)
|
||||||
Kalman
|
gamma, dt
|
||||||
|
Kalman sigma_x, dt One-dimensional steady-state Kalman filter
|
||||||
|
sigma_v
|
||||||
|
|
||||||
N Order/smoothing factor/number of previous samples
|
N Order/smoothing factor/number of previous samples
|
||||||
alpha Damping term
|
alpha Damping term
|
||||||
P, P1, P2 Cut-off/critical period (50% power loss, -3 dB)
|
P, P1, P2 Cut-off/critical period (50% power loss, -3 dB)
|
||||||
delta Band centered in P and in fraction
|
delta Semi-band centered in P
|
||||||
(30% of P => 0.3, = 0.3*P, if P = 12 => 0.3*12 = 4)
|
|
||||||
K Coefficient/gain
|
K Coefficient/gain
|
||||||
Vn Look back sample (for the momentum)
|
Vn Look-back sample
|
||||||
|
nel Number of frequencies in the sinc function
|
||||||
correction = update = measurement
|
alpha Parameter(s) to correct the position in the ABG filter
|
||||||
prediction = motion
|
beta Parameter(s) to correct the velocity in the ABG filter
|
||||||
|
gamma Parameter(s) to correct the acceleration in the ABG filter
|
||||||
X (n_states, 1) State estimate
|
dt Sampling interval in the ABG and Kalman filters
|
||||||
P (n_states, n_states) Covariance estimate
|
sigma_x Process variance in the Kalman filter
|
||||||
F (n_states, n_states) State transition model
|
sigma_v Noise variance in the Kalman filter
|
||||||
Z (n_obs, 1) Observations
|
|
||||||
H (n_obs, n_states) Observation model
|
|
||||||
R (n_obs, n_obs) Covariance of the observation noise
|
|
||||||
S (n_obs, n_obs) Covariance of the observation residual
|
|
||||||
K (n_states, n_obs) Optimal Kalman gain
|
|
||||||
Q (n_states, n_states) Covariance of the process noise matrix
|
|
||||||
Y (n_obs, 1) Observation residual (innovation)
|
|
||||||
"""
|
"""
|
||||||
|
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
@ -73,19 +66,32 @@ from scipy import signal
|
||||||
import matplotlib.pyplot as plt
|
import matplotlib.pyplot as plt
|
||||||
|
|
||||||
|
|
||||||
def plot_signals(signals, start=0):
|
def plot_signals(signals, names=None, start=0):
|
||||||
"""
|
"""
|
||||||
signals must be a list
|
Plot the signals specified in list <signals> with their names specified in
|
||||||
|
list <names>. Each signal is plotted in its full length.
|
||||||
"""
|
"""
|
||||||
|
# Identify the signals by index if their name is not specified
|
||||||
|
if (names is None):
|
||||||
legend = []
|
legend = []
|
||||||
count = 0
|
count = 0
|
||||||
|
else:
|
||||||
|
legend = names
|
||||||
|
|
||||||
|
# Loop over the signals
|
||||||
for signal in signals:
|
for signal in signals:
|
||||||
|
|
||||||
signal = signal.flatten()
|
signal = signal.flatten()
|
||||||
end = len(signal)
|
end = len(signal)
|
||||||
t = np.arange(start, end)
|
t = np.arange(start, end)
|
||||||
plt.plot(t, signal[start:end])
|
plt.plot(t, signal[start:end])
|
||||||
|
|
||||||
|
# If no name is given use the list index to identify the signals
|
||||||
|
if (names is None):
|
||||||
legend.append('Signal [' + str(count) + ']')
|
legend.append('Signal [' + str(count) + ']')
|
||||||
count += 1
|
count += 1
|
||||||
|
|
||||||
|
# Plot and format
|
||||||
plt.xlabel('Index')
|
plt.xlabel('Index')
|
||||||
plt.ylabel('Value')
|
plt.ylabel('Value')
|
||||||
plt.grid(b=True)
|
plt.grid(b=True)
|
||||||
|
@ -95,23 +101,27 @@ def plot_signals(signals, start=0):
|
||||||
|
|
||||||
def filter_data(data, b, a):
|
def filter_data(data, b, a):
|
||||||
"""
|
"""
|
||||||
Applies a filter with transfer response coefficients <a> and <b>.
|
Applies a filter with transfer response coefficients <b> (numerator) and
|
||||||
|
<a> (denominator).
|
||||||
"""
|
"""
|
||||||
n_samples = len(data)
|
n_samples = len(data)
|
||||||
nb = len(b)
|
nb = len(b)
|
||||||
na = len(a)
|
na = len(a)
|
||||||
idx = np.amax([0, nb-1, na-1])
|
idx = np.amax([0, nb-1, na-1]) # Index of the 1st filtered sample
|
||||||
Y = data.copy()
|
Y = data.copy()
|
||||||
|
|
||||||
|
# Loop over the samples
|
||||||
for i in range(idx, n_samples):
|
for i in range(idx, n_samples):
|
||||||
|
|
||||||
tmp = 0
|
tmp = 0
|
||||||
|
|
||||||
|
# Contribution from the numerator term (input samples)
|
||||||
for j in range(nb):
|
for j in range(nb):
|
||||||
tmp += b[j] * data[i-j] # Numerator term
|
tmp += b[j] * data[i-j]
|
||||||
|
|
||||||
|
# Contribution from the denominator term (previous output samples)
|
||||||
for j in range(1, na):
|
for j in range(1, na):
|
||||||
tmp -= a[j] * Y[i-j] # Denominator term
|
tmp -= a[j] * Y[i-j]
|
||||||
|
|
||||||
Y[i] = tmp / a[0]
|
Y[i] = tmp / a[0]
|
||||||
|
|
||||||
|
@ -122,13 +132,16 @@ class Filter:
|
||||||
|
|
||||||
def __init__(self, data):
|
def __init__(self, data):
|
||||||
"""
|
"""
|
||||||
|
Initialize the filter object.
|
||||||
"""
|
"""
|
||||||
self.data = np.asarray(data)
|
self.data = np.asarray(data).flatten()
|
||||||
self.n_samples = len(data)
|
self.idx = 0
|
||||||
|
self.b = 0.0
|
||||||
|
self.a = 0.0
|
||||||
|
|
||||||
def Generic(self, b=1.0, a=1.0):
|
def Generic(self, b=1.0, a=1.0):
|
||||||
"""
|
"""
|
||||||
Filter with generic transfer response coefficients <a> and <b>.
|
Filter with generic transfer response coefficients <b> and <a>.
|
||||||
"""
|
"""
|
||||||
self.b = np.asarray(b)
|
self.b = np.asarray(b)
|
||||||
self.a = np.asarray(a)
|
self.a = np.asarray(a)
|
||||||
|
@ -136,9 +149,9 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def SMA(self, N=10):
|
def SMA(self, N=5):
|
||||||
"""
|
"""
|
||||||
Simple moving average (?? order, FIR, ?? band).
|
Simple moving average.
|
||||||
"""
|
"""
|
||||||
self.b = np.ones(N) / N
|
self.b = np.ones(N) / N
|
||||||
self.a = np.array([1.0])
|
self.a = np.array([1.0])
|
||||||
|
@ -146,9 +159,10 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def EMA(self, N=10, alpha=None):
|
def EMA(self, N=5, alpha=None):
|
||||||
"""
|
"""
|
||||||
Exponential moving average (?? order, IIR, pass ??).
|
Exponential moving average.
|
||||||
|
|
||||||
If not given, <alpha> is determined as equivalent to a N-SMA.
|
If not given, <alpha> is determined as equivalent to a N-SMA.
|
||||||
"""
|
"""
|
||||||
if (alpha is None):
|
if (alpha is None):
|
||||||
|
@ -160,10 +174,11 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def WMA(self, N=10):
|
def WMA(self, N=5):
|
||||||
"""
|
"""
|
||||||
Weighted moving average (?? order, FIR, pass ??).
|
Weighted moving average.
|
||||||
Example: N = 5 --> [5.0, 4.0, 3.0, 2.0, 1.0] / 15.0
|
|
||||||
|
Example: N = 5 --> [5, 4, 3, 2, 1] / 15.
|
||||||
"""
|
"""
|
||||||
w = np.arange(N, 0, -1)
|
w = np.arange(N, 0, -1)
|
||||||
|
|
||||||
|
@ -173,10 +188,11 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def MSMA(self, N=10):
|
def MSMA(self, N=5):
|
||||||
"""
|
"""
|
||||||
Modified simple moving average (?? order, FIR, pass ??).
|
Modified simple moving average.
|
||||||
Example: N = 4 --> [0.5, 1.0, 1.0, 1.0, 0.5] / 4.0
|
|
||||||
|
Example: N = 5 --> [1/2, 1, 1, 1, 1, 1/2] / 5.
|
||||||
"""
|
"""
|
||||||
w = np.ones(N+1)
|
w = np.ones(N+1)
|
||||||
w[0] = 0.5
|
w[0] = 0.5
|
||||||
|
@ -190,26 +206,26 @@ class Filter:
|
||||||
|
|
||||||
def MLSQ(self, N=5):
|
def MLSQ(self, N=5):
|
||||||
"""
|
"""
|
||||||
Modified simple moving average (?? order, FIR, pass ??).
|
Modified least-squares quadratic.
|
||||||
Only N = 5, 7, 9, and 11 are implemented. If not returns the unfiltered
|
|
||||||
dataset.
|
Must be N = 5, 7, 9, or 11. If wrong N, prints a warning and returns
|
||||||
|
the unfiltered dataset.
|
||||||
"""
|
"""
|
||||||
if (N == 5):
|
if (N == 5):
|
||||||
w = np.array([7.0, 24.0, 34.0, 24.0, 7.0]) / 96.0
|
w = np.array([7, 24, 34, 24, 7]) / 96
|
||||||
|
|
||||||
elif (N == 7):
|
elif (N == 7):
|
||||||
w = np.array([1.0, 6.0, 12.0, 14.0, 12.0, 6.0, 1.0]) / 52.0
|
w = np.array([1, 6, 12, 14, 12, 6, 1]) / 52
|
||||||
|
|
||||||
elif (N == 9):
|
elif (N == 9):
|
||||||
w = np.array([-1.0, 28.0, 78.0, 108.0, 118.0, 108.0, 78.0, 28.0,
|
w = np.array([-1, 28, 78, 108, 118, 108, 78, 28, -1]) / 544
|
||||||
-1.0]) / 544.0
|
|
||||||
|
|
||||||
elif (N == 11):
|
elif (N == 11):
|
||||||
w = np.array([-11.0, 18.0, 88.0, 138.0, 168.0, 178.0, 168.0, 138.0,
|
w = np.array([-11, 18, 88, 138, 168, 178, 168, 138, 88, 18,
|
||||||
88.0, 18.0, -11.0]) / 980.0
|
-11]) / 980
|
||||||
|
|
||||||
else:
|
else:
|
||||||
print("Warning: data returned unfiltered (wrong N)")
|
print("Warning: data returned unfiltered (MLSQ - Wrong N)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -221,8 +237,10 @@ class Filter:
|
||||||
|
|
||||||
def ButterOrig(self, N=2, P=10):
|
def ButterOrig(self, N=2, P=10):
|
||||||
"""
|
"""
|
||||||
Butterworth original version (?? order, IIR, pass ??).
|
Butterworth original filter.
|
||||||
Only N = 2 and 3 are implemented. If not returns the unfiltered dataset.
|
|
||||||
|
Must be N = 2 or 3. If wrong N, prints a warning and returns the
|
||||||
|
unfiltered dataset.
|
||||||
"""
|
"""
|
||||||
if (N == 2):
|
if (N == 2):
|
||||||
beta = np.exp(-np.sqrt(2.0) * np.pi / P)
|
beta = np.exp(-np.sqrt(2.0) * np.pi / P)
|
||||||
|
@ -239,7 +257,7 @@ class Filter:
|
||||||
(1.0 + alpha) * beta ** 2.0, - beta ** 4.0])
|
(1.0 + alpha) * beta ** 2.0, - beta ** 4.0])
|
||||||
|
|
||||||
else:
|
else:
|
||||||
print("Warning: data returned unfiltered (wrong N)")
|
print("Warning: data returned unfiltered (ButterOrig - Wrong N)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -251,8 +269,11 @@ class Filter:
|
||||||
|
|
||||||
def ButterMod(self, N=2, P=10):
|
def ButterMod(self, N=2, P=10):
|
||||||
"""
|
"""
|
||||||
Butterworth modified version (?? order, IIR, pass ??).
|
Butterworth modified filter. It is derived from the Butterworth original
|
||||||
Only N = 2 and 3 are implemented. If not returns the unfiltered dataset.
|
filter deleting all but the constant term at the numerator.
|
||||||
|
|
||||||
|
Must be N = 2 or 3. If wrong N, prints a warning and returns the
|
||||||
|
unfiltered dataset.
|
||||||
"""
|
"""
|
||||||
if (N == 2):
|
if (N == 2):
|
||||||
beta = np.exp(-np.sqrt(2.0) * np.pi / P)
|
beta = np.exp(-np.sqrt(2.0) * np.pi / P)
|
||||||
|
@ -268,7 +289,7 @@ class Filter:
|
||||||
(1.0 + alpha) * beta ** 2.0, - beta ** 4.0])
|
(1.0 + alpha) * beta ** 2.0, - beta ** 4.0])
|
||||||
|
|
||||||
else:
|
else:
|
||||||
print("Warning: data returned unfiltered (wrong N)")
|
print("Warning: data returned unfiltered (ButterMod - Wrong N)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -280,8 +301,11 @@ class Filter:
|
||||||
|
|
||||||
def SuperSmooth(self, N=2, P=10):
|
def SuperSmooth(self, N=2, P=10):
|
||||||
"""
|
"""
|
||||||
SuperSmooth (?? order, IIR, pass ??).
|
Supersmoother filter. It is derived from the Butterworth modified
|
||||||
Only N = 2 and 3 are implemented. If not returns the unfiltered dataset.
|
filter adding a two-element moving average at the numerator.
|
||||||
|
|
||||||
|
Must be N = 2 or 3. If wrong N, prints a warning and returns the
|
||||||
|
unfiltered dataset.
|
||||||
"""
|
"""
|
||||||
if (N == 2):
|
if (N == 2):
|
||||||
beta = np.exp(-np.sqrt(2.0) * np.pi / P)
|
beta = np.exp(-np.sqrt(2.0) * np.pi / P)
|
||||||
|
@ -298,7 +322,7 @@ class Filter:
|
||||||
(1.0 + alpha) * beta ** 2.0, - beta ** 4.0])
|
(1.0 + alpha) * beta ** 2.0, - beta ** 4.0])
|
||||||
|
|
||||||
else:
|
else:
|
||||||
print("Warning: data returned unfiltered (wrong N)")
|
print("Warning: data returned unfiltered (SuperSmooth - Wrong N)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -308,13 +332,15 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def GaussLow(self, N=1, P=2):
|
def GaussLow(self, N=1, P=10):
|
||||||
"""
|
"""
|
||||||
Gauss low pass (IIR, N-th order, low pass).
|
Gauss low pass filter.
|
||||||
Must be P > 1. If not returns the unfiltered dataset.
|
|
||||||
|
Must be P > 1. If wrong P, prints a warning and returns the unfiltered
|
||||||
|
dataset.
|
||||||
"""
|
"""
|
||||||
if (P < 2):
|
if (P <= 1):
|
||||||
print("Warning: data returned unfiltered (P < 2)")
|
print("Warning: data returned unfiltered (GaussLow - Wrong P)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -331,13 +357,15 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def GaussHigh(self, N=1, P=5):
|
def GaussHigh(self, N=1, P=10):
|
||||||
"""
|
"""
|
||||||
Gauss high pass (IIR, Nth order, high pass).
|
Gauss high pass filter.
|
||||||
Must be P > 4. If not returns the unfiltered dataset.
|
|
||||||
|
Must be P > 4. If wrong P, prints a warning and returns the unfiltered
|
||||||
|
dataset.
|
||||||
"""
|
"""
|
||||||
if (P < 5):
|
if (P <= 4):
|
||||||
print("Warning: data returned unfiltered (P < 5)")
|
print("Warning: data returned unfiltered (GaussHigh - Wrong P)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -354,11 +382,11 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def BandPass(self, P=5, delta=0.3):
|
def BandPass(self, P=10, delta=0.3):
|
||||||
"""
|
"""
|
||||||
Band-pass (type, order, IIR).
|
Band-pass filter.
|
||||||
Example: delta = 0.3, P = 12
|
|
||||||
(30% of P => 0.3, = 0.3*P, if P = 12 => 0.3*12 = 4)
|
Example: delta = 0.3, P = 10 --> 0.3 * 10 = 3 --> band is [7, 13]
|
||||||
"""
|
"""
|
||||||
beta = np.cos(2.0 * np.pi / P)
|
beta = np.cos(2.0 * np.pi / P)
|
||||||
gamma = np.cos(4.0 * np.pi * delta / P)
|
gamma = np.cos(4.0 * np.pi * delta / P)
|
||||||
|
@ -372,9 +400,9 @@ class Filter:
|
||||||
|
|
||||||
def BandStop(self, P=5, delta=0.3):
|
def BandStop(self, P=5, delta=0.3):
|
||||||
"""
|
"""
|
||||||
Band-stop (type, order, IIR)
|
Band-stop filter.
|
||||||
Example: delta = 0.3, P = 12
|
|
||||||
(30% of P => 0.3, = 0.3*P, if P = 12 => 0.3*12 = 4)
|
Example: delta = 0.3, P = 10 --> 0.3 * 10 = 3 --> band is [7, 13]
|
||||||
"""
|
"""
|
||||||
beta = np.cos(2.0 * np.pi / P)
|
beta = np.cos(2.0 * np.pi / P)
|
||||||
gamma = np.cos(4.0 * np.pi * delta / P)
|
gamma = np.cos(4.0 * np.pi * delta / P)
|
||||||
|
@ -388,13 +416,15 @@ class Filter:
|
||||||
|
|
||||||
def ZEMA1(self, N=10, alpha=None, K=1.0, Vn=5):
|
def ZEMA1(self, N=10, alpha=None, K=1.0, Vn=5):
|
||||||
"""
|
"""
|
||||||
Zero lag Exponential Moving Average (type 1).
|
Zero-lag EMA (type 1). It is an alpha-beta type filter with sub-optimal
|
||||||
|
parameters.
|
||||||
|
|
||||||
If not given, <alpha> is determined as equivalent to a N-SMA.
|
If not given, <alpha> is determined as equivalent to a N-SMA.
|
||||||
"""
|
"""
|
||||||
if (alpha is None):
|
if (alpha is None):
|
||||||
alpha = 2.0 / (N + 1.0)
|
alpha = 2.0 / (N + 1.0)
|
||||||
|
|
||||||
w = np.zeros(Vn+1)
|
w = np.zeros(Vn + 1)
|
||||||
w[0] = alpha * (1.0 + K)
|
w[0] = alpha * (1.0 + K)
|
||||||
w[Vn] = - alpha * K
|
w[Vn] = - alpha * K
|
||||||
|
|
||||||
|
@ -406,7 +436,9 @@ class Filter:
|
||||||
|
|
||||||
def ZEMA2(self, N=10, alpha=None, K=1.0):
|
def ZEMA2(self, N=10, alpha=None, K=1.0):
|
||||||
"""
|
"""
|
||||||
Zero lag Exponential Moving Average (type 2).
|
Zero-lag EMA (type 2). It is derived from the type 1 ZEMA removing the
|
||||||
|
look-back term Vn.
|
||||||
|
|
||||||
If not given, <alpha> is determined as equivalent to a N-SMA.
|
If not given, <alpha> is determined as equivalent to a N-SMA.
|
||||||
"""
|
"""
|
||||||
if (alpha is None):
|
if (alpha is None):
|
||||||
|
@ -420,7 +452,9 @@ class Filter:
|
||||||
|
|
||||||
def InstTrend(self, N=10, alpha=None):
|
def InstTrend(self, N=10, alpha=None):
|
||||||
"""
|
"""
|
||||||
Instantaneous Trendline (2nd order, IIR, low pass).
|
Instantaneous Trendline. It is created by removing the dominant cycle
|
||||||
|
from the signal.
|
||||||
|
|
||||||
If not given, <alpha> is determined as equivalent to a N-SMA.
|
If not given, <alpha> is determined as equivalent to a N-SMA.
|
||||||
"""
|
"""
|
||||||
if (alpha is None):
|
if (alpha is None):
|
||||||
|
@ -433,15 +467,22 @@ class Filter:
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def SincFunction(self, N=10, nel=10):
|
def SincFilter(self, P=10, nel=10):
|
||||||
"""
|
"""
|
||||||
Sinc function (order, FIR, pass).
|
Sinc function filter. The cut off point is at 0.5/P.
|
||||||
(N > 1, cut off at 0.5/N)
|
|
||||||
|
Must be P > 1. If wrong P, prints a warning and returns the unfiltered
|
||||||
|
dataset.
|
||||||
"""
|
"""
|
||||||
|
if (P <= 1):
|
||||||
|
print("Warning: data returned unfiltered (SincFilter - Wrong P)")
|
||||||
|
self.idx = 0
|
||||||
|
return self.data
|
||||||
|
|
||||||
K = np.arange(1, nel)
|
K = np.arange(1, nel)
|
||||||
w = np.zeros(nel)
|
w = np.zeros(nel)
|
||||||
w[0] = 1.0 / N
|
w[0] = 1.0 / N
|
||||||
w[1:] = np.sin(np.pi * K / N) / (np.pi * K)
|
w[1:] = np.sin(np.pi * K / P) / (np.pi * K)
|
||||||
|
|
||||||
self.b = w
|
self.b = w
|
||||||
self.a = np.array([1.0])
|
self.a = np.array([1.0])
|
||||||
|
@ -451,12 +492,14 @@ class Filter:
|
||||||
|
|
||||||
def Decycler(self, P=10):
|
def Decycler(self, P=10):
|
||||||
"""
|
"""
|
||||||
Decycler (?? order, IIR, pass ??). Gauss,HP,1st,P
|
De-cycler filter. It is derived subtracting a 1st order high pass Gauss
|
||||||
Built subtracting high pass Gauss filter from 1 (order 1)
|
filter from 1.
|
||||||
Must be P > 4. If not returns the unfiltered dataset.
|
|
||||||
|
Must be P > 4. If wrong P, prints a warning and returns the unfiltered
|
||||||
|
dataset.
|
||||||
"""
|
"""
|
||||||
if (P < 5):
|
if (P <= 4):
|
||||||
print("Warning: data returned unfiltered (P < 5)")
|
print("Warning: data returned unfiltered (Decycler - Wrong P)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -465,17 +508,18 @@ class Filter:
|
||||||
|
|
||||||
def DecyclerOsc(self, P1=5, P2=10):
|
def DecyclerOsc(self, P1=5, P2=10):
|
||||||
"""
|
"""
|
||||||
DecyclerOsc (?? order 2, IIR, pass ??).
|
De-cycler oscillator. It is derived subtracting a 2nd order high pass
|
||||||
|
Gauss filter with higher cut-off period from a 2nd order high pass Gauss
|
||||||
|
filter with higher cut-off period.
|
||||||
|
|
||||||
(Gauss, HP, 2nd order, Pmax - Gauss, HP, 2nd order, Pmin)
|
Must be P > 4. If wrong P, prints a warning and returns the unfiltered
|
||||||
P1 = 1st cut off period, P2 = 2nd cut off period. Automatically fixed.
|
dataset.
|
||||||
Must be P1, P2 > 4. If not returns the unfiltered dataset.
|
|
||||||
"""
|
"""
|
||||||
P_low = np.amin([P1, P2])
|
P_low = np.amin([P1, P2])
|
||||||
P_high = np.amax([P1, P2])
|
P_high = np.amax([P1, P2])
|
||||||
|
|
||||||
if (P_low < 5):
|
if (P_low <= 4):
|
||||||
print("Warning: data returned unfiltered (P_low < 5)")
|
print("Warning: data returned unfiltered (DecyclerOsc - Wrong P)")
|
||||||
self.idx = 0
|
self.idx = 0
|
||||||
return self.data
|
return self.data
|
||||||
|
|
||||||
|
@ -484,34 +528,38 @@ class Filter:
|
||||||
|
|
||||||
def ABG(self, alpha=0.0, beta=0.0, gamma=0.0, dt=1.0):
|
def ABG(self, alpha=0.0, beta=0.0, gamma=0.0, dt=1.0):
|
||||||
"""
|
"""
|
||||||
alpha-beta-gamma
|
Alpha-beta-gamma filter. It is a predictor-corrector type of filter.
|
||||||
For numerical stability: 0 < alpha, beta < 1
|
|
||||||
|
Arguments alpha, beta, and gamma can be a scalar (used for all samples)
|
||||||
|
or an array with one value for each sample. For numerical stability it
|
||||||
|
should be 0 < alpha, beta < 1.
|
||||||
"""
|
"""
|
||||||
# If necessary change scalars to arrays
|
n_samples = len(data)
|
||||||
|
Y = np.zeros(n_samples)
|
||||||
|
|
||||||
|
# Change scalar arguments to arrays if necessary
|
||||||
if (np.ndim(alpha) == 0):
|
if (np.ndim(alpha) == 0):
|
||||||
alpha = np.ones(self.n_samples) * alpha
|
alpha = np.ones(n_samples) * alpha
|
||||||
if (np.ndim(beta) == 0):
|
if (np.ndim(beta) == 0):
|
||||||
beta = np.ones(self.n_samples) * beta
|
beta = np.ones(n_samples) * beta
|
||||||
if (np.ndim(gamma) == 0):
|
if (np.ndim(gamma) == 0):
|
||||||
gamma = np.ones(self.n_samples) * gamma
|
gamma = np.ones(n_samples) * gamma
|
||||||
|
|
||||||
# Initialize
|
# Initialize
|
||||||
Y_corr = self.data.copy()
|
x0 = self.data[0]
|
||||||
Y_pred = self.data.copy()
|
v0 = 0.0
|
||||||
x0 = self.data[0, :]
|
a0 = 0.0
|
||||||
v0 = np.zeros(self.n_series)
|
Y[0] = x0
|
||||||
a0 = np.zeros(self.n_series)
|
|
||||||
|
|
||||||
for i in range(1, self.n_samples):
|
for i in range(1, n_samples):
|
||||||
|
|
||||||
# Predictor (predicts state in <i>)
|
# Predictor (predicts state in <i>)
|
||||||
x_pred = x0 + dt * v0 + 0.5 * a0 * dt ** 2.0
|
x_pred = x0 + dt * v0 + 0.5 * a0 * dt ** 2.0
|
||||||
v_pred = v0 + dt * a0
|
v_pred = v0 + dt * a0
|
||||||
a_pred = a0
|
a_pred = a0
|
||||||
Y_pred[i, :] = x_pred
|
|
||||||
|
|
||||||
# Residual (innovation)
|
# Residual (innovation)
|
||||||
r = self.data[i, :] - x_pred
|
r = self.data[i] - x_pred
|
||||||
|
|
||||||
# Corrector (corrects state in <i>)
|
# Corrector (corrects state in <i>)
|
||||||
x_corr = x_pred + alpha[i] * r
|
x_corr = x_pred + alpha[i] * r
|
||||||
|
@ -519,18 +567,23 @@ class Filter:
|
||||||
a_corr = a_pred + (2.0 * gamma[i] / dt ** 2.0) * r
|
a_corr = a_pred + (2.0 * gamma[i] / dt ** 2.0) * r
|
||||||
|
|
||||||
# Save value and prepare next iteration
|
# Save value and prepare next iteration
|
||||||
Y_corr[i, :] = x_corr
|
|
||||||
x0 = x_corr
|
x0 = x_corr
|
||||||
v0 = v_corr
|
v0 = v_corr
|
||||||
a0 = a_corr
|
a0 = a_corr
|
||||||
|
Y[i] = x_corr
|
||||||
|
|
||||||
self.idx = 1
|
self.idx = 1
|
||||||
|
|
||||||
return Y_corr, Y_pred
|
return Y
|
||||||
|
|
||||||
def Kalman(self, sigma_x, sigma_v, dt, abg_type="abg"):
|
def Kalman(self, sigma_x, sigma_v, dt, abg_type="abg"):
|
||||||
"""
|
"""
|
||||||
Steady-state Kalman filter (also limited to one-dimension)
|
One-dimensional steady-state Kalman filter. It is obtained from the
|
||||||
|
alpha-beta-gamma filter using the process variance, the noise variance
|
||||||
|
and optimizing the three parameters.
|
||||||
|
|
||||||
|
Arguments sigma_x and sigma_v can be a scalar (used for all samples) or
|
||||||
|
an array with one value for each sample.
|
||||||
"""
|
"""
|
||||||
L = (sigma_x / sigma_v) * dt ** 2.0
|
L = (sigma_x / sigma_v) * dt ** 2.0
|
||||||
|
|
||||||
|
@ -561,19 +614,23 @@ class Filter:
|
||||||
beta = 2.0 * (1 - s) ** 2.0
|
beta = 2.0 * (1 - s) ** 2.0
|
||||||
gamma = (beta ** 2.0) / (2.0 * alpha)
|
gamma = (beta ** 2.0) / (2.0 * alpha)
|
||||||
|
|
||||||
# Apply filter
|
# Apply the alpha-beta-gamma filter
|
||||||
Y = self.abg(alpha=alpha, beta=beta, gamma=gamma, dt=dt)
|
Y = self.abg(alpha=alpha, beta=beta, gamma=gamma, dt=dt)
|
||||||
|
|
||||||
return Y
|
return Y
|
||||||
|
|
||||||
def plot_frequency(self):
|
def plot_frequency(self):
|
||||||
"""
|
"""
|
||||||
|
Plots the frequency response (in decibels) of the filter with transfer
|
||||||
|
response coefficients <b> and <a>.
|
||||||
"""
|
"""
|
||||||
w, h = signal.freqz(self.b, self.a)
|
w, h = signal.freqz(self.b, self.a)
|
||||||
h_db = 20.0 * np.log10(np.abs(h))
|
h_db = 20.0 * np.log10(np.abs(h)) # Convert to decibels
|
||||||
wf = w / (2.0 * np.pi)
|
wf = w / (2.0 * np.pi) # Scale to [0, 0.5]
|
||||||
|
|
||||||
|
# Plot and format
|
||||||
plt.plot(wf, h_db)
|
plt.plot(wf, h_db)
|
||||||
plt.axhline(-3.0, lw=1.5, ls='--', C='r')
|
plt.axhline(-3.0, lw=1.5, ls='--', C='r') # -3 dB (50% power loss)
|
||||||
plt.grid(b=True)
|
plt.grid(b=True)
|
||||||
plt.xlim(np.amin(wf), np.amax(wf))
|
plt.xlim(np.amin(wf), np.amax(wf))
|
||||||
plt.xlabel(r'$\omega$ [rad/sample]')
|
plt.xlabel(r'$\omega$ [rad/sample]')
|
||||||
|
@ -584,9 +641,13 @@ class Filter:
|
||||||
|
|
||||||
def plot_lag(self):
|
def plot_lag(self):
|
||||||
"""
|
"""
|
||||||
|
Plots the lag (group delay) of the filter with transfer response
|
||||||
|
coefficients <b> and <a>.
|
||||||
"""
|
"""
|
||||||
w, gd = signal.group_delay((self.b, self.a))
|
w, gd = signal.group_delay((self.b, self.a))
|
||||||
wf = w / (2.0 * np.pi)
|
wf = w / (2.0 * np.pi) # Scale to [0, 0.5]
|
||||||
|
|
||||||
|
# Plot and format
|
||||||
plt.plot(wf, gd)
|
plt.plot(wf, gd)
|
||||||
plt.grid(b=True)
|
plt.grid(b=True)
|
||||||
plt.xlim(np.amin(wf), np.amax(wf))
|
plt.xlim(np.amin(wf), np.amax(wf))
|
||||||
|
|
|
@ -49,14 +49,12 @@ def synthetic_wave(P, A=None, phi=None, num=1000):
|
||||||
def synthetic_FFT(X, n_reps=1):
|
def synthetic_FFT(X, n_reps=1):
|
||||||
"""
|
"""
|
||||||
Generates surrogates of the time-serie X using the phase-randomized
|
Generates surrogates of the time-serie X using the phase-randomized
|
||||||
Fourier-transform algorithm.
|
Fourier-transform algorithm. Input X needs to be a 1D array.
|
||||||
|
|
||||||
X (n, ) Original time-series
|
X (n, ) Original time-series
|
||||||
X_fft (n, ) FFT of the original time-series
|
X_fft (n, ) FFT of the original time-series
|
||||||
X_synt_fft (n_reps, n) FFT of the synthetic time-series
|
X_synt_fft (n_reps, n) FFT of the synthetic time-series
|
||||||
X_synt (n_reps, n) Synthetic time-series
|
X_synt (n_reps, n) Synthetic time-series
|
||||||
|
|
||||||
Input array X needs to be a 1D array (of any shape).
|
|
||||||
"""
|
"""
|
||||||
X = X.flatten() # Reshape to (n, )
|
X = X.flatten() # Reshape to (n, )
|
||||||
n = len(X)
|
n = len(X)
|
||||||
|
@ -95,13 +93,11 @@ def synthetic_FFT(X, n_reps=1):
|
||||||
def synthetic_sampling(X, n_reps=1, replace=True):
|
def synthetic_sampling(X, n_reps=1, replace=True):
|
||||||
"""
|
"""
|
||||||
Generates surrogates of the time-serie X using randomized-sampling
|
Generates surrogates of the time-serie X using randomized-sampling
|
||||||
(bootstrap) with or without replacement.
|
(bootstrap) with or without replacement. Input X needs to be a 1D array.
|
||||||
|
|
||||||
X (n, ) Original time-series
|
X (n, ) Original time-series
|
||||||
idx (n_reps, n) Random index of X
|
idx (n_reps, n) Random index of X
|
||||||
X_synt (n_reps, n) Synthetic time-series
|
X_synt (n_reps, n) Synthetic time-series
|
||||||
|
|
||||||
Input array X needs to be a 1D array (of any shape).
|
|
||||||
"""
|
"""
|
||||||
X = X.flatten() # Reshape to (n, )
|
X = X.flatten() # Reshape to (n, )
|
||||||
n = len(X)
|
n = len(X)
|
||||||
|
@ -123,7 +119,7 @@ def synthetic_sampling(X, n_reps=1, replace=True):
|
||||||
def synthetic_MEboot(X, n_reps=1, alpha=0.1, bounds=False, scale=False):
|
def synthetic_MEboot(X, n_reps=1, alpha=0.1, bounds=False, scale=False):
|
||||||
"""
|
"""
|
||||||
Generates surrogates of the time-serie X using the maximum entropy
|
Generates surrogates of the time-serie X using the maximum entropy
|
||||||
bootstrap algorithm.
|
bootstrap algorithm. Input X needs to be a 1D array.
|
||||||
|
|
||||||
X (n, ) Original time-series
|
X (n, ) Original time-series
|
||||||
idx (n, ) Original order of X
|
idx (n, ) Original order of X
|
||||||
|
@ -135,8 +131,6 @@ def synthetic_MEboot(X, n_reps=1, alpha=0.1, bounds=False, scale=False):
|
||||||
w_corr (n_reps, n) Interpolated new points with corrections for first
|
w_corr (n_reps, n) Interpolated new points with corrections for first
|
||||||
and last interval
|
and last interval
|
||||||
X_synt (n_reps, n) Synthetic time-series
|
X_synt (n_reps, n) Synthetic time-series
|
||||||
|
|
||||||
Input array X needs to be a 1D array (of any shape).
|
|
||||||
"""
|
"""
|
||||||
X = X.flatten() # Reshape to (n, )
|
X = X.flatten() # Reshape to (n, )
|
||||||
n = len(X)
|
n = len(X)
|
||||||
|
@ -259,13 +253,13 @@ def value2diff(X, percent=True):
|
||||||
"""
|
"""
|
||||||
Returns the 1st discrete difference of array X.
|
Returns the 1st discrete difference of array X.
|
||||||
|
|
||||||
X (n, ) Original dataset
|
X (n, ) Input dataset
|
||||||
dX (n-1, ) 1st discrete differences
|
dX (n-1, ) 1st discrete differences
|
||||||
|
|
||||||
Notes:
|
Notes:
|
||||||
- the discrete difference can be calculated in percent or in value.
|
- the discrete difference can be calculated in percent or in value.
|
||||||
- array dX is one element shorter than array X.
|
- dX is one element shorter than X.
|
||||||
- array X needs to be a 1D array (of any shape).
|
- X needs to be a 1D array.
|
||||||
"""
|
"""
|
||||||
X = X.flatten() # Reshape to (n, )
|
X = X.flatten() # Reshape to (n, )
|
||||||
|
|
||||||
|
@ -282,36 +276,39 @@ def value2diff(X, percent=True):
|
||||||
|
|
||||||
def diff2value(dX, X0, percent=True):
|
def diff2value(dX, X0, percent=True):
|
||||||
"""
|
"""
|
||||||
Returns array X from the 1st discrete difference.
|
Returns array X from the 1st discrete difference using X0 as initial value.
|
||||||
|
|
||||||
dX (n, ) Discrete differences
|
dX (n, ) Discrete differences
|
||||||
X0 scalar Initial value
|
X0 scalar Initial value
|
||||||
X (n+1, ) Original dataset ?????
|
X (n+1, ) Output dataset
|
||||||
|
|
||||||
Notes:
|
Notes:
|
||||||
- the discrete difference can be in percent or in value.
|
- the discrete difference can be in percent or in value.
|
||||||
- array X is one element longer than array dX.
|
- X is one element longer than dX.
|
||||||
- array dX needs to be a 1D array (of any shape).
|
- dX needs to be a 1D array.
|
||||||
|
|
||||||
If the discrete difference is in percent, the first column of X is set to
|
If the discrete difference is in percent:
|
||||||
one. The original array is X[0] * X
|
X[0] = X0
|
||||||
|
X[1] = X[0] * (1 + dX[0])
|
||||||
|
X[2] = X[1] * (1 + dX[1]) = X[0] * (1 + dX[0]) * (1 + dX[1])
|
||||||
|
....
|
||||||
|
|
||||||
If the discrete difference is in value, the first column of X is set to
|
If the discrete difference is in value:
|
||||||
zero. X0+X
|
X[0] = X0
|
||||||
- array X needs to be a 1D array (of any shape).
|
X[1] = X[0] + dX[0]
|
||||||
|
X[2] = X[1] + dX[1] = X[0] + dX[0] + dX[1]
|
||||||
|
....
|
||||||
"""
|
"""
|
||||||
dX = dX.flatten() # Reshape to (n, )
|
dX = dX.flatten() # Reshape to (n, )
|
||||||
n = len(dX)
|
X = np.zeros(len(dX) + 1)
|
||||||
X = np.zeros(n+1)
|
X[0] = X0 # Initial value
|
||||||
|
|
||||||
# Discrete difference in percent
|
# Discrete difference in percent
|
||||||
if (percent):
|
if (percent):
|
||||||
X[0] = X0
|
X[1:] = X0 * np.cumprod(1.0 + dX)
|
||||||
X[1:] = np.cumprod(1.0 + dX) * X0 # !!!! check
|
|
||||||
|
|
||||||
# Discrete difference in value
|
# Discrete difference in value
|
||||||
else:
|
else:
|
||||||
X[0] = X0
|
X[1:] = X0 + np.cumsum(dX)
|
||||||
X[1:] = np.cumsum(dX) + X0 # !!!!!c check
|
|
||||||
|
|
||||||
return X
|
return X
|
||||||
|
|
|
@ -5,10 +5,10 @@ Copyright (c) 2020 Gabriele Gilardi
|
||||||
|
|
||||||
|
|
||||||
ToDo:
|
ToDo:
|
||||||
- add comments to the code
|
- do examples and check type (low pass, etc)
|
||||||
- in comments write what filters do
|
- put type in description
|
||||||
- is necessary to copy X for Y untouched?
|
- see paper ZEMA for tests
|
||||||
- decide default values in functions
|
- check results with ML
|
||||||
"""
|
"""
|
||||||
|
|
||||||
import sys
|
import sys
|
||||||
|
@ -31,18 +31,21 @@ data_file = sys.argv[1] + '.csv'
|
||||||
|
|
||||||
# Read data from a csv file (one time-series each column)
|
# Read data from a csv file (one time-series each column)
|
||||||
data = np.loadtxt(data_file, delimiter=',')
|
data = np.loadtxt(data_file, delimiter=',')
|
||||||
|
print(data.shape)
|
||||||
|
|
||||||
t, f = syn.synthetic_wave([1., 2., 3.], A=None, phi=None, num=1000)
|
# t, f = syn.synthetic_wave([1., 2., 3.], A=None, phi=None, num=1000)
|
||||||
plt.plot(t,f)
|
# plt.plot(t,f)
|
||||||
plt.show()
|
# plt.show()
|
||||||
|
|
||||||
# spx = flt.Filter(data)
|
spx = flt.Filter(data)
|
||||||
# res = spx.EMA(N=10)
|
# res = spx.EMA(N=10)
|
||||||
# signals = [spx.data, res[0:400]]
|
# signals = [spx.data, res[0:400]]
|
||||||
# flt.plot_signals(signals, start=100)
|
# flt.plot_signals(signals, ['SPX', 'SMA'])
|
||||||
|
|
||||||
# spx.plot_frequency()
|
res = spx.BandPass(P=10, delta=0.3)
|
||||||
# spx.plot_lag()
|
|
||||||
|
spx.plot_frequency()
|
||||||
|
spx.plot_lag()
|
||||||
|
|
||||||
# sigma_x = 0.1
|
# sigma_x = 0.1
|
||||||
# sigma_v = 0.1 * np.ones(n_samples)
|
# sigma_v = 0.1 * np.ones(n_samples)
|
||||||
|
|
10
README.md
10
README.md
|
@ -1,14 +1,12 @@
|
||||||
# Signal Filtering
|
# Signal Smoothing / Filtering and Generation of Synthetic Time-Series
|
||||||
|
|
||||||
## Reference
|
## Reference
|
||||||
|
|
||||||
- John F. Ehlers, "[Cycle Analytics for Traders: Advanced Technical Trading Concepts](http://www.mesasoftware.com/ehlers_books.htm)."
|
- John F. Ehlers, "[Cycle Analytics for Traders: Advanced Technical Trading Concepts](http://www.mesasoftware.com/ehlers_books.htm)".
|
||||||
|
|
||||||
- Wikipedia, "[Alpha beta filter](https://en.wikipedia.org/wiki/Alpha_beta_filter)."
|
- D. Prichard and J. Theiler, "[Generating surrogate data for time series with several simultaneously measured variables](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.73.951)".
|
||||||
|
|
||||||
- D. Prichard, and J. Theiler, "[Generating surrogate data for time series with several simultaneously measured variables](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.73.951)."
|
- H. Vinod and J. Lopez-de-Lacalle, "[Maximum entropy bootstrap for time series: the meboot R package](https://www.jstatsoft.org/article/view/v029i05)".
|
||||||
|
|
||||||
- H. Vinod, and J. Lopez-de-Lacalle, "[Maximum entropy bootstrap for time series: the meboot R package](https://www.jstatsoft.org/article/view/v029i05)."
|
|
||||||
|
|
||||||
## Characteristics
|
## Characteristics
|
||||||
|
|
||||||
|
|
Ładowanie…
Reference in New Issue