{{ message }}

Instantly share code, notes, and snippets.

# tartakynov/fourex.py

Last active Dec 26, 2020
Fourier Extrapolation in Python
 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 * 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 * 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()

### Vayel commented Jun 7, 2015

 Line 15: `AttributeError: 'range' object has no attribute 'sort'`. Use `indexes = list(range(n))` instead.

### treviesweets commented Jun 11, 2015

 Thanks!

### etamponi commented Jul 6, 2015

 Why do you apply detrending on the original signal?

### tartakynov commented Jul 22, 2015

 @etamponi if you don't do detrending and your signal goes gradually up or down (e.g. has trend) then you will get zigzags on extrapolation, just try it yourself :)

### LetterRip commented Jun 1, 2016 • edited

 Hi, found a bug, but not clear on the solution. If I generate this synthetic series and use it with your code above, the prediction can be excellent or awful depending on when I extrapolate from. `````` X = [i for i in range(365*5)] slope = 5./365 b = 0.0 trend_Y = [i * slope + b for i in X] seasonal_Y = [10*cos(2*pi * i /365. ) for i in X] monthly_Y = [5*cos(2*pi * i /30. ) for i in X] combined_Y = [a+b+c for (a,b,c) in zip(trend_Y, seasonal_Y, monthly_Y)] x = np.array(combined_Y) n_predict = 180 #355 almost perfect, 180 awful extrapolation = fourierExtrapolation(x[:-n_predict], n_predict) pl.plot(np.arange(0, x.size), x, 'b', label = 'x', linewidth = 3) pl.plot(np.arange(0, extrapolation.size), extrapolation, 'r', label = 'extrapolation') `````` 180 day extrapolation http://postimg.org/image/5cuepv0h7/ 355 day extrapolation http://postimg.org/image/krodvv7l7/

### gbrand-salesforce commented Nov 15, 2016

 Why don't you zero out the irrelevant component in the frequency domain and inverse fft instead of iterating over the relevant frequencies and constructing the signal by adding cosines? Your code runs at O(kn), where k is the number of relevant frequencies while ifft is O(nlog(n)), so if the number of frequencies your interested in is higher the log(n), it'd be better to do an ifft instead.

### bqzz commented Jan 5, 2018

 @LetterRip, it may relate to spectrum leakage because you didn't include integer multiples samples of your signal's periods in the FFT. I guess the jump effect you watched should be more obvious when there're only a few frequencies in your signal.

### grx7 commented May 2, 2018

 @gbrand-salesforce what would be the better version?

### keithdlandry commented Nov 15, 2018

 Nice work. Why don't you sort by the amplitudes instead of the frequencies? It would be better to take the `n_harm` harmonics with the largest amplitudes instead of the `n_harm` harmonics with the lowest frequencies. Basically just change line fifteen to the following: ``````indexes.sort(key=lambda i: np.absolute(x_freqdom[i])) indexes.reverse() `````` You can see the difference in the following examples. 40 harmonics It's even more clear when you try to capture trends with a minimum number of harmonics. With amplitude sorting just 4 harmonics can fit the data nicely. The other sorting methods don't capture anything at all. 4 harmonics

### liuguiyangnwpu commented Dec 11, 2018

 @etamponi DC component needs to be removed during Fourier transform

### PushkhallaChandramoulli commented Jan 18, 2019

 Hi, Can you please enlighten me on how did you assign no. of harmonics (n_harm ) as 10 ?

### samy-n commented Jan 11, 2020

 On the x axis, how can I have the timestamps column values instead of numbers?

### deathlyhallows010 commented Feb 10, 2020

 Can we extend the band-limited amplitude spectrum of a zero-phase wavelet reliably using this method

### PRANEETH-ALURU commented Feb 12, 2020

 indexes = np.array(range(0,n)) sorted(indexes,key = lambda i: np.absolute(f[i]))

### Okeus commented May 1, 2020 • edited

 #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) # 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 * 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 * t def main(): x = np.loadtxt(sys.argv,delimiter=',') x=x.ravel() n_predict = int(sys.argv) 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() ``````

### ZainabL commented May 4, 2020

 Nice coding !

### FaridAbrari commented Aug 17, 2020

 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.

### alik604 commented Aug 24, 2020

 Line #13 should `indexes = list(range(n))`