Skip to content

Instantly share code, notes, and snippets.

@tartakynov
Last active July 15, 2024 16:24
Show Gist options
  • Save tartakynov/83f3cd8f44208a1856ce to your computer and use it in GitHub Desktop.
Save tartakynov/83f3cd8f44208a1856ce to your computer and use it in GitHub Desktop.
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[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()
@Vayel
Copy link

Vayel commented Jun 7, 2015

Line 15: AttributeError: 'range' object has no attribute 'sort'.

Use indexes = list(range(n)) instead.

@treviesweets
Copy link

Thanks!

@etamponi
Copy link

etamponi commented Jul 6, 2015

Why do you apply detrending on the original signal?

@tartakynov
Copy link
Author

@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 :)

Copy link

ghost commented Jun 1, 2016

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
Copy link

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
Copy link

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
Copy link

grx7 commented May 2, 2018

@gbrand-salesforce what would be the better version?

@keithdlandry
Copy link

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

@guiyang882
Copy link

@etamponi
DC component needs to be removed during Fourier transform

@guiyang882
Copy link

not Well for this signal
image

@PushkhallaChandramoulli

Hi,

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

@samy-n
Copy link

samy-n commented Jan 11, 2020

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

@deathlyhallows010
Copy link

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

@PRANEETH-ALURU
Copy link

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

Copy link

ghost commented May 1, 2020

#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()

@SheCoderr
Copy link

Nice coding !

@FaridAbrari
Copy link

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
Copy link

alik604 commented Aug 24, 2020

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

@Satyaki44
Copy link

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..

@endolith
Copy link

endolith commented Feb 26, 2021

@unixsysdev
Copy link

@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()

download

@DanielPark827
Copy link

@keithdlandry
Can you explain Why sorting by highest amplitude is better than frequency?

@ChrisCalvin
Copy link

@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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment