-
Star
(116)
You must be signed in to star a gist -
Fork
(33)
You must be signed in to fork a gist
-
-
Save tartakynov/83f3cd8f44208a1856ce to your computer and use it in GitHub Desktop.
| import numpy as np | |
| import pylab as pl | |
| from numpy import fft | |
| def fourierExtrapolation(x, n_predict): | |
| n = x.size | |
| n_harm = 10 # number of harmonics in model | |
| t = np.arange(0, n) | |
| p = np.polyfit(t, x, 1) # find linear trend in x | |
| x_notrend = x - p[0] * t # detrended x | |
| x_freqdom = fft.fft(x_notrend) # detrended x in frequency domain | |
| f = fft.fftfreq(n) # frequencies | |
| indexes = range(n) | |
| # sort indexes by frequency, lower -> higher | |
| indexes.sort(key = lambda i: np.absolute(f[i])) | |
| t = np.arange(0, n + n_predict) | |
| restored_sig = np.zeros(t.size) | |
| for i in indexes[:1 + n_harm * 2]: | |
| ampli = np.absolute(x_freqdom[i]) / n # amplitude | |
| phase = np.angle(x_freqdom[i]) # phase | |
| restored_sig += ampli * np.cos(2 * np.pi * f[i] * t + phase) | |
| return restored_sig + p[0] * t | |
| def main(): | |
| x = np.array([669, 592, 664, 1005, 699, 401, 646, 472, 598, 681, 1126, 1260, 562, 491, 714, 530, 521, 687, 776, 802, 499, 536, 871, 801, 965, 768, 381, 497, 458, 699, 549, 427, 358, 219, 635, 756, 775, 969, 598, 630, 649, 722, 835, 812, 724, 966, 778, 584, 697, 737, 777, 1059, 1218, 848, 713, 884, 879, 1056, 1273, 1848, 780, 1206, 1404, 1444, 1412, 1493, 1576, 1178, 836, 1087, 1101, 1082, 775, 698, 620, 651, 731, 906, 958, 1039, 1105, 620, 576, 707, 888, 1052, 1072, 1357, 768, 986, 816, 889, 973, 983, 1351, 1266, 1053, 1879, 2085, 2419, 1880, 2045, 2212, 1491, 1378, 1524, 1231, 1577, 2459, 1848, 1506, 1589, 1386, 1111, 1180, 1075, 1595, 1309, 2092, 1846, 2321, 2036, 3587, 1637, 1416, 1432, 1110, 1135, 1233, 1439, 894, 628, 967, 1176, 1069, 1193, 1771, 1199, 888, 1155, 1254, 1403, 1502, 1692, 1187, 1110, 1382, 1808, 2039, 1810, 1819, 1408, 803, 1568, 1227, 1270, 1268, 1535, 873, 1006, 1328, 1733, 1352, 1906, 2029, 1734, 1314, 1810, 1540, 1958, 1420, 1530, 1126, 721, 771, 874, 997, 1186, 1415, 973, 1146, 1147, 1079, 3854, 3407, 2257, 1200, 734, 1051, 1030, 1370, 2422, 1531, 1062, 530, 1030, 1061, 1249, 2080, 2251, 1190, 756, 1161, 1053, 1063, 932, 1604, 1130, 744, 930, 948, 1107, 1161, 1194, 1366, 1155, 785, 602, 903, 1142, 1410, 1256, 742, 985, 1037, 1067, 1196, 1412, 1127, 779, 911, 989, 946, 888, 1349, 1124, 761, 994, 1068, 971, 1157, 1558, 1223, 782, 2790, 1835, 1444, 1098, 1399, 1255, 950, 1110, 1345, 1224, 1092, 1446, 1210, 1122, 1259, 1181, 1035, 1325, 1481, 1278, 769, 911, 876, 877, 950, 1383, 980, 705, 888, 877, 638, 1065, 1142, 1090, 1316, 1270, 1048, 1256, 1009, 1175, 1176, 870, 856, 860]) | |
| n_predict = 100 | |
| extrapolation = fourierExtrapolation(x, n_predict) | |
| pl.plot(np.arange(0, extrapolation.size), extrapolation, 'r', label = 'extrapolation') | |
| pl.plot(np.arange(0, x.size), x, 'b', label = 'x', linewidth = 3) | |
| pl.legend() | |
| pl.show() | |
| if __name__ == "__main__": | |
| main() |
Can we extend the band-limited amplitude spectrum of a zero-phase wavelet reliably using this method
indexes = np.array(range(0,n))
sorted(indexes,key = lambda i: np.absolute(f[i]))
#Great script. Below is code to load the input array from a csv file and adjust the number of predictions and number of harmonics from the command line.
import numpy as np
import pylab as pl
from numpy import fft
import sys
#Example Usage: python fourex.py inputFile.csv numberOfPredictions numberOfHarmonics
#Example Usage: python fourex.py inputFile.csv 100 10
def fourierExtrapolation(x, n_predict):
n = x.size
n_harm = int(sys.argv[3]) # number of harmonics in model
t = np.arange(0, n)
p = np.polyfit(t, x, 1) # find linear trend in x
x_notrend = x - p[0] * t # detrended x
x_freqdom = fft.fft(x_notrend) # detrended x in frequency domain
f = fft.fftfreq(n) # frequencies
indexes = range(n)
# sort indexes by frequency, lower -> higher
indexes.sort(key = lambda i: np.absolute(f[i]))
t = np.arange(0, n + n_predict)
restored_sig = np.zeros(t.size)
for i in indexes[:1 + n_harm * 2]:
ampli = np.absolute(x_freqdom[i]) / n # amplitude
phase = np.angle(x_freqdom[i]) # phase
restored_sig += ampli * np.cos(2 * np.pi * f[i] * t + phase)
return restored_sig + p[0] * t
def main():
x = np.loadtxt(sys.argv[1],delimiter=',')
x=x.ravel()
n_predict = int(sys.argv[2])
extrapolation = fourierExtrapolation(x, n_predict)
pl.plot(np.arange(0, x.size), x, 'b', label = 'inputFile', linewidth = 3)
pl.plot(np.arange(0, extrapolation.size), extrapolation, 'r', label = 'extrapolation')
pl.legend()
pl.show()
if __name__ == "__main__":
main()
Nice coding !
Hello - nice code indeed. But when I try the code in Python v3.8.5 I get a value error like this:
“ValueError: Operands could not be broadcast together with shapes (256,) (251,)”
The operands referring to in the code are: f and t. The size of f is n, whereas the size of t is n+n_predict.
The code runs with n_predict=0, but that defeats the purpose of the code to extrapolate the time series.
Has anyone come across with this error in past? Any suggestion to resolve the error?
Thank you for your help.
Line #13 should indexes = list(range(n))
Hello - nice code indeed. But when I try the code in Python v3.8.5 I get a value error like this:
“ValueError: Operands could not be broadcast together with shapes (256,) (251,)”
The operands referring to in the code are: f and t. The size of f is n, whereas the size of t is n+n_predict.
The code runs with n_predict=0, but that defeats the purpose of the code to extrapolate the time series.
Has anyone come across with this error in past? Any suggestion to resolve the error?
Thank you for your help.
Facing the same issue..
@bqzz, @LetterRip - wouldn't that be Gibb's phenomenon ?
@grx7, @gbrand-salesforce - somethink like this I imagine:
from math import cos,pi
import numpy as np
import matplotlib.pyplot as plt
# Generate Seasonal Data
X = [i for i in range(360 * 3)]
slope = 1 / 365
t_Y = [i * slope for i in X]
s_Y = [60 * cos(2 * pi * i /365) for i in X]
m_Y = [30 * cos(2 * pi * i /30) for i in X]
c_Y = [a + b + c for (a,b,c) in zip(t_Y, s_Y, m_Y)]
x = np.array(c_Y)
# Fast Fourier Transform
fft = np.fft.fft(x)
plt.figure(figsize=(14, 7), dpi=100)
for num_ in [6]:
fft_list = np.copy(fft)
fft_list[num_:-num_] = 0
# Inverse Fast Fourier transform
t = np.fft.ifft(fft_list)
# The trend is your friend
plt.plot(np.concatenate([t,t]), color = 'red')
plt.plot(np.arange(0, x.size), x, 'b', label = 'x', linewidth = 1)
plt.show()
@keithdlandry
Can you explain Why sorting by highest amplitude is better than frequency?
@DanielPark827
Sorting the largest amplitude will give you a better frequency selection because the amplitude in the frequency domain indicates how much each frequency contributes to the value of the original data. So selecting larger amplitudes will give you frequencies that predict more of the pattern in your data.
On the x axis, how can I have the timestamps column values instead of numbers?