kopia lustrzana https://github.com/gabrielegilardi/SignalFilters
add comments
rodzic
afad3296ed
commit
753245f285
|
@ -1,14 +1,9 @@
|
||||||
"""
|
"""
|
||||||
Class for filter/smooth data.
|
Signal Filtering/Smoothing and Generation of Synthetic Time-Series.
|
||||||
|
|
||||||
Copyright (c) 2020 Gabriele Gilardi
|
Copyright (c) 2020 Gabriele Gilardi
|
||||||
|
|
||||||
|
|
||||||
References (both from John F. Ehlers):
|
|
||||||
[1] "Cycle Analytics for Traders: Advanced Technical Trading Concepts".
|
|
||||||
[2] "Signal Analysis, Filters And Trading Strategies".
|
|
||||||
|
|
||||||
|
|
||||||
X (n_samples, n_series) Dataset to filter
|
X (n_samples, n_series) Dataset to filter
|
||||||
b (n_b, ) Numerator coefficients
|
b (n_b, ) Numerator coefficients
|
||||||
a (n_a, ) Denominator coefficients
|
a (n_a, ) Denominator coefficients
|
||||||
|
@ -75,23 +70,26 @@ Y (n_obs, 1) Observation residual (innovation)
|
||||||
|
|
||||||
import numpy as np
|
import numpy as np
|
||||||
from scipy import signal
|
from scipy import signal
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
|
||||||
|
|
||||||
def plot_signals(signals, idx_start=0, idx_end=None):
|
def plot_signals(signals, start=0):
|
||||||
"""
|
"""
|
||||||
signals must be all same lenght
|
signals must be a list
|
||||||
"""
|
"""
|
||||||
if (idx_end is None):
|
legend = []
|
||||||
idx_end = len(signals[0])
|
|
||||||
t = np.arange(idx_start, idx_end)
|
|
||||||
names = []
|
|
||||||
count = 0
|
count = 0
|
||||||
for signal in signals:
|
for signal in signals:
|
||||||
plt.plot(t, signal[idx_start:idx_end])
|
signal = signal.flatten()
|
||||||
names.append(str(count))
|
end = len(signal)
|
||||||
|
t = np.arange(start, end)
|
||||||
|
plt.plot(t, signal[start:end])
|
||||||
|
legend.append('Signal [' + str(count) + ']')
|
||||||
count += 1
|
count += 1
|
||||||
|
plt.xlabel('Index')
|
||||||
|
plt.ylabel('Value')
|
||||||
plt.grid(b=True)
|
plt.grid(b=True)
|
||||||
plt.legend(names)
|
plt.legend(legend)
|
||||||
plt.show()
|
plt.show()
|
||||||
|
|
||||||
|
|
||||||
|
@ -99,22 +97,23 @@ def filter_data(data, b, a):
|
||||||
"""
|
"""
|
||||||
Applies a filter with transfer response coefficients <a> and <b>.
|
Applies a filter with transfer response coefficients <a> and <b>.
|
||||||
"""
|
"""
|
||||||
n_samples, n_series = data.shape
|
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])
|
||||||
Y = data.copy()
|
Y = data.copy()
|
||||||
|
|
||||||
for i in range(idx, n_samples):
|
for i in range(idx, n_samples):
|
||||||
tmp = np.zeros(n_series)
|
|
||||||
|
tmp = 0
|
||||||
|
|
||||||
for j in range(nb):
|
for j in range(nb):
|
||||||
tmp = tmp + b[j] * data[i-j, :] # Numerator term
|
tmp += b[j] * data[i-j] # Numerator term
|
||||||
|
|
||||||
for j in range(1, na):
|
for j in range(1, na):
|
||||||
tmp = tmp - a[j] * Y[i-j, :] # Denominator term
|
tmp -= a[j] * Y[i-j] # Denominator term
|
||||||
|
|
||||||
Y[i, :] = tmp / a[0]
|
Y[i] = tmp / a[0]
|
||||||
|
|
||||||
return Y, idx
|
return Y, idx
|
||||||
|
|
||||||
|
@ -125,7 +124,7 @@ class Filter:
|
||||||
"""
|
"""
|
||||||
"""
|
"""
|
||||||
self.data = np.asarray(data)
|
self.data = np.asarray(data)
|
||||||
self.n_samples, self.n_series = data.shape
|
self.n_samples = len(data)
|
||||||
|
|
||||||
def Generic(self, b=1.0, a=1.0):
|
def Generic(self, b=1.0, a=1.0):
|
||||||
"""
|
"""
|
||||||
|
@ -577,8 +576,10 @@ class Filter:
|
||||||
plt.axhline(-3.0, lw=1.5, ls='--', C='r')
|
plt.axhline(-3.0, lw=1.5, ls='--', C='r')
|
||||||
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('$omega$ [rad/sample]')
|
plt.xlabel(r'$\omega$ [rad/sample]')
|
||||||
plt.ylabel('$h$ [db]')
|
plt.ylabel('$h$ [db]')
|
||||||
|
legend = ['Filter', '-3dB']
|
||||||
|
plt.legend(legend)
|
||||||
plt.show()
|
plt.show()
|
||||||
|
|
||||||
def plot_lag(self):
|
def plot_lag(self):
|
||||||
|
@ -589,6 +590,6 @@ class Filter:
|
||||||
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))
|
||||||
plt.xlabel('$omega$ [rad/sample]')
|
plt.xlabel(r'$\omega$ [rad/sample]')
|
||||||
plt.ylabel('$gd$ [samples]')
|
plt.ylabel('$gd$ [samples]')
|
||||||
plt.show()
|
plt.show()
|
||||||
|
|
|
@ -1,500 +0,0 @@
|
||||||
1399.48,1399.48,1399.48
|
|
||||||
1410.49,1410.49,1410.49
|
|
||||||
1409.3,1409.3,1409.3
|
|
||||||
1410.44,1410.44,1410.44
|
|
||||||
1411.13,1411.13,1411.13
|
|
||||||
1402.08,1402.08,1402.08
|
|
||||||
1413.49,1413.49,1413.49
|
|
||||||
1413.17,1413.17,1413.17
|
|
||||||
1418.13,1418.13,1418.13
|
|
||||||
1418.16,1418.16,1418.16
|
|
||||||
1415.51,1415.51,1415.51
|
|
||||||
1405.53,1405.53,1405.53
|
|
||||||
1403.93,1403.93,1403.93
|
|
||||||
1404.11,1404.11,1404.11
|
|
||||||
1405.87,1405.87,1405.87
|
|
||||||
1402.8,1402.8,1402.8
|
|
||||||
1402.22,1402.22,1402.22
|
|
||||||
1401.35,1401.35,1401.35
|
|
||||||
1394.23,1394.23,1394.23
|
|
||||||
1390.99,1390.99,1390.99
|
|
||||||
1365,1365,1365
|
|
||||||
1375.32,1375.32,1375.32
|
|
||||||
1379.32,1379.32,1379.32
|
|
||||||
1385.3,1385.3,1385.3
|
|
||||||
1385.97,1385.97,1385.97
|
|
||||||
1360.02,1360.02,1360.02
|
|
||||||
1337.89,1337.89,1337.89
|
|
||||||
1338.31,1338.31,1338.31
|
|
||||||
1350.52,1350.52,1350.52
|
|
||||||
1362.66,1362.66,1362.66
|
|
||||||
1376.51,1376.51,1376.51
|
|
||||||
1372.78,1372.78,1372.78
|
|
||||||
1363.67,1363.67,1363.67
|
|
||||||
1353.64,1353.64,1353.64
|
|
||||||
1356.78,1356.78,1356.78
|
|
||||||
1334.76,1334.76,1334.76
|
|
||||||
1341.45,1341.45,1341.45
|
|
||||||
1341.47,1341.47,1341.47
|
|
||||||
1352.46,1352.46,1352.46
|
|
||||||
1354.68,1354.68,1354.68
|
|
||||||
1367.58,1367.58,1367.58
|
|
||||||
1374.02,1374.02,1374.02
|
|
||||||
1365.51,1365.51,1365.51
|
|
||||||
1362.16,1362.16,1362.16
|
|
||||||
1329.04,1329.04,1329.04
|
|
||||||
1331.85,1331.85,1331.85
|
|
||||||
1319.99,1319.99,1319.99
|
|
||||||
1313.72,1313.72,1313.72
|
|
||||||
1335.02,1335.02,1335.02
|
|
||||||
1325.51,1325.51,1325.51
|
|
||||||
1355.69,1355.69,1355.69
|
|
||||||
1357.98,1357.98,1357.98
|
|
||||||
1344.78,1344.78,1344.78
|
|
||||||
1342.84,1342.84,1342.84
|
|
||||||
1329.1,1329.1,1329.1
|
|
||||||
1314.88,1314.88,1314.88
|
|
||||||
1324.18,1324.18,1324.18
|
|
||||||
1308.93,1308.93,1308.93
|
|
||||||
1325.66,1325.66,1325.66
|
|
||||||
1314.99,1314.99,1314.99
|
|
||||||
1315.13,1315.13,1315.13
|
|
||||||
1285.5,1285.5,1285.5
|
|
||||||
1278.18,1278.18,1278.18
|
|
||||||
1278.04,1278.04,1278.04
|
|
||||||
1310.33,1310.33,1310.33
|
|
||||||
1313.32,1313.32,1313.32
|
|
||||||
1332.42,1332.42,1332.42
|
|
||||||
1317.82,1317.82,1317.82
|
|
||||||
1320.68,1320.68,1320.68
|
|
||||||
1318.86,1318.86,1318.86
|
|
||||||
1316.63,1316.63,1316.63
|
|
||||||
1315.99,1315.99,1315.99
|
|
||||||
1295.22,1295.22,1295.22
|
|
||||||
1304.86,1304.86,1304.86
|
|
||||||
1324.8,1324.8,1324.8
|
|
||||||
1330.66,1330.66,1330.66
|
|
||||||
1338.35,1338.35,1338.35
|
|
||||||
1353.39,1353.39,1353.39
|
|
||||||
1357.99,1357.99,1357.99
|
|
||||||
1354.58,1354.58,1354.58
|
|
||||||
1363.72,1363.72,1363.72
|
|
||||||
1369.58,1369.58,1369.58
|
|
||||||
1369.1,1369.1,1369.1
|
|
||||||
1391.57,1391.57,1391.57
|
|
||||||
1402.31,1402.31,1402.31
|
|
||||||
1405.82,1405.82,1405.82
|
|
||||||
1397.91,1397.91,1397.91
|
|
||||||
1403.36,1403.36,1403.36
|
|
||||||
1399.98,1399.98,1399.98
|
|
||||||
1390.69,1390.69,1390.69
|
|
||||||
1371.97,1371.97,1371.97
|
|
||||||
1366.94,1366.94,1366.94
|
|
||||||
1378.53,1378.53,1378.53
|
|
||||||
1376.92,1376.92,1376.92
|
|
||||||
1385.14,1385.14,1385.14
|
|
||||||
1390.78,1390.78,1390.78
|
|
||||||
1369.57,1369.57,1369.57
|
|
||||||
1370.26,1370.26,1370.26
|
|
||||||
1387.57,1387.57,1387.57
|
|
||||||
1368.71,1368.71,1368.71
|
|
||||||
1358.59,1358.59,1358.59
|
|
||||||
1382.2,1382.2,1382.2
|
|
||||||
1398.08,1398.08,1398.08
|
|
||||||
1398.96,1398.96,1398.96
|
|
||||||
1413.38,1413.38,1413.38
|
|
||||||
1419.04,1419.04,1419.04
|
|
||||||
1408.47,1408.47,1408.47
|
|
||||||
1403.28,1403.28,1403.28
|
|
||||||
1405.54,1405.54,1405.54
|
|
||||||
1412.52,1412.52,1412.52
|
|
||||||
1416.51,1416.51,1416.51
|
|
||||||
1397.11,1397.11,1397.11
|
|
||||||
1392.78,1392.78,1392.78
|
|
||||||
1402.89,1402.89,1402.89
|
|
||||||
1405.52,1405.52,1405.52
|
|
||||||
1409.75,1409.75,1409.75
|
|
||||||
1404.17,1404.17,1404.17
|
|
||||||
1402.6,1402.6,1402.6
|
|
||||||
1394.28,1394.28,1394.28
|
|
||||||
1395.95,1395.95,1395.95
|
|
||||||
1371.09,1371.09,1371.09
|
|
||||||
1370.87,1370.87,1370.87
|
|
||||||
1365.91,1365.91,1365.91
|
|
||||||
1352.63,1352.63,1352.63
|
|
||||||
1343.36,1343.36,1343.36
|
|
||||||
1364.33,1364.33,1364.33
|
|
||||||
1369.63,1369.63,1369.63
|
|
||||||
1374.09,1374.09,1374.09
|
|
||||||
1365.68,1365.68,1365.68
|
|
||||||
1372.18,1372.18,1372.18
|
|
||||||
1367.59,1367.59,1367.59
|
|
||||||
1365.74,1365.74,1365.74
|
|
||||||
1363.46,1363.46,1363.46
|
|
||||||
1357.66,1357.66,1357.66
|
|
||||||
1362.21,1362.21,1362.21
|
|
||||||
1361.23,1361.23,1361.23
|
|
||||||
1358.04,1358.04,1358.04
|
|
||||||
1343.23,1343.23,1343.23
|
|
||||||
1350.5,1350.5,1350.5
|
|
||||||
1351.77,1351.77,1351.77
|
|
||||||
1342.64,1342.64,1342.64
|
|
||||||
1351.95,1351.95,1351.95
|
|
||||||
1349.96,1349.96,1349.96
|
|
||||||
1347.05,1347.05,1347.05
|
|
||||||
1344.33,1344.33,1344.33
|
|
||||||
1344.9,1344.9,1344.9
|
|
||||||
1325.54,1325.54,1325.54
|
|
||||||
1324.09,1324.09,1324.09
|
|
||||||
1312.41,1312.41,1312.41
|
|
||||||
1313.01,1313.01,1313.01
|
|
||||||
1316.33,1316.33,1316.33
|
|
||||||
1318.43,1318.43,1318.43
|
|
||||||
1326.06,1326.06,1326.06
|
|
||||||
1314.65,1314.65,1314.65
|
|
||||||
1316,1316,1316
|
|
||||||
1315.38,1315.38,1315.38
|
|
||||||
1314.5,1314.5,1314.5
|
|
||||||
1308.04,1308.04,1308.04
|
|
||||||
1293.67,1293.67,1293.67
|
|
||||||
1289.09,1289.09,1289.09
|
|
||||||
1295.5,1295.5,1295.5
|
|
||||||
1292.48,1292.48,1292.48
|
|
||||||
1292.08,1292.08,1292.08
|
|
||||||
1280.7,1280.7,1280.7
|
|
||||||
1277.81,1277.81,1277.81
|
|
||||||
1281.06,1281.06,1281.06
|
|
||||||
1277.3,1277.3,1277.3
|
|
||||||
1277.06,1277.06,1277.06
|
|
||||||
1257.6,1257.6,1257.6
|
|
||||||
1263.02,1263.02,1263.02
|
|
||||||
1249.64,1249.64,1249.64
|
|
||||||
1265.43,1265.43,1265.43
|
|
||||||
1265.33,1265.33,1265.33
|
|
||||||
1254,1254,1254
|
|
||||||
1243.72,1243.72,1243.72
|
|
||||||
1241.3,1241.3,1241.3
|
|
||||||
1205.35,1205.35,1205.35
|
|
||||||
1219.66,1219.66,1219.66
|
|
||||||
1215.75,1215.75,1215.75
|
|
||||||
1211.82,1211.82,1211.82
|
|
||||||
1225.73,1225.73,1225.73
|
|
||||||
1236.47,1236.47,1236.47
|
|
||||||
1255.19,1255.19,1255.19
|
|
||||||
1234.35,1234.35,1234.35
|
|
||||||
1261.01,1261.01,1261.01
|
|
||||||
1258.47,1258.47,1258.47
|
|
||||||
1257.08,1257.08,1257.08
|
|
||||||
1244.28,1244.28,1244.28
|
|
||||||
1244.58,1244.58,1244.58
|
|
||||||
1246.96,1246.96,1246.96
|
|
||||||
1195.19,1195.19,1195.19
|
|
||||||
1192.55,1192.55,1192.55
|
|
||||||
1158.67,1158.67,1158.67
|
|
||||||
1161.79,1161.79,1161.79
|
|
||||||
1188.04,1188.04,1188.04
|
|
||||||
1192.98,1192.98,1192.98
|
|
||||||
1215.65,1215.65,1215.65
|
|
||||||
1216.13,1216.13,1216.13
|
|
||||||
1236.91,1236.91,1236.91
|
|
||||||
1257.81,1257.81,1257.81
|
|
||||||
1251.78,1251.78,1251.78
|
|
||||||
1263.85,1263.85,1263.85
|
|
||||||
1239.7,1239.7,1239.7
|
|
||||||
1229.1,1229.1,1229.1
|
|
||||||
1275.92,1275.92,1275.92
|
|
||||||
1261.12,1261.12,1261.12
|
|
||||||
1253.23,1253.23,1253.23
|
|
||||||
1261.15,1261.15,1261.15
|
|
||||||
1237.9,1237.9,1237.9
|
|
||||||
1218.28,1218.28,1218.28
|
|
||||||
1253.3,1253.3,1253.3
|
|
||||||
1285.09,1285.09,1285.09
|
|
||||||
1284.59,1284.59,1284.59
|
|
||||||
1242,1242,1242
|
|
||||||
1229.05,1229.05,1229.05
|
|
||||||
1254.19,1254.19,1254.19
|
|
||||||
1238.25,1238.25,1238.25
|
|
||||||
1215.39,1215.39,1215.39
|
|
||||||
1209.88,1209.88,1209.88
|
|
||||||
1225.38,1225.38,1225.38
|
|
||||||
1200.86,1200.86,1200.86
|
|
||||||
1224.58,1224.58,1224.58
|
|
||||||
1203.66,1203.66,1203.66
|
|
||||||
1207.25,1207.25,1207.25
|
|
||||||
1195.54,1195.54,1195.54
|
|
||||||
1194.89,1194.89,1194.89
|
|
||||||
1155.46,1155.46,1155.46
|
|
||||||
1164.97,1164.97,1164.97
|
|
||||||
1144.03,1144.03,1144.03
|
|
||||||
1123.95,1123.95,1123.95
|
|
||||||
1099.23,1099.23,1099.23
|
|
||||||
1131.42,1131.42,1131.42
|
|
||||||
1160.4,1160.4,1160.4
|
|
||||||
1151.06,1151.06,1151.06
|
|
||||||
1175.38,1175.38,1175.38
|
|
||||||
1162.95,1162.95,1162.95
|
|
||||||
1136.43,1136.43,1136.43
|
|
||||||
1129.56,1129.56,1129.56
|
|
||||||
1166.76,1166.76,1166.76
|
|
||||||
1202.09,1202.09,1202.09
|
|
||||||
1204.09,1204.09,1204.09
|
|
||||||
1216.01,1216.01,1216.01
|
|
||||||
1209.11,1209.11,1209.11
|
|
||||||
1188.68,1188.68,1188.68
|
|
||||||
1172.87,1172.87,1172.87
|
|
||||||
1162.27,1162.27,1162.27
|
|
||||||
1154.23,1154.23,1154.23
|
|
||||||
1185.9,1185.9,1185.9
|
|
||||||
1198.62,1198.62,1198.62
|
|
||||||
1165.24,1165.24,1165.24
|
|
||||||
1173.97,1173.97,1173.97
|
|
||||||
1204.42,1204.42,1204.42
|
|
||||||
1218.89,1218.89,1218.89
|
|
||||||
1212.92,1212.92,1212.92
|
|
||||||
1210.08,1210.08,1210.08
|
|
||||||
1176.8,1176.8,1176.8
|
|
||||||
1159.27,1159.27,1159.27
|
|
||||||
1177.6,1177.6,1177.6
|
|
||||||
1162.35,1162.35,1162.35
|
|
||||||
1123.82,1123.82,1123.82
|
|
||||||
1123.53,1123.53,1123.53
|
|
||||||
1140.65,1140.65,1140.65
|
|
||||||
1193.89,1193.89,1193.89
|
|
||||||
1192.76,1192.76,1192.76
|
|
||||||
1204.49,1204.49,1204.49
|
|
||||||
1178.81,1178.81,1178.81
|
|
||||||
1172.64,1172.64,1172.64
|
|
||||||
1120.76,1120.76,1120.76
|
|
||||||
1172.53,1172.53,1172.53
|
|
||||||
1119.46,1119.46,1119.46
|
|
||||||
1199.38,1199.38,1199.38
|
|
||||||
1200.07,1200.07,1200.07
|
|
||||||
1260.34,1260.34,1260.34
|
|
||||||
1254.05,1254.05,1254.05
|
|
||||||
1286.94,1286.94,1286.94
|
|
||||||
1292.28,1292.28,1292.28
|
|
||||||
1300.67,1300.67,1300.67
|
|
||||||
1304.89,1304.89,1304.89
|
|
||||||
1331.94,1331.94,1331.94
|
|
||||||
1337.43,1337.43,1337.43
|
|
||||||
1345.02,1345.02,1345.02
|
|
||||||
1343.8,1343.8,1343.8
|
|
||||||
1325.84,1325.84,1325.84
|
|
||||||
1326.73,1326.73,1326.73
|
|
||||||
1305.44,1305.44,1305.44
|
|
||||||
1316.14,1316.14,1316.14
|
|
||||||
1308.87,1308.87,1308.87
|
|
||||||
1317.72,1317.72,1317.72
|
|
||||||
1313.64,1313.64,1313.64
|
|
||||||
1319.49,1319.49,1319.49
|
|
||||||
1343.8,1343.8,1343.8
|
|
||||||
1353.22,1353.22,1353.22
|
|
||||||
1339.22,1339.22,1339.22
|
|
||||||
1337.88,1337.88,1337.88
|
|
||||||
1339.67,1339.67,1339.67
|
|
||||||
1320.64,1320.64,1320.64
|
|
||||||
1307.41,1307.41,1307.41
|
|
||||||
1296.67,1296.67,1296.67
|
|
||||||
1280.1,1280.1,1280.1
|
|
||||||
1268.45,1268.45,1268.45
|
|
||||||
1283.5,1283.5,1283.5
|
|
||||||
1287.14,1287.14,1287.14
|
|
||||||
1295.52,1295.52,1295.52
|
|
||||||
1278.36,1278.36,1278.36
|
|
||||||
1271.5,1271.5,1271.5
|
|
||||||
1267.64,1267.64,1267.64
|
|
||||||
1265.42,1265.42,1265.42
|
|
||||||
1287.87,1287.87,1287.87
|
|
||||||
1271.83,1271.83,1271.83
|
|
||||||
1270.98,1270.98,1270.98
|
|
||||||
1289,1289,1289
|
|
||||||
1279.56,1279.56,1279.56
|
|
||||||
1284.94,1284.94,1284.94
|
|
||||||
1286.17,1286.17,1286.17
|
|
||||||
1300.16,1300.16,1300.16
|
|
||||||
1312.94,1312.94,1312.94
|
|
||||||
1314.55,1314.55,1314.55
|
|
||||||
1345.2,1345.2,1345.2
|
|
||||||
1331.1,1331.1,1331.1
|
|
||||||
1325.69,1325.69,1325.69
|
|
||||||
1320.47,1320.47,1320.47
|
|
||||||
1316.28,1316.28,1316.28
|
|
||||||
1317.37,1317.37,1317.37
|
|
||||||
1333.27,1333.27,1333.27
|
|
||||||
1343.6,1343.6,1343.6
|
|
||||||
1340.68,1340.68,1340.68
|
|
||||||
1328.98,1328.98,1328.98
|
|
||||||
1329.47,1329.47,1329.47
|
|
||||||
1337.77,1337.77,1337.77
|
|
||||||
1348.65,1348.65,1348.65
|
|
||||||
1342.08,1342.08,1342.08
|
|
||||||
1357.16,1357.16,1357.16
|
|
||||||
1346.29,1346.29,1346.29
|
|
||||||
1340.2,1340.2,1340.2
|
|
||||||
1335.1,1335.1,1335.1
|
|
||||||
1347.32,1347.32,1347.32
|
|
||||||
1356.62,1356.62,1356.62
|
|
||||||
1361.22,1361.22,1361.22
|
|
||||||
1363.61,1363.61,1363.61
|
|
||||||
1360.48,1360.48,1360.48
|
|
||||||
1355.66,1355.66,1355.66
|
|
||||||
1347.24,1347.24,1347.24
|
|
||||||
1335.25,1335.25,1335.25
|
|
||||||
1337.38,1337.38,1337.38
|
|
||||||
1330.36,1330.36,1330.36
|
|
||||||
1312.62,1312.62,1312.62
|
|
||||||
1305.14,1305.14,1305.14
|
|
||||||
1319.68,1319.68,1319.68
|
|
||||||
1314.52,1314.52,1314.52
|
|
||||||
1314.41,1314.41,1314.41
|
|
||||||
1314.16,1314.16,1314.16
|
|
||||||
1324.46,1324.46,1324.46
|
|
||||||
1328.17,1328.17,1328.17
|
|
||||||
1333.51,1333.51,1333.51
|
|
||||||
1335.54,1335.54,1335.54
|
|
||||||
1332.63,1332.63,1332.63
|
|
||||||
1332.87,1332.87,1332.87
|
|
||||||
1332.41,1332.41,1332.41
|
|
||||||
1325.83,1325.83,1325.83
|
|
||||||
1328.26,1328.26,1328.26
|
|
||||||
1319.44,1319.44,1319.44
|
|
||||||
1310.19,1310.19,1310.19
|
|
||||||
1313.8,1313.8,1313.8
|
|
||||||
1309.66,1309.66,1309.66
|
|
||||||
1297.54,1297.54,1297.54
|
|
||||||
1293.77,1293.77,1293.77
|
|
||||||
1298.38,1298.38,1298.38
|
|
||||||
1279.2,1279.2,1279.2
|
|
||||||
1273.72,1273.72,1273.72
|
|
||||||
1256.88,1256.88,1256.88
|
|
||||||
1281.87,1281.87,1281.87
|
|
||||||
1296.39,1296.39,1296.39
|
|
||||||
1304.28,1304.28,1304.28
|
|
||||||
1295.11,1295.11,1295.11
|
|
||||||
1320.02,1320.02,1320.02
|
|
||||||
1321.82,1321.82,1321.82
|
|
||||||
1310.13,1310.13,1310.13
|
|
||||||
1321.15,1321.15,1321.15
|
|
||||||
1330.97,1330.97,1330.97
|
|
||||||
1308.44,1308.44,1308.44
|
|
||||||
1306.33,1306.33,1306.33
|
|
||||||
1327.22,1327.22,1327.22
|
|
||||||
1319.88,1319.88,1319.88
|
|
||||||
1306.1,1306.1,1306.1
|
|
||||||
1307.4,1307.4,1307.4
|
|
||||||
1315.44,1315.44,1315.44
|
|
||||||
1343.01,1343.01,1343.01
|
|
||||||
1340.43,1340.43,1340.43
|
|
||||||
1336.32,1336.32,1336.32
|
|
||||||
1328.01,1328.01,1328.01
|
|
||||||
1332.32,1332.32,1332.32
|
|
||||||
1329.15,1329.15,1329.15
|
|
||||||
1321.87,1321.87,1321.87
|
|
||||||
1320.88,1320.88,1320.88
|
|
||||||
1324.57,1324.57,1324.57
|
|
||||||
1319.05,1319.05,1319.05
|
|
||||||
1310.87,1310.87,1310.87
|
|
||||||
1307.1,1307.1,1307.1
|
|
||||||
1304.03,1304.03,1304.03
|
|
||||||
1307.59,1307.59,1307.59
|
|
||||||
1286.12,1286.12,1286.12
|
|
||||||
1276.34,1276.34,1276.34
|
|
||||||
1299.54,1299.54,1299.54
|
|
||||||
1296.63,1296.63,1296.63
|
|
||||||
1291.18,1291.18,1291.18
|
|
||||||
1290.84,1290.84,1290.84
|
|
||||||
1283.35,1283.35,1283.35
|
|
||||||
1280.26,1280.26,1280.26
|
|
||||||
1281.92,1281.92,1281.92
|
|
||||||
1295.02,1295.02,1295.02
|
|
||||||
1293.24,1293.24,1293.24
|
|
||||||
1283.76,1283.76,1283.76
|
|
||||||
1285.96,1285.96,1285.96
|
|
||||||
1274.48,1274.48,1274.48
|
|
||||||
1269.75,1269.75,1269.75
|
|
||||||
1271.5,1271.5,1271.5
|
|
||||||
1273.85,1273.85,1273.85
|
|
||||||
1276.56,1276.56,1276.56
|
|
||||||
1270.2,1270.2,1270.2
|
|
||||||
1271.89,1271.89,1271.89
|
|
||||||
1257.64,1257.64,1257.64
|
|
||||||
1257.88,1257.88,1257.88
|
|
||||||
1259.78,1259.78,1259.78
|
|
||||||
1258.51,1258.51,1258.51
|
|
||||||
1257.54,1257.54,1257.54
|
|
||||||
1256.77,1256.77,1256.77
|
|
||||||
1258.84,1258.84,1258.84
|
|
||||||
1254.6,1254.6,1254.6
|
|
||||||
1247.08,1247.08,1247.08
|
|
||||||
1243.91,1243.91,1243.91
|
|
||||||
1242.87,1242.87,1242.87
|
|
||||||
1235.23,1235.23,1235.23
|
|
||||||
1241.59,1241.59,1241.59
|
|
||||||
1240.46,1240.46,1240.46
|
|
||||||
1240.4,1240.4,1240.4
|
|
||||||
1233,1233,1233
|
|
||||||
1228.28,1228.28,1228.28
|
|
||||||
1223.75,1223.75,1223.75
|
|
||||||
1223.12,1223.12,1223.12
|
|
||||||
1224.71,1224.71,1224.71
|
|
||||||
1221.53,1221.53,1221.53
|
|
||||||
1206.07,1206.07,1206.07
|
|
||||||
1180.55,1180.55,1180.55
|
|
||||||
1187.76,1187.76,1187.76
|
|
||||||
1189.4,1189.4,1189.4
|
|
||||||
1198.35,1198.35,1198.35
|
|
||||||
1180.73,1180.73,1180.73
|
|
||||||
1197.84,1197.84,1197.84
|
|
||||||
1199.73,1199.73,1199.73
|
|
||||||
1196.69,1196.69,1196.69
|
|
||||||
1178.59,1178.59,1178.59
|
|
||||||
1178.34,1178.34,1178.34
|
|
||||||
1197.75,1197.75,1197.75
|
|
||||||
1199.21,1199.21,1199.21
|
|
||||||
1213.54,1213.54,1213.54
|
|
||||||
1218.71,1218.71,1218.71
|
|
||||||
1213.4,1213.4,1213.4
|
|
||||||
1223.25,1223.25,1223.25
|
|
||||||
1225.85,1225.85,1225.85
|
|
||||||
1221.06,1221.06,1221.06
|
|
||||||
1197.96,1197.96,1197.96
|
|
||||||
1193.57,1193.57,1193.57
|
|
||||||
1184.38,1184.38,1184.38
|
|
||||||
1183.26,1183.26,1183.26
|
|
||||||
1183.78,1183.78,1183.78
|
|
||||||
1182.45,1182.45,1182.45
|
|
||||||
1185.64,1185.64,1185.64
|
|
||||||
1185.62,1185.62,1185.62
|
|
||||||
1183.08,1183.08,1183.08
|
|
||||||
1180.26,1180.26,1180.26
|
|
||||||
1178.17,1178.17,1178.17
|
|
||||||
1165.9,1165.9,1165.9
|
|
||||||
1184.71,1184.71,1184.71
|
|
||||||
1176.19,1176.19,1176.19
|
|
||||||
1173.81,1173.81,1173.81
|
|
||||||
1178.1,1178.1,1178.1
|
|
||||||
1169.77,1169.77,1169.77
|
|
||||||
1165.32,1165.32,1165.32
|
|
||||||
1165.15,1165.15,1165.15
|
|
||||||
1158.06,1158.06,1158.06
|
|
||||||
1159.97,1159.97,1159.97
|
|
||||||
1160.75,1160.75,1160.75
|
|
||||||
1137.03,1137.03,1137.03
|
|
||||||
1146.24,1146.24,1146.24
|
|
||||||
1141.2,1141.2,1141.2
|
|
||||||
1144.73,1144.73,1144.73
|
|
||||||
1147.7,1147.7,1147.7
|
|
||||||
1142.16,1142.16,1142.16
|
|
||||||
1148.67,1148.67,1148.67
|
|
||||||
1124.83,1124.83,1124.83
|
|
||||||
1134.28,1134.28,1134.28
|
|
||||||
1139.78,1139.78,1139.78
|
|
||||||
1142.71,1142.71,1142.71
|
|
||||||
1125.59,1125.59,1125.59
|
|
||||||
1124.66,1124.66,1124.66
|
|
||||||
1125.07,1125.07,1125.07
|
|
||||||
1121.1,1121.1,1121.1
|
|
||||||
1121.9,1121.9,1121.9
|
|
||||||
1109.55,1109.55,1109.55
|
|
||||||
1104.18,1104.18,1104.18
|
|
|
|
@ -1,5 +1,5 @@
|
||||||
"""
|
"""
|
||||||
Utility functions for ????.
|
Signal Filtering/Smoothing and Generation of Synthetic Time-Series.
|
||||||
|
|
||||||
Copyright (c) 2020 Gabriele Gilardi
|
Copyright (c) 2020 Gabriele Gilardi
|
||||||
"""
|
"""
|
||||||
|
@ -7,6 +7,200 @@ Copyright (c) 2020 Gabriele Gilardi
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
|
||||||
|
|
||||||
|
def synthetic_wave(P, A=None, phi=None, num=1000):
|
||||||
|
"""
|
||||||
|
Generates a multi-sine wave given a periods, amplitudes, and phases.
|
||||||
|
|
||||||
|
P (n, ) Periods
|
||||||
|
A (n, ) Amplitudes
|
||||||
|
phi (n, ) Phases (rad)
|
||||||
|
t (num, ) Time
|
||||||
|
f (num, ) Multi-sine wave
|
||||||
|
|
||||||
|
The default value for the amplitudes is 1 and for the phases is zero. The
|
||||||
|
time goes from zero to the largest period.
|
||||||
|
"""
|
||||||
|
n_waves = len(P) # Number of waves
|
||||||
|
P = np.asarray(P)
|
||||||
|
|
||||||
|
# Amplitudes
|
||||||
|
if (A is None):
|
||||||
|
A = np.ones(n_waves) # Default is 1
|
||||||
|
else:
|
||||||
|
A = np.asarray(A)
|
||||||
|
|
||||||
|
# Phases
|
||||||
|
if (phi is None):
|
||||||
|
phi = np.zeros(n_waves) # Default is 0
|
||||||
|
else:
|
||||||
|
phi = np.asarray(phi)
|
||||||
|
|
||||||
|
# Time
|
||||||
|
t = np.linspace(0.0, np.amax(P), num=num)
|
||||||
|
|
||||||
|
# Add up all the sine waves
|
||||||
|
f = np.zeros(len(t))
|
||||||
|
for i in range(n_waves):
|
||||||
|
f = f + A[i] * np.sin(2.0 * np.pi * t / P[i] + phi[i])
|
||||||
|
|
||||||
|
return t, f
|
||||||
|
|
||||||
|
|
||||||
|
def synthetic_FFT(X, n_reps=1):
|
||||||
|
"""
|
||||||
|
Generates surrogates of the time-serie X using the phase-randomized
|
||||||
|
Fourier-transform algorithm.
|
||||||
|
|
||||||
|
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)
|
||||||
|
|
||||||
|
# The number of samples must be odd
|
||||||
|
if ((n % 2) == 0):
|
||||||
|
print("Warning: data reduced by one (even number of samples)")
|
||||||
|
n = n - 1
|
||||||
|
X = X[0:n, :]
|
||||||
|
|
||||||
|
# FFT of the original time-serie
|
||||||
|
X_fft = np.fft.fft(X)
|
||||||
|
|
||||||
|
# Parameters
|
||||||
|
half_len = (n - 1) // 2
|
||||||
|
idx1 = np.arange(1, half_len+1, dtype=int) # 1st half
|
||||||
|
idx2 = np.arange(half_len+1, n, dtype=int) # 2nd half
|
||||||
|
|
||||||
|
# Generate the random phases
|
||||||
|
phases = np.random.rand(n_reps, half_len)
|
||||||
|
phases1 = np.exp(2.0 * np.pi * 1j * phases)
|
||||||
|
phases2 = np.conj(np.flipud(phases1))
|
||||||
|
|
||||||
|
# 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[:, idx1] = X_fft[idx1] * phases1 # 1st half
|
||||||
|
X_synt_fft[:, idx2] = X_fft[idx2] * phases2 # 2nd half
|
||||||
|
|
||||||
|
# Synthetic time-series
|
||||||
|
X_synt = np.real(np.fft.ifft(X_synt_fft, axis=1))
|
||||||
|
|
||||||
|
return X_synt
|
||||||
|
|
||||||
|
|
||||||
|
def synthetic_sampling(X, n_reps=1, replace=True):
|
||||||
|
"""
|
||||||
|
Generates surrogates of the time-serie X using randomized-sampling
|
||||||
|
(bootstrap) with or without replacement.
|
||||||
|
|
||||||
|
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)
|
||||||
|
|
||||||
|
# Sampling with replacement
|
||||||
|
if (replace):
|
||||||
|
idx = np.random.randint(0, n, size=(n_reps, n))
|
||||||
|
|
||||||
|
# Sampling without replacement
|
||||||
|
else:
|
||||||
|
idx = np.argsort(np.random.rand(n_reps, n), axis=1)
|
||||||
|
|
||||||
|
# Synthetic time-series
|
||||||
|
X_synt = X[idx]
|
||||||
|
|
||||||
|
return X_synt
|
||||||
|
|
||||||
|
|
||||||
|
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.
|
||||||
|
|
||||||
|
X (n, ) Original time-series
|
||||||
|
idx (n, ) Original order of X
|
||||||
|
y (n, ) Sorted original time-series
|
||||||
|
z (n+1, ) Intermediate points
|
||||||
|
mt (n, ) Interval means
|
||||||
|
t_w (n_reps, n) Random new points
|
||||||
|
w_int (n_reps, n) Interpolated new points
|
||||||
|
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)
|
||||||
|
|
||||||
|
# Sort the time series keeping track of the original order
|
||||||
|
idx = np.argsort(X)
|
||||||
|
y = X[idx]
|
||||||
|
|
||||||
|
# Trimmed mean
|
||||||
|
g = int(np.floor(n * alpha))
|
||||||
|
r = n * alpha - g
|
||||||
|
m_trm = ((1.0 - r) * (y[g] + y[n-g-1]) + y[g+1:n-g-1].sum()) \
|
||||||
|
/ (n * (1.0 - 2.0 * alpha))
|
||||||
|
|
||||||
|
# Intermediate points
|
||||||
|
z = np.zeros(n+1)
|
||||||
|
z[0] = y[0] - m_trm
|
||||||
|
z[1:-1] = (y[0:-1] + y[1:]) / 2.0
|
||||||
|
z[n] = y[n-1] + m_trm
|
||||||
|
|
||||||
|
# Interval means
|
||||||
|
mt = np.zeros(n)
|
||||||
|
mt[0] = 0.75 * y[0] + 0.25 * y[1]
|
||||||
|
mt[1:n-1] = 0.25 * y[0:n-2] + 0.5 * y[1:n-1] + 0.25 * y[2:n]
|
||||||
|
mt[n-1] = 0.25 * y[n-2] + 0.75 * y[n-1]
|
||||||
|
|
||||||
|
# Generate randomly new points and sort them
|
||||||
|
t_w = np.random.rand(n_reps, n)
|
||||||
|
t_w = np.sort(t_w, axis=1)
|
||||||
|
|
||||||
|
# Interpolate the new points
|
||||||
|
t = np.linspace(0.0, 1.0, num=n+1)
|
||||||
|
w_int = np.interp(t_w, t, z)
|
||||||
|
|
||||||
|
# Correct the new points in the first and last interval to satisfy
|
||||||
|
# the mass constraint
|
||||||
|
idw = (np.floor_divide(t_w, 1.0 / n)).astype(int)
|
||||||
|
corr = np.where(idw == 0, mt[idw] - (z[idw] + z[idw+1]) / 2.0, 0.0)
|
||||||
|
w_corr = w_int + corr
|
||||||
|
if (n > 1):
|
||||||
|
corr = np.where(idw == n-1, mt[idw] - (z[idw] + z[idw+1]) / 2.0, 0.0)
|
||||||
|
w_corr += corr
|
||||||
|
|
||||||
|
# Enforce limits (if requested)
|
||||||
|
if (bounds):
|
||||||
|
w_corr = np.fmin(np.fmax(w_corr, z[0]), z[n])
|
||||||
|
|
||||||
|
# Recovery the time-dependency of the original time-series
|
||||||
|
X_synt = np.zeros((n_reps, n))
|
||||||
|
X_synt[:, idx] = w_corr
|
||||||
|
|
||||||
|
# Scale to force equal variance (if requested)
|
||||||
|
if (scale):
|
||||||
|
var_z = np.diff(z) ** 2.0 / 12.0
|
||||||
|
X_mean = X.mean(axis=0)
|
||||||
|
var_ME = (((mt - X_mean) ** 2).sum() + var_z.sum()) / n
|
||||||
|
std_X = X.std(ddof=1)
|
||||||
|
std_ME = np.sqrt(var_ME)
|
||||||
|
k_scale = std_X / std_ME - 1.0
|
||||||
|
X_synt = X_synt + k_scale * (X_synt - X_mean)
|
||||||
|
|
||||||
|
return X_synt
|
||||||
|
|
||||||
|
|
||||||
def normalize_data(X, param=(), ddof=0):
|
def normalize_data(X, param=(), ddof=0):
|
||||||
"""
|
"""
|
||||||
If mu and sigma are not defined, returns a column-normalized version of
|
If mu and sigma are not defined, returns a column-normalized version of
|
||||||
|
@ -63,235 +257,61 @@ def scale_data(X, param=()):
|
||||||
|
|
||||||
def value2diff(X, percent=True):
|
def value2diff(X, percent=True):
|
||||||
"""
|
"""
|
||||||
from value to difference in abs or %
|
Returns the 1st discrete difference of array X.
|
||||||
dx is reduced by 1 row
|
|
||||||
|
X (n, ) Original 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).
|
||||||
"""
|
"""
|
||||||
# Difference in percent
|
X = X.flatten() # Reshape to (n, )
|
||||||
|
|
||||||
|
# Discrete difference in percent
|
||||||
if (percent):
|
if (percent):
|
||||||
dX = X[1:, :] / X[:-1, :] - 1.0
|
dX = X[1:] / X[:-1] - 1.0
|
||||||
|
|
||||||
# Difference in value
|
# Discrete difference in value
|
||||||
else:
|
else:
|
||||||
dX = X[1:, :] - X[:-1, :]
|
dX = X[1:] - X[:-1]
|
||||||
|
|
||||||
return dX
|
return dX
|
||||||
|
|
||||||
|
|
||||||
def diff2value(dX, percent=True):
|
def diff2value(dX, X0, percent=True):
|
||||||
"""
|
"""
|
||||||
from difference in abs or % to value
|
Returns array X from the 1st discrete difference.
|
||||||
X is increased by one row
|
|
||||||
|
|
||||||
Value from percent: first row set to one. If X0 defines the starting
|
dX (n, ) Discrete differences
|
||||||
values, then X0*X would be the ???
|
X0 scalar Initial value
|
||||||
|
X (n+1, ) Original 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).
|
||||||
|
|
||||||
Value from difference: first row set to zero. If X0 defines the starting
|
If the discrete difference is in percent, the first column of X is set to
|
||||||
values, then X0+X would be the ???
|
one. The original array is X[0] * X
|
||||||
|
|
||||||
|
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).
|
||||||
"""
|
"""
|
||||||
n_rows, n_cols = dX.shape
|
dX = dX.flatten() # Reshape to (n, )
|
||||||
X = np.zeros((n_rows+1, n_cols))
|
n = len(dX)
|
||||||
|
X = np.zeros(n+1)
|
||||||
|
|
||||||
# Value from percent
|
# Discrete difference in percent
|
||||||
# X[0, :] = 1
|
|
||||||
# X[1, :] = X[0, :] * (1 + dX[1, :]) = (1 + dX[1, :])
|
|
||||||
# X[2, :] = X[1, :] * (1 + dX[2, :])
|
|
||||||
# = X[0, :] * (1 + dX[1, :]) * (1 + dX[2, :])
|
|
||||||
# = (1 + dX[1, :]) * (1 + dX[2, :])
|
|
||||||
# ....
|
|
||||||
if (percent):
|
if (percent):
|
||||||
X[0, :] = 1.0
|
X[0] = X0
|
||||||
X[1:, :] = np.cumprod((1.0 + dX), axis=0)
|
X[1:] = np.cumprod(1.0 + dX) * X0 # !!!! check
|
||||||
|
|
||||||
# Value from difference
|
# Discrete difference in value
|
||||||
# X[0, :] = 0
|
|
||||||
# X[1, :] = X[0, :] + dX[1, :] = dX[1, :]
|
|
||||||
# X[2, :] = X[0, :] + dX[1, :] + dX[2, :] = dX[1, :] + dX[2, :]
|
|
||||||
# ....
|
|
||||||
else:
|
else:
|
||||||
# First row already set to zero
|
X[0] = X0
|
||||||
X[1:, :] = np.cumsum(dX, axis=0)
|
X[1:] = np.cumsum(dX) + X0 # !!!!!c check
|
||||||
|
|
||||||
return X
|
return X
|
||||||
|
|
||||||
|
|
||||||
def synthetic_wave(P, A=None, phi=None, num=1000):
|
|
||||||
"""
|
|
||||||
Generates a multi-sinewave.
|
|
||||||
P = [P_1, P_2, ... P_n] Periods
|
|
||||||
A = [A_1, A_2, ... A_n] Amplitudes
|
|
||||||
phi = [phi_1, phi_2, ... phi_n] Phases (rad)
|
|
||||||
|
|
||||||
Default amplitudes are ones
|
|
||||||
Default phases are zeros
|
|
||||||
Time is from 0 to largest period (default 1000 steps)
|
|
||||||
"""
|
|
||||||
n_waves = len(P)
|
|
||||||
P = np.asarray(P)
|
|
||||||
|
|
||||||
# Define amplitudes
|
|
||||||
if (A is None):
|
|
||||||
A = np.ones(n_waves)
|
|
||||||
else:
|
|
||||||
A = np.asarray(A)
|
|
||||||
|
|
||||||
# Define phases
|
|
||||||
if (phi is None):
|
|
||||||
phi = np.zeros(n_waves)
|
|
||||||
else:
|
|
||||||
phi = np.asarray(phi)
|
|
||||||
|
|
||||||
# Add all the waves
|
|
||||||
t = np.linspace(0.0, np.amax(P), num=num)
|
|
||||||
f = np.zeros(len(t))
|
|
||||||
for i in range(n_waves):
|
|
||||||
f = f + A[i] * np.sin(2.0 * np.pi * t / P[i] + phi[i])
|
|
||||||
|
|
||||||
return t, f
|
|
||||||
|
|
||||||
|
|
||||||
def synthetic_FFT(X, multiv=False):
|
|
||||||
"""
|
|
||||||
- univariate and single time-series
|
|
||||||
- univariate and multi-time series (can be used to generate multi from same)
|
|
||||||
- multi-variate multi-time series
|
|
||||||
|
|
||||||
generate more than n_samples, remome multi-time series for uni-variate
|
|
||||||
use 3d matrix looping on the time-series for multi-variate
|
|
||||||
|
|
||||||
"""
|
|
||||||
n_samples, n_series = X.shape
|
|
||||||
|
|
||||||
# The number of samples must be odd (if the number is even drop the last value)
|
|
||||||
if ((n_samples % 2) == 0):
|
|
||||||
print("Warning: data reduced by one (even number of samples)")
|
|
||||||
n_samples = n_samples - 1
|
|
||||||
X = X[0:n_samples, :]
|
|
||||||
|
|
||||||
# FFT of the original data
|
|
||||||
X_fft = np.fft.fft(X, axis=0)
|
|
||||||
|
|
||||||
# Parameters
|
|
||||||
half_len = (n_samples - 1) // 2
|
|
||||||
idx1 = np.arange(1, half_len+1, dtype=int)
|
|
||||||
idx2 = np.arange(half_len+1, n_samples, dtype=int)
|
|
||||||
|
|
||||||
# If multivariate the random phases is the same
|
|
||||||
if (multiv):
|
|
||||||
phases = np.random.rand(half_len, 1)
|
|
||||||
phases1 = np.tile(np.exp(2.0 * np.pi * 1j * phases), (1, n_series))
|
|
||||||
phases2 = np.conj(np.flipud(phases1))
|
|
||||||
|
|
||||||
# If univariate the random phases is different
|
|
||||||
else:
|
|
||||||
phases = np.random.rand(half_len, n_series)
|
|
||||||
phases1 = np.exp(2.0 * np.pi * 1j * phases)
|
|
||||||
phases2 = np.conj(np.flipud(phases1))
|
|
||||||
|
|
||||||
# FFT of the synthetic data
|
|
||||||
synt_fft = X_fft.copy()
|
|
||||||
synt_fft[idx1, :] = X_fft[idx1, :] * phases1
|
|
||||||
synt_fft[idx2, :] = X_fft[idx2, :] * phases2
|
|
||||||
|
|
||||||
# Inverse FFT of the synthetic data
|
|
||||||
X_synt = np.real(np.fft.ifft(synt_fft, axis=0))
|
|
||||||
|
|
||||||
return X_synt
|
|
||||||
|
|
||||||
|
|
||||||
def synthetic_sampling(X, n_reps=1, replace=True):
|
|
||||||
"""
|
|
||||||
X must be (n_samples, 1)
|
|
||||||
"""
|
|
||||||
n_samples = X.shape[0]
|
|
||||||
|
|
||||||
# Sampling with replacement
|
|
||||||
if (replace):
|
|
||||||
idx = np.random.randint(0, n_samples, size=(n_samples, n_reps))
|
|
||||||
i = np.arange(n_reps)
|
|
||||||
Xe = np.tile(X,(1, n_reps))
|
|
||||||
X_synt = Xe[idx, i]
|
|
||||||
|
|
||||||
# Sampling without replacement
|
|
||||||
else:
|
|
||||||
idx = np.zeros(n_samples, n_reps)
|
|
||||||
for j in range(n_reps):
|
|
||||||
idx[:, j] = np.random.permutation(n_samples)
|
|
||||||
i = np.arange(n_reps)
|
|
||||||
Xe = np.tile(X,(1, n_reps))
|
|
||||||
X_synt = Xe[idx, i]
|
|
||||||
|
|
||||||
return X_synt
|
|
||||||
|
|
||||||
|
|
||||||
def synthetic_MEboot(X, alpha=0.1, bounds=True, scale=False):
|
|
||||||
"""
|
|
||||||
"""
|
|
||||||
n_samples, n_series = X.shape
|
|
||||||
X_synt = np.zeros_like(X)
|
|
||||||
|
|
||||||
# Loop over time-series iterate over number of desired series or set parallel
|
|
||||||
# if possile
|
|
||||||
n = n_samples
|
|
||||||
for ts in range(n_series):
|
|
||||||
|
|
||||||
# Sort the time series keeping track of the original position
|
|
||||||
idx = np.argsort(X[:, ts])
|
|
||||||
Y = X[idx, ts]
|
|
||||||
|
|
||||||
# Compute the trimmed mean
|
|
||||||
# trimmed_mean = np.abs(np.diff(x)).mean()
|
|
||||||
g = int(np.floor(n * alpha))
|
|
||||||
r = n * alpha - g
|
|
||||||
m_trm = ((1.0 - r) * (Y[g] + Y[n-g-1]) + Y[g+1:n-g-1].sum()) \
|
|
||||||
/ (n * (1.0 - 2.0 * alpha))
|
|
||||||
|
|
||||||
# Compute the intermediate points
|
|
||||||
Z = np.zeros(n+1)
|
|
||||||
Z[0] = Y[0] - m_trm
|
|
||||||
Z[1:-1] = (Y[0:-1] + Y[1:]) / 2.0
|
|
||||||
Z[n] = Y[n-1] + m_trm
|
|
||||||
|
|
||||||
# Compute the interval means
|
|
||||||
mt = np.zeros(n)
|
|
||||||
mt[0] = 0.75 * Y[0] + 0.25 * Y[1]
|
|
||||||
mt[1:n-1] = 0.25 * Y[0:n-2] + 0.5 * Y[1:n-1] + 0.25 * Y[2:n]
|
|
||||||
mt[n-1] = 0.25 * Y[n-2] + 0.75 * Y[n-1]
|
|
||||||
|
|
||||||
# Randomly generate new points
|
|
||||||
# t_w = np.random.rand(n)
|
|
||||||
t_w = np.array([0.12, 0.83, 0.53, 0.59, 0.11])
|
|
||||||
t_w = np.sort(t_w)
|
|
||||||
|
|
||||||
# Interpolate new points
|
|
||||||
t = np.linspace(0.0, 1.0, num=n+1)
|
|
||||||
w_int = np.interp(t_w, t, Z)
|
|
||||||
|
|
||||||
# Correct the new points (first and last interval) to satisfy mass constraint
|
|
||||||
idw = (np.floor_divide(t_w, 1.0 / n)).astype(int)
|
|
||||||
corr = np.where(idw == 0, mt[idw] - (Z[idw] + Z[idw+1]) / 2.0, 0.0)
|
|
||||||
w_corr = w_int + corr
|
|
||||||
if (n > 1):
|
|
||||||
corr = np.where(idw == n-1, mt[idw] - (Z[idw] + Z[idw+1]) / 2.0, 0.0)
|
|
||||||
w_corr += corr
|
|
||||||
|
|
||||||
# Enforce limits
|
|
||||||
if (bounds):
|
|
||||||
w_corr = np.fmin(np.fmax(w_corr, Z[0]), Z[n])
|
|
||||||
|
|
||||||
# Recovery the time-dependencies (done all togheter?)
|
|
||||||
X_synt = np.zeros(n)
|
|
||||||
X_synt[idx] = w_corr
|
|
||||||
|
|
||||||
# Scale (done all together?)
|
|
||||||
if (scale):
|
|
||||||
var_Z = np.diff(Z) ** 2.0 / 12.0
|
|
||||||
X_mean = X.mean(axis=0)
|
|
||||||
var_ME = (((mt - X_mean) ** 2).sum() + var_Z.sum()) / n
|
|
||||||
std_X= X.std(axis=0, ddof=1)
|
|
||||||
std_ME = np.sqrt(var_ME)
|
|
||||||
k_scale = std_X / std_ME - 1.0
|
|
||||||
X_synt = X_synt + k_scale * (X_synt - X_mean)
|
|
||||||
|
|
||||||
return X_synt
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -1,20 +1,14 @@
|
||||||
"""
|
"""
|
||||||
Filters for time series.
|
Signal Filtering/Smoothing and Generation of Synthetic Time-Series.
|
||||||
|
|
||||||
Copyright (c) 2020 Gabriele Gilardi
|
Copyright (c) 2020 Gabriele Gilardi
|
||||||
|
|
||||||
|
|
||||||
ToDo:
|
ToDo:
|
||||||
- add comments to the code
|
- add comments to the code
|
||||||
- in comments write what filters do
|
- in comments write what filters do
|
||||||
- is necessary to copy X for Y untouched?
|
- is necessary to copy X for Y untouched?
|
||||||
- decide default values in functions
|
- decide default values in functions
|
||||||
- why lag plot gives errors
|
|
||||||
- fix plotting function
|
|
||||||
- example for alpha-beta-gamma using variable sigma as in financial time series
|
|
||||||
(see Ehler)
|
|
||||||
- example using noisy multi-sine-waves
|
|
||||||
- vectors must be ( .., 1)
|
|
||||||
- synt remove multi-series unless multi-variate and add number of istances (fft 3 D)
|
|
||||||
"""
|
"""
|
||||||
|
|
||||||
import sys
|
import sys
|
||||||
|
@ -24,26 +18,32 @@ import matplotlib.pyplot as plt
|
||||||
import filters as flt
|
import filters as flt
|
||||||
import synthetic as syn
|
import synthetic as syn
|
||||||
|
|
||||||
|
# Added to avoid message warning "The group delay is singular at frequencies
|
||||||
|
# [....], setting to 0" when plotting the lag for SMA filters.
|
||||||
|
np.seterr(all='ignore')
|
||||||
|
np.random.seed(1294404794)
|
||||||
|
|
||||||
# Read data to filter
|
# Read data to filter
|
||||||
if len(sys.argv) != 2:
|
if len(sys.argv) != 2:
|
||||||
print("Usage: python test.py <data_file>")
|
print("Usage: python test.py <data_file>")
|
||||||
sys.exit(1)
|
sys.exit(1)
|
||||||
data_file = sys.argv[1] + '.csv'
|
data_file = sys.argv[1] + '.csv'
|
||||||
|
|
||||||
# Read data from a csv file
|
# Read data from a csv file (one time-series each column)
|
||||||
data = np.loadtxt(data_file, delimiter=',')
|
data = np.loadtxt(data_file, delimiter=',')
|
||||||
n_samples = data.shape[0]
|
|
||||||
data = data.reshape(n_samples, -1)
|
|
||||||
|
|
||||||
np.random.seed(1294404794)
|
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)
|
||||||
|
|
||||||
|
# spx.plot_frequency()
|
||||||
|
# spx.plot_lag()
|
||||||
|
|
||||||
# res, bb, aa = spx.SincFunction(2, 50)
|
|
||||||
# print(bb)
|
|
||||||
# print(aa)
|
|
||||||
# utl.plot_frequency_response(bb, aa)
|
|
||||||
# utl.plot_lag_response(bb, aa)
|
|
||||||
# sigma_x = 0.1
|
# sigma_x = 0.1
|
||||||
# sigma_v = 0.1 * np.ones(n_samples)
|
# sigma_v = 0.1 * np.ones(n_samples)
|
||||||
# res = spx.Kalman(sigma_x=sigma_x, sigma_v=sigma_v, dt=1.0, abg_type="abg")
|
# res = spx.Kalman(sigma_x=sigma_x, sigma_v=sigma_v, dt=1.0, abg_type="abg")
|
||||||
|
@ -71,8 +71,9 @@ np.random.seed(1294404794)
|
||||||
# [-0.8236, -0.8314],
|
# [-0.8236, -0.8314],
|
||||||
# [-1.5771, -0.9792],
|
# [-1.5771, -0.9792],
|
||||||
# [ 0.5080, -1.1564]])
|
# [ 0.5080, -1.1564]])
|
||||||
# synt_data1 = syn.synthetic_FFT(data, False)
|
# synt_data = syn.synthetic_FFT(aa[0:5,0], n_reps=1)
|
||||||
# synt_data2 = syn.synthetic_FFT(data, False)
|
# print(synt_data)
|
||||||
|
|
||||||
# plt.plot(synt_data1)
|
# plt.plot(synt_data1)
|
||||||
# plt.plot(synt_data2)
|
# plt.plot(synt_data2)
|
||||||
# plt.plot(data)
|
# plt.plot(data)
|
||||||
|
@ -85,30 +86,14 @@ np.random.seed(1294404794)
|
||||||
# print(bb[0:10, :])
|
# print(bb[0:10, :])
|
||||||
# cc = syn.diff2value(bb, percent)
|
# cc = syn.diff2value(bb, percent)
|
||||||
# print(cc[0:10, :]+1399.48)
|
# print(cc[0:10, :]+1399.48)
|
||||||
# aa = np.arange(10,28).reshape(6,3)
|
|
||||||
# print(aa)
|
|
||||||
# idx = np.zeros_like(aa)
|
|
||||||
# bb = np.zeros_like(aa)
|
|
||||||
# for i in range(aa.shape[1]):
|
|
||||||
# idx[:, i] = np.random.permutation(aa.shape[0])
|
|
||||||
# print(idx)
|
|
||||||
# i = np.arange(aa.shape[1])
|
# i = np.arange(aa.shape[1])
|
||||||
# bb[:, i] = aa[idx[:, i], i]
|
# bb[:, i] = aa[idx[:, i], i]
|
||||||
# bb = syn.synthetic_boot(aa, replace=False)
|
# aa = np.array([4, 12, 36, 20, 8])
|
||||||
|
# bb = syn.synthetic_sampling(aa, n_reps=2, replace=True)
|
||||||
# print(bb)
|
# print(bb)
|
||||||
# aa = np.array([4, 12, 36, 20, 8]).reshape(5, 1)
|
# aa = np.array([4, 12, 36, 20, 8])
|
||||||
# W = syn.synthetic_MEboot(aa, alpha=0.1, bounds=False, scale=True)
|
# W = syn.synthetic_MEboot(aa, n_reps=1, alpha=0.1, bounds=False, scale=False)
|
||||||
# print(bb.sum())
|
# print('W=')
|
||||||
# print('W=', W)
|
# print(W)
|
||||||
|
|
||||||
n = 8
|
|
||||||
aa = np.arange(n).reshape(n,1) * 1.1
|
|
||||||
print(aa)
|
|
||||||
idx = np.random.randint(0, n, size=(n, 3))
|
|
||||||
print(idx)
|
|
||||||
i = np.arange(3)
|
|
||||||
# print(i)
|
|
||||||
bb = np.tile(aa,(1, 3))
|
|
||||||
# print(bb)
|
|
||||||
cc = bb[idx, i]
|
|
||||||
print(cc)
|
|
Ładowanie…
Reference in New Issue