diff --git a/Code_Python/filters.py b/Code_Python/filters.py index bff0ef0..d98e912 100644 --- a/Code_Python/filters.py +++ b/Code_Python/filters.py @@ -4,68 +4,61 @@ Signal Filtering/Smoothing and Generation of Synthetic Time-Series. Copyright (c) 2020 Gabriele Gilardi -X (n_samples, n_series) Dataset to filter -b (n_b, ) Numerator coefficients -a (n_a, ) Denominator coefficients -Y (n_samples, n_series) Filtered dataset -idx scalar First filtered element in Y +X (n_samples, ) Dataset to filter (input) +b (n_b, ) Transfer response coefficients (numerator) +a (n_a, ) Transfer response coefficients (denominator) +Y (n_samples, ) Filtered dataset (output) +idx scalar First filtered element in Y -n_samples Number of data to filter -n_series Number of series to filter -nb Number of coefficients (numerator) -na Number of coefficients (denominator) +n_samples Number of samples in the input dataset +nb Number of coefficients in array +na Number of coefficients in array Notes: -- the filter is applied starting from index. -- non filtered data are set equal to the original input, i.e. - Y[0:idx-1,:] = X[0:idx-1,:] -- if n_series = 1 then must be ( ..., 1) +- the filter is applied starting from index idx = MAX(0, nb-1, na-1). +- non filtered data are set equal to the input, i.e. Y[0:idx-1] = X[0:idx-1] +- X needs to be a 1D array. -Filters: -Generic b,a Generic case -SMA N Simple Moving Average -EMA N/alpha Exponential Moving Average +Filter list: +----------- +Generic b, Generic +SMA N Simple moving average +EMA N/alpha Exponential moving average WMA N Weighted moving average -MSMA N Modified Simple Moving Average -MLSQ N Modified Least-Squares Quadratic (N=5,7,9,11) -ButterOrig P,N Butterworth original (N=2,3) -ButterMod P,N Butterworth modified (N=2,3) -SuperSmooth P,N Super smoother (N=2,3) -GaussLow P,N Gauss low pass (P>=2) -GaussHigh P,N Gauss high pass (P>=5) -BandPass P,delta Band-pass filter -BandStop P,delta Band-stop filter -ZEMA1 N/alpha,K,Vn Zero-lag EMA (type 1) -ZEMA2 N/alpha,K Zero-lag EMA (type 2) +MSMA N Modified simple moving average +MLSQ N Modified least-squares quadratic (N = 5, 7, 9, 11) +ButterOrig P, N Butterworth original filter (N = 2, 3) +ButterMod P, N Butterworth modified filter (N = 2, 3) +SuperSmooth P, N Supersmoother filter (N = 2, 3) +GaussLow P, N Gauss low pass filter (P > 1) +GaussHigh P, N Gauss high pass filter (P > 4) +BandPass P, delta Band-pass filter +BandStop P, delta Band-stop filter +ZEMA1 N/alpha, K, Vn Zero-lag EMA (type 1) +ZEMA2 N/alpha, K Zero-lag EMA (type 2) InstTrend N/alpha Instantaneous trendline -SincFunction N Sinc function -Decycler P Decycler, 1-GaussHigh (P>=5) -DecyclerOsc P1,P2 Decycle oscillator, GH(P1) - GH(P2), (P1>=5) -ABG -Kalman +SincFilter P, nel Sinc function filter (N > 1) +Decycler P De-cycler filter (P >= 4) +DecyclerOsc P1, P2 De-cycle oscillator (P >= 4) +ABG alpha, beta, Alpha-beta-gamma filter (0 < alpha, beta < 1) + gamma, dt +Kalman sigma_x, dt One-dimensional steady-state Kalman filter + sigma_v N Order/smoothing factor/number of previous samples alpha Damping term P, P1, P2 Cut-off/critical period (50% power loss, -3 dB) -delta Band centered in P and in fraction - (30% of P => 0.3, = 0.3*P, if P = 12 => 0.3*12 = 4) +delta Semi-band centered in P K Coefficient/gain -Vn Look back sample (for the momentum) - -correction = update = measurement -prediction = motion - -X (n_states, 1) State estimate -P (n_states, n_states) Covariance estimate -F (n_states, n_states) State transition model -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) +Vn Look-back sample +nel Number of frequencies in the sinc function +alpha Parameter(s) to correct the position in the ABG filter +beta Parameter(s) to correct the velocity in the ABG filter +gamma Parameter(s) to correct the acceleration in the ABG filter +dt Sampling interval in the ABG and Kalman filters +sigma_x Process variance in the Kalman filter +sigma_v Noise variance in the Kalman filter """ import numpy as np @@ -73,19 +66,32 @@ from scipy import signal 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 with their names specified in + list . Each signal is plotted in its full length. """ - legend = [] - count = 0 + # Identify the signals by index if their name is not specified + if (names is None): + legend = [] + count = 0 + else: + legend = names + + # Loop over the signals for signal in signals: + signal = signal.flatten() end = len(signal) t = np.arange(start, end) plt.plot(t, signal[start:end]) - legend.append('Signal [' + str(count) + ']') - count += 1 + + # If no name is given use the list index to identify the signals + if (names is None): + legend.append('Signal [' + str(count) + ']') + count += 1 + + # Plot and format plt.xlabel('Index') plt.ylabel('Value') plt.grid(b=True) @@ -95,23 +101,27 @@ def plot_signals(signals, start=0): def filter_data(data, b, a): """ - Applies a filter with transfer response coefficients and . + Applies a filter with transfer response coefficients (numerator) and + (denominator). """ n_samples = len(data) nb = len(b) 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() + # Loop over the samples for i in range(idx, n_samples): tmp = 0 + # Contribution from the numerator term (input samples) 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): - tmp -= a[j] * Y[i-j] # Denominator term + tmp -= a[j] * Y[i-j] Y[i] = tmp / a[0] @@ -122,13 +132,16 @@ class Filter: def __init__(self, data): """ + Initialize the filter object. """ - self.data = np.asarray(data) - self.n_samples = len(data) + self.data = np.asarray(data).flatten() + self.idx = 0 + self.b = 0.0 + self.a = 0.0 def Generic(self, b=1.0, a=1.0): """ - Filter with generic transfer response coefficients and . + Filter with generic transfer response coefficients and . """ self.b = np.asarray(b) self.a = np.asarray(a) @@ -136,9 +149,9 @@ class Filter: 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.a = np.array([1.0]) @@ -146,9 +159,10 @@ class Filter: 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, is determined as equivalent to a N-SMA. """ if (alpha is None): @@ -160,10 +174,11 @@ class Filter: return Y - def WMA(self, N=10): + def WMA(self, N=5): """ - Weighted moving average (?? order, FIR, pass ??). - Example: N = 5 --> [5.0, 4.0, 3.0, 2.0, 1.0] / 15.0 + Weighted moving average. + + Example: N = 5 --> [5, 4, 3, 2, 1] / 15. """ w = np.arange(N, 0, -1) @@ -173,10 +188,11 @@ class Filter: return Y - def MSMA(self, N=10): + def MSMA(self, N=5): """ - Modified simple moving average (?? order, FIR, pass ??). - Example: N = 4 --> [0.5, 1.0, 1.0, 1.0, 0.5] / 4.0 + Modified simple moving average. + + Example: N = 5 --> [1/2, 1, 1, 1, 1, 1/2] / 5. """ w = np.ones(N+1) w[0] = 0.5 @@ -190,26 +206,26 @@ class Filter: def MLSQ(self, N=5): """ - Modified simple moving average (?? order, FIR, pass ??). - Only N = 5, 7, 9, and 11 are implemented. If not returns the unfiltered - dataset. + Modified least-squares quadratic. + + Must be N = 5, 7, 9, or 11. If wrong N, prints a warning and returns + the unfiltered dataset. """ 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): - 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): - w = np.array([-1.0, 28.0, 78.0, 108.0, 118.0, 108.0, 78.0, 28.0, - -1.0]) / 544.0 + w = np.array([-1, 28, 78, 108, 118, 108, 78, 28, -1]) / 544 elif (N == 11): - w = np.array([-11.0, 18.0, 88.0, 138.0, 168.0, 178.0, 168.0, 138.0, - 88.0, 18.0, -11.0]) / 980.0 + w = np.array([-11, 18, 88, 138, 168, 178, 168, 138, 88, 18, + -11]) / 980 else: - print("Warning: data returned unfiltered (wrong N)") + print("Warning: data returned unfiltered (MLSQ - Wrong N)") self.idx = 0 return self.data @@ -221,8 +237,10 @@ class Filter: def ButterOrig(self, N=2, P=10): """ - Butterworth original version (?? order, IIR, pass ??). - Only N = 2 and 3 are implemented. If not returns the unfiltered dataset. + Butterworth original filter. + + Must be N = 2 or 3. If wrong N, prints a warning and returns the + unfiltered dataset. """ if (N == 2): 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]) else: - print("Warning: data returned unfiltered (wrong N)") + print("Warning: data returned unfiltered (ButterOrig - Wrong N)") self.idx = 0 return self.data @@ -251,8 +269,11 @@ class Filter: def ButterMod(self, N=2, P=10): """ - Butterworth modified version (?? order, IIR, pass ??). - Only N = 2 and 3 are implemented. If not returns the unfiltered dataset. + Butterworth modified filter. It is derived from the Butterworth original + 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): 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]) else: - print("Warning: data returned unfiltered (wrong N)") + print("Warning: data returned unfiltered (ButterMod - Wrong N)") self.idx = 0 return self.data @@ -280,8 +301,11 @@ class Filter: def SuperSmooth(self, N=2, P=10): """ - SuperSmooth (?? order, IIR, pass ??). - Only N = 2 and 3 are implemented. If not returns the unfiltered dataset. + Supersmoother filter. It is derived from the Butterworth modified + 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): 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]) else: - print("Warning: data returned unfiltered (wrong N)") + print("Warning: data returned unfiltered (SuperSmooth - Wrong N)") self.idx = 0 return self.data @@ -308,13 +332,15 @@ class Filter: 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). - Must be P > 1. If not returns the unfiltered dataset. + Gauss low pass filter. + + Must be P > 1. If wrong P, prints a warning and returns the unfiltered + dataset. """ - if (P < 2): - print("Warning: data returned unfiltered (P < 2)") + if (P <= 1): + print("Warning: data returned unfiltered (GaussLow - Wrong P)") self.idx = 0 return self.data @@ -331,13 +357,15 @@ class Filter: return Y - def GaussHigh(self, N=1, P=5): + def GaussHigh(self, N=1, P=10): """ - Gauss high pass (IIR, Nth order, high pass). - Must be P > 4. If not returns the unfiltered dataset. + Gauss high pass filter. + + Must be P > 4. If wrong P, prints a warning and returns the unfiltered + dataset. """ - if (P < 5): - print("Warning: data returned unfiltered (P < 5)") + if (P <= 4): + print("Warning: data returned unfiltered (GaussHigh - Wrong P)") self.idx = 0 return self.data @@ -354,11 +382,11 @@ class Filter: return Y - def BandPass(self, P=5, delta=0.3): + def BandPass(self, P=10, delta=0.3): """ - Band-pass (type, order, IIR). - Example: delta = 0.3, P = 12 - (30% of P => 0.3, = 0.3*P, if P = 12 => 0.3*12 = 4) + Band-pass filter. + + Example: delta = 0.3, P = 10 --> 0.3 * 10 = 3 --> band is [7, 13] """ beta = np.cos(2.0 * np.pi / P) gamma = np.cos(4.0 * np.pi * delta / P) @@ -372,9 +400,9 @@ class Filter: def BandStop(self, P=5, delta=0.3): """ - Band-stop (type, order, IIR) - Example: delta = 0.3, P = 12 - (30% of P => 0.3, = 0.3*P, if P = 12 => 0.3*12 = 4) + Band-stop filter. + + Example: delta = 0.3, P = 10 --> 0.3 * 10 = 3 --> band is [7, 13] """ beta = np.cos(2.0 * np.pi / 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): """ - 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, is determined as equivalent to a N-SMA. """ if (alpha is None): alpha = 2.0 / (N + 1.0) - w = np.zeros(Vn+1) + w = np.zeros(Vn + 1) w[0] = alpha * (1.0 + K) w[Vn] = - alpha * K @@ -406,7 +436,9 @@ class Filter: 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, is determined as equivalent to a N-SMA. """ if (alpha is None): @@ -420,7 +452,9 @@ class Filter: 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, is determined as equivalent to a N-SMA. """ if (alpha is None): @@ -433,15 +467,22 @@ class Filter: return Y - def SincFunction(self, N=10, nel=10): + def SincFilter(self, P=10, nel=10): """ - Sinc function (order, FIR, pass). - (N > 1, cut off at 0.5/N) + Sinc function filter. The cut off point is at 0.5/P. + + 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) w = np.zeros(nel) 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.a = np.array([1.0]) @@ -451,12 +492,14 @@ class Filter: def Decycler(self, P=10): """ - Decycler (?? order, IIR, pass ??). Gauss,HP,1st,P - Built subtracting high pass Gauss filter from 1 (order 1) - Must be P > 4. If not returns the unfiltered dataset. + De-cycler filter. It is derived subtracting a 1st order high pass Gauss + filter from 1. + + Must be P > 4. If wrong P, prints a warning and returns the unfiltered + dataset. """ - if (P < 5): - print("Warning: data returned unfiltered (P < 5)") + if (P <= 4): + print("Warning: data returned unfiltered (Decycler - Wrong P)") self.idx = 0 return self.data @@ -465,17 +508,18 @@ class Filter: 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) - P1 = 1st cut off period, P2 = 2nd cut off period. Automatically fixed. - Must be P1, P2 > 4. If not returns the unfiltered dataset. + Must be P > 4. If wrong P, prints a warning and returns the unfiltered + dataset. """ P_low = np.amin([P1, P2]) P_high = np.amax([P1, P2]) - if (P_low < 5): - print("Warning: data returned unfiltered (P_low < 5)") + if (P_low <= 4): + print("Warning: data returned unfiltered (DecyclerOsc - Wrong P)") self.idx = 0 return self.data @@ -484,34 +528,38 @@ class Filter: def ABG(self, alpha=0.0, beta=0.0, gamma=0.0, dt=1.0): """ - alpha-beta-gamma - For numerical stability: 0 < alpha, beta < 1 + Alpha-beta-gamma filter. It is a predictor-corrector type of filter. + + 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): - alpha = np.ones(self.n_samples) * alpha + alpha = np.ones(n_samples) * alpha if (np.ndim(beta) == 0): - beta = np.ones(self.n_samples) * beta + beta = np.ones(n_samples) * beta if (np.ndim(gamma) == 0): - gamma = np.ones(self.n_samples) * gamma + gamma = np.ones(n_samples) * gamma # Initialize - Y_corr = self.data.copy() - Y_pred = self.data.copy() - x0 = self.data[0, :] - v0 = np.zeros(self.n_series) - a0 = np.zeros(self.n_series) + x0 = self.data[0] + v0 = 0.0 + a0 = 0.0 + Y[0] = x0 - for i in range(1, self.n_samples): + for i in range(1, n_samples): # Predictor (predicts state in ) x_pred = x0 + dt * v0 + 0.5 * a0 * dt ** 2.0 v_pred = v0 + dt * a0 a_pred = a0 - Y_pred[i, :] = x_pred # Residual (innovation) - r = self.data[i, :] - x_pred + r = self.data[i] - x_pred # Corrector (corrects state in ) 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 # Save value and prepare next iteration - Y_corr[i, :] = x_corr x0 = x_corr v0 = v_corr a0 = a_corr + Y[i] = x_corr self.idx = 1 - return Y_corr, Y_pred + return Y 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 @@ -561,19 +614,23 @@ class Filter: beta = 2.0 * (1 - s) ** 2.0 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) return Y def plot_frequency(self): """ + Plots the frequency response (in decibels) of the filter with transfer + response coefficients and . """ w, h = signal.freqz(self.b, self.a) - h_db = 20.0 * np.log10(np.abs(h)) - wf = w / (2.0 * np.pi) + h_db = 20.0 * np.log10(np.abs(h)) # Convert to decibels + wf = w / (2.0 * np.pi) # Scale to [0, 0.5] + + # Plot and format 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.xlim(np.amin(wf), np.amax(wf)) plt.xlabel(r'$\omega$ [rad/sample]') @@ -584,9 +641,13 @@ class Filter: def plot_lag(self): """ + Plots the lag (group delay) of the filter with transfer response + coefficients and . """ 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.grid(b=True) plt.xlim(np.amin(wf), np.amax(wf)) diff --git a/Code_Python/synthetic.py b/Code_Python/synthetic.py index 2401e0f..6079f08 100644 --- a/Code_Python/synthetic.py +++ b/Code_Python/synthetic.py @@ -49,14 +49,12 @@ def synthetic_wave(P, A=None, phi=None, num=1000): def synthetic_FFT(X, n_reps=1): """ 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_fft (n, ) FFT of the original time-series X_synt_fft (n_reps, n) FFT of the 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, ) n = len(X) @@ -82,7 +80,7 @@ def synthetic_FFT(X, n_reps=1): # FFT of the synthetic time-series (1st sample is unchanged) X_synt_fft = np.zeros((n_reps, n), dtype=complex) - X_synt_fft[:, 0] = X_fft[0] + X_synt_fft[:, 0] = X_fft[0] X_synt_fft[:, idx1] = X_fft[idx1] * phases1 # 1st half X_synt_fft[:, idx2] = X_fft[idx2] * phases2 # 2nd half @@ -95,13 +93,11 @@ def synthetic_FFT(X, n_reps=1): def synthetic_sampling(X, n_reps=1, replace=True): """ 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 idx (n_reps, n) Random index of X 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, ) 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): """ 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 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 and last interval 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, ) n = len(X) @@ -259,16 +253,16 @@ def value2diff(X, percent=True): """ Returns the 1st discrete difference of array X. - X (n, ) Original dataset + X (n, ) Input dataset dX (n-1, ) 1st discrete differences - + Notes: - the discrete difference can be calculated in percent or in value. - - array dX is one element shorter than array X. - - array X needs to be a 1D array (of any shape). + - dX is one element shorter than X. + - X needs to be a 1D array. """ X = X.flatten() # Reshape to (n, ) - + # Discrete difference in percent if (percent): dX = X[1:] / X[:-1] - 1.0 @@ -282,36 +276,39 @@ def value2diff(X, 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 X0 scalar Initial value - X (n+1, ) Original dataset ????? - + X (n+1, ) Output dataset + Notes: - the discrete difference can be in percent or in value. - - array X is one element longer than array dX. - - array dX needs to be a 1D array (of any shape). + - X is one element longer than dX. + - dX needs to be a 1D array. - If the discrete difference is in percent, the first column of X is set to - one. The original array is X[0] * X + If the discrete difference is in percent: + 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 - zero. X0+X - - array X needs to be a 1D array (of any shape). + If the discrete difference is in value: + X[0] = X0 + X[1] = X[0] + dX[0] + X[2] = X[1] + dX[1] = X[0] + dX[0] + dX[1] + .... """ - dX = dX.flatten() # Reshape to (n, ) - n = len(dX) - X = np.zeros(n+1) + dX = dX.flatten() # Reshape to (n, ) + X = np.zeros(len(dX) + 1) + X[0] = X0 # Initial value # Discrete difference in percent if (percent): - X[0] = X0 - X[1:] = np.cumprod(1.0 + dX) * X0 # !!!! check + X[1:] = X0 * np.cumprod(1.0 + dX) # Discrete difference in value else: - X[0] = X0 - X[1:] = np.cumsum(dX) + X0 # !!!!!c check + X[1:] = X0 + np.cumsum(dX) return X diff --git a/Code_Python/test.py b/Code_Python/test.py index ed32bc1..5c6b82f 100644 --- a/Code_Python/test.py +++ b/Code_Python/test.py @@ -5,10 +5,10 @@ Copyright (c) 2020 Gabriele Gilardi ToDo: -- add comments to the code -- in comments write what filters do -- is necessary to copy X for Y untouched? -- decide default values in functions +- do examples and check type (low pass, etc) +- put type in description +- see paper ZEMA for tests +- check results with ML """ import sys @@ -31,18 +31,21 @@ data_file = sys.argv[1] + '.csv' # Read data from a csv file (one time-series each column) data = np.loadtxt(data_file, delimiter=',') +print(data.shape) -t, f = syn.synthetic_wave([1., 2., 3.], A=None, phi=None, num=1000) -plt.plot(t,f) -plt.show() +# t, f = syn.synthetic_wave([1., 2., 3.], A=None, phi=None, num=1000) +# plt.plot(t,f) +# plt.show() -# spx = flt.Filter(data) +spx = flt.Filter(data) # res = spx.EMA(N=10) # signals = [spx.data, res[0:400]] -# flt.plot_signals(signals, start=100) +# flt.plot_signals(signals, ['SPX', 'SMA']) -# spx.plot_frequency() -# spx.plot_lag() +res = spx.BandPass(P=10, delta=0.3) + +spx.plot_frequency() +spx.plot_lag() # sigma_x = 0.1 # sigma_v = 0.1 * np.ones(n_samples) diff --git a/README.md b/README.md index ecbbc45..17f6fb8 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,12 @@ -# Signal Filtering +# Signal Smoothing / Filtering and Generation of Synthetic Time-Series ## 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