From 753245f28579dcf838e0fa8fde12e7490744dbe3 Mon Sep 17 00:00:00 2001 From: Gabriele Gilardi Date: Thu, 25 Jun 2020 21:44:24 +0900 Subject: [PATCH] add comments --- Code_Python/filters.py | 47 ++-- Code_Python/spx3.csv | 500 --------------------------------------- Code_Python/synthetic.py | 458 ++++++++++++++++++----------------- Code_Python/test.py | 69 +++--- 4 files changed, 290 insertions(+), 784 deletions(-) delete mode 100644 Code_Python/spx3.csv diff --git a/Code_Python/filters.py b/Code_Python/filters.py index de79e20..bff0ef0 100644 --- a/Code_Python/filters.py +++ b/Code_Python/filters.py @@ -1,14 +1,9 @@ """ -Class for filter/smooth data. +Signal Filtering/Smoothing and Generation of Synthetic Time-Series. 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 b (n_b, ) Numerator coefficients a (n_a, ) Denominator coefficients @@ -75,23 +70,26 @@ Y (n_obs, 1) Observation residual (innovation) import numpy as np 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): - idx_end = len(signals[0]) - t = np.arange(idx_start, idx_end) - names = [] + legend = [] count = 0 for signal in signals: - plt.plot(t, signal[idx_start:idx_end]) - names.append(str(count)) + signal = signal.flatten() + end = len(signal) + t = np.arange(start, end) + plt.plot(t, signal[start:end]) + legend.append('Signal [' + str(count) + ']') count += 1 + plt.xlabel('Index') + plt.ylabel('Value') plt.grid(b=True) - plt.legend(names) + plt.legend(legend) plt.show() @@ -99,22 +97,23 @@ def filter_data(data, b, a): """ Applies a filter with transfer response coefficients and . """ - n_samples, n_series = data.shape + n_samples = len(data) nb = len(b) na = len(a) idx = np.amax([0, nb-1, na-1]) Y = data.copy() for i in range(idx, n_samples): - tmp = np.zeros(n_series) + + tmp = 0 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): - 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 @@ -125,7 +124,7 @@ class Filter: """ """ 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): """ @@ -577,8 +576,10 @@ class Filter: plt.axhline(-3.0, lw=1.5, ls='--', C='r') plt.grid(b=True) plt.xlim(np.amin(wf), np.amax(wf)) - plt.xlabel('$omega$ [rad/sample]') + plt.xlabel(r'$\omega$ [rad/sample]') plt.ylabel('$h$ [db]') + legend = ['Filter', '-3dB'] + plt.legend(legend) plt.show() def plot_lag(self): @@ -589,6 +590,6 @@ class Filter: plt.plot(wf, gd) plt.grid(b=True) plt.xlim(np.amin(wf), np.amax(wf)) - plt.xlabel('$omega$ [rad/sample]') + plt.xlabel(r'$\omega$ [rad/sample]') plt.ylabel('$gd$ [samples]') plt.show() diff --git a/Code_Python/spx3.csv b/Code_Python/spx3.csv deleted file mode 100644 index ed0bd7b..0000000 --- a/Code_Python/spx3.csv +++ /dev/null @@ -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 diff --git a/Code_Python/synthetic.py b/Code_Python/synthetic.py index 55c7b2c..2401e0f 100644 --- a/Code_Python/synthetic.py +++ b/Code_Python/synthetic.py @@ -1,5 +1,5 @@ """ -Utility functions for ????. +Signal Filtering/Smoothing and Generation of Synthetic Time-Series. Copyright (c) 2020 Gabriele Gilardi """ @@ -7,6 +7,200 @@ Copyright (c) 2020 Gabriele Gilardi 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): """ 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): """ - from value to difference in abs or % - dx is reduced by 1 row + Returns the 1st discrete difference of array X. + + 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): - dX = X[1:, :] / X[:-1, :] - 1.0 - - # Difference in value + dX = X[1:] / X[:-1] - 1.0 + + # Discrete difference in value else: - dX = X[1:, :] - X[:-1, :] - + dX = X[1:] - X[:-1] + return dX -def diff2value(dX, percent=True): +def diff2value(dX, X0, percent=True): """ - from difference in abs or % to value - X is increased by one row + Returns array X from the 1st discrete difference. - Value from percent: first row set to one. If X0 defines the starting - values, then X0*X would be the ??? + dX (n, ) Discrete differences + 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 - values, then X0+X would be the ??? + 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 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 - X = np.zeros((n_rows+1, n_cols)) + dX = dX.flatten() # Reshape to (n, ) + n = len(dX) + X = np.zeros(n+1) - # Value from 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, :]) - # .... + # Discrete difference in percent if (percent): - X[0, :] = 1.0 - X[1:, :] = np.cumprod((1.0 + dX), axis=0) - - # Value from difference - # X[0, :] = 0 - # X[1, :] = X[0, :] + dX[1, :] = dX[1, :] - # X[2, :] = X[0, :] + dX[1, :] + dX[2, :] = dX[1, :] + dX[2, :] - # .... + X[0] = X0 + X[1:] = np.cumprod(1.0 + dX) * X0 # !!!! check + + # Discrete difference in value else: - # First row already set to zero - X[1:, :] = np.cumsum(dX, axis=0) - + X[0] = X0 + X[1:] = np.cumsum(dX) + X0 # !!!!!c check + 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 - - diff --git a/Code_Python/test.py b/Code_Python/test.py index cf800d4..ed32bc1 100644 --- a/Code_Python/test.py +++ b/Code_Python/test.py @@ -1,20 +1,14 @@ """ -Filters for time series. +Signal Filtering/Smoothing and Generation of Synthetic Time-Series. 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 -- 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 @@ -24,26 +18,32 @@ import matplotlib.pyplot as plt import filters as flt 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 if len(sys.argv) != 2: print("Usage: python test.py ") sys.exit(1) 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=',') -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) +# 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_v = 0.1 * np.ones(n_samples) # 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], # [-1.5771, -0.9792], # [ 0.5080, -1.1564]]) -# synt_data1 = syn.synthetic_FFT(data, False) -# synt_data2 = syn.synthetic_FFT(data, False) +# synt_data = syn.synthetic_FFT(aa[0:5,0], n_reps=1) +# print(synt_data) + # plt.plot(synt_data1) # plt.plot(synt_data2) # plt.plot(data) @@ -85,30 +86,14 @@ np.random.seed(1294404794) # print(bb[0:10, :]) # cc = syn.diff2value(bb, percent) # 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]) # 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) -# aa = np.array([4, 12, 36, 20, 8]).reshape(5, 1) -# W = syn.synthetic_MEboot(aa, alpha=0.1, bounds=False, scale=True) -# print(bb.sum()) -# print('W=', W) +# aa = np.array([4, 12, 36, 20, 8]) +# W = syn.synthetic_MEboot(aa, n_reps=1, alpha=0.1, bounds=False, scale=False) +# print('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) \ No newline at end of file