import numpy as np
from scipy.fftpack import fft, ifft, rfft
from scipy.signal import fftconvolve
import matplotlib.pyplot as plt
#-----------------------------------------
 
# cyclic convolution
def convolute_cyc (f, g):
    Nf = len(f)
    Ng = len(g)
    Nc = Nf
    c = np.zeros_like(f, dtype= float)
 
    for n in xrange(Nc):
        for m in xrange(Nf):
            c[n] += f[m] * g[n-m]
    #end for
 
    return c
 
# acyclic convolution
def convolute_acyc (f, g):
    Nf = len(f)
    Ng = len(g)
    Nc = Nf + Ng - 1
    c = np.zeros(shape= (Nc), dtype= float)
 
    for n in xrange(Nc):
        for m in xrange(Nf):
            if m <= n and m > n - Ng:
                c[n] += f[m] * g[n-m]
            #endif
        #endfor
    #endfor
 
    return c
#----------------------------------------
# convolution example
 
T0 = 1.
a0 = 1.
 
def F (t):
    if np.abs(t) <= T0:
        return 1.
    #else:
    return 0.
 
def G (t):
    if t > 0:
        return np.exp(-a0 * t)
    # else:
    return 0.
 
# exact convolution
def C (t):
    if -T0 < t <= T0:
        return ( 1. - np.exp(-a0 * (t + T0)) ) / a0
    elif t > T0:
        return ( np.exp(a0 * T0) - np.exp(-a0 * T0) ) * np.exp(-a0*t) / a0
    # else:
    return 0.
 
LL = 6*T0  # boundary of x-axis
Npt = 2**10  # nr of points on x-axis
X = np.linspace(-LL,LL, Npt)  # x-axis
 
# arrays
f = np.asarray( map(F, X), dtype= float)
g = np.asarray( map(G, X), dtype= float)
c_exact = np.asarray( map(C, X), dtype= float)
 
# cyclic convolution
c_cyc = convolute_cyc(f, g)
 
# non-cyclic convolution
c_acyc = convolute_acyc(f, g)
 
# non-cyclic convolution by zero-padding
NPAD = Npt-1  # length of zero-padding
f_pad = np.concatenate( (f, [0.]*NPAD ) )  # (f | zero-pad of length NPAD)
g_pad = np.concatenate( (g, [0.]*NPAD ) ) #  (g | zero-pad of length NPAD)
c_fft = ifft( fft(f_pad) * fft(g_pad) )  # convolution using fftpack
 
#----------------------------------------
# plotting
 
YMAX = max(c_exact)
plt.plot(X, f, 'g-', lw= 1, label= "F")
plt.plot(X, g, 'r-', lw= 1, label= "G")
plt.plot(X, c_exact, 'b-', lw= 2, label= "exact conv.")
plt.plot((X[0], X[-1]), (YMAX, YMAX), 'k--') 
plt.ylim([0, 1.5*YMAX])
plt.legend()
plt.show()
 
plt.plot(c_fft.real, 'g-', lw= 1, label= "FFT conv.")
plt.plot(c_exact, 'b-', lw= 2, label= "exact conv.")
plt.legend()
plt.show()
 
plt.plot(c_acyc.real, 'g-', lw= 1, label= "conv. acyclic") 
plt.plot(c_exact, 'b-', lw= 2, label= "exact conv.")
plt.legend()
plt.show()

aW1wb3J0IG51bXB5IGFzIG5wCmZyb20gc2NpcHkuZmZ0cGFjayBpbXBvcnQgZmZ0LCBpZmZ0LCByZmZ0CmZyb20gc2NpcHkuc2lnbmFsIGltcG9ydCBmZnRjb252b2x2ZQppbXBvcnQgbWF0cGxvdGxpYi5weXBsb3QgYXMgcGx0CiMtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQoKIyBjeWNsaWMgY29udm9sdXRpb24KZGVmIGNvbnZvbHV0ZV9jeWMgKGYsIGcpOgogICAgTmYgPSBsZW4oZikKICAgIE5nID0gbGVuKGcpCiAgICBOYyA9IE5mCiAgICBjID0gbnAuemVyb3NfbGlrZShmLCBkdHlwZT0gZmxvYXQpCgogICAgZm9yIG4gaW4geHJhbmdlKE5jKToKICAgICAgICBmb3IgbSBpbiB4cmFuZ2UoTmYpOgogICAgICAgICAgICBjW25dICs9IGZbbV0gKiBnW24tbV0KICAgICNlbmQgZm9yCgogICAgcmV0dXJuIGMKCiMgYWN5Y2xpYyBjb252b2x1dGlvbgpkZWYgY29udm9sdXRlX2FjeWMgKGYsIGcpOgogICAgTmYgPSBsZW4oZikKICAgIE5nID0gbGVuKGcpCiAgICBOYyA9IE5mICsgTmcgLSAxCiAgICBjID0gbnAuemVyb3Moc2hhcGU9IChOYyksIGR0eXBlPSBmbG9hdCkKCiAgICBmb3IgbiBpbiB4cmFuZ2UoTmMpOgogICAgICAgIGZvciBtIGluIHhyYW5nZShOZik6CiAgICAgICAgICAgIGlmIG0gPD0gbiBhbmQgbSA+IG4gLSBOZzoKICAgICAgICAgICAgICAgIGNbbl0gKz0gZlttXSAqIGdbbi1tXQogICAgICAgICAgICAjZW5kaWYKICAgICAgICAjZW5kZm9yCiAgICAjZW5kZm9yCiAgICAKICAgIHJldHVybiBjCiMtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiMgY29udm9sdXRpb24gZXhhbXBsZQoKVDAgPSAxLgphMCA9IDEuCgpkZWYgRiAodCk6CiAgICBpZiBucC5hYnModCkgPD0gVDA6CiAgICAgICAgcmV0dXJuIDEuCiAgICAjZWxzZToKICAgIHJldHVybiAwLgoKZGVmIEcgKHQpOgogICAgaWYgdCA+IDA6CiAgICAgICAgcmV0dXJuIG5wLmV4cCgtYTAgKiB0KQogICAgIyBlbHNlOgogICAgcmV0dXJuIDAuCgojIGV4YWN0IGNvbnZvbHV0aW9uCmRlZiBDICh0KToKICAgIGlmIC1UMCA8IHQgPD0gVDA6CiAgICAgICAgcmV0dXJuICggMS4gLSBucC5leHAoLWEwICogKHQgKyBUMCkpICkgLyBhMAogICAgZWxpZiB0ID4gVDA6CiAgICAgICAgcmV0dXJuICggbnAuZXhwKGEwICogVDApIC0gbnAuZXhwKC1hMCAqIFQwKSApICogbnAuZXhwKC1hMCp0KSAvIGEwCiAgICAjIGVsc2U6CiAgICByZXR1cm4gMC4KICAgIApMTCA9IDYqVDAgICMgYm91bmRhcnkgb2YgeC1heGlzCk5wdCA9IDIqKjEwICAjIG5yIG9mIHBvaW50cyBvbiB4LWF4aXMKWCA9IG5wLmxpbnNwYWNlKC1MTCxMTCwgTnB0KSAgIyB4LWF4aXMKCiMgYXJyYXlzCmYgPSBucC5hc2FycmF5KCBtYXAoRiwgWCksIGR0eXBlPSBmbG9hdCkKZyA9IG5wLmFzYXJyYXkoIG1hcChHLCBYKSwgZHR5cGU9IGZsb2F0KQpjX2V4YWN0ID0gbnAuYXNhcnJheSggbWFwKEMsIFgpLCBkdHlwZT0gZmxvYXQpCgojIGN5Y2xpYyBjb252b2x1dGlvbgpjX2N5YyA9IGNvbnZvbHV0ZV9jeWMoZiwgZykKCiMgbm9uLWN5Y2xpYyBjb252b2x1dGlvbgpjX2FjeWMgPSBjb252b2x1dGVfYWN5YyhmLCBnKQoKIyBub24tY3ljbGljIGNvbnZvbHV0aW9uIGJ5IHplcm8tcGFkZGluZwpOUEFEID0gTnB0LTEgICMgbGVuZ3RoIG9mIHplcm8tcGFkZGluZwpmX3BhZCA9IG5wLmNvbmNhdGVuYXRlKCAoZiwgWzAuXSpOUEFEICkgKSAgIyAoZiB8IHplcm8tcGFkIG9mIGxlbmd0aCBOUEFEKQpnX3BhZCA9IG5wLmNvbmNhdGVuYXRlKCAoZywgWzAuXSpOUEFEICkgKSAjICAoZyB8IHplcm8tcGFkIG9mIGxlbmd0aCBOUEFEKQpjX2ZmdCA9IGlmZnQoIGZmdChmX3BhZCkgKiBmZnQoZ19wYWQpICkgICMgY29udm9sdXRpb24gdXNpbmcgZmZ0cGFjawoKIy0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KIyBwbG90dGluZwoKWU1BWCA9IG1heChjX2V4YWN0KQpwbHQucGxvdChYLCBmLCAnZy0nLCBsdz0gMSwgbGFiZWw9ICJGIikKcGx0LnBsb3QoWCwgZywgJ3ItJywgbHc9IDEsIGxhYmVsPSAiRyIpCnBsdC5wbG90KFgsIGNfZXhhY3QsICdiLScsIGx3PSAyLCBsYWJlbD0gImV4YWN0IGNvbnYuIikKcGx0LnBsb3QoKFhbMF0sIFhbLTFdKSwgKFlNQVgsIFlNQVgpLCAnay0tJykgCnBsdC55bGltKFswLCAxLjUqWU1BWF0pCnBsdC5sZWdlbmQoKQpwbHQuc2hvdygpCgpwbHQucGxvdChjX2ZmdC5yZWFsLCAnZy0nLCBsdz0gMSwgbGFiZWw9ICJGRlQgY29udi4iKQpwbHQucGxvdChjX2V4YWN0LCAnYi0nLCBsdz0gMiwgbGFiZWw9ICJleGFjdCBjb252LiIpCnBsdC5sZWdlbmQoKQpwbHQuc2hvdygpCgpwbHQucGxvdChjX2FjeWMucmVhbCwgJ2ctJywgbHc9IDEsIGxhYmVsPSAiY29udi4gYWN5Y2xpYyIpIApwbHQucGxvdChjX2V4YWN0LCAnYi0nLCBsdz0gMiwgbGFiZWw9ICJleGFjdCBjb252LiIpCnBsdC5sZWdlbmQoKQpwbHQuc2hvdygp