import numpy
import matplotlib.pyplot as plt

N=2048

def dft_forw(f):
    N = len(f)
    fhat = numpy.zeros(N, dtype=numpy.complex)
    fac = -2.0j*numpy.pi/N
    for n in range(N):
        for k in range(N):
            fhat[k] += f[n]*numpy.exp(fac*k*n)
    return fhat

def dft_back(fhat):
    N = len(fhat)
    f = numpy.zeros(N, dtype=numpy.complex)
    fac = 2.0j*numpy.pi/N
    for n in range(N):
        for k in range(N):
            f[n] += fhat[k]*numpy.exp(fac*k*n)
    f /= N
    return f

def rdft_forw(f):
    N = len(f)
    Nk = N/2+1
    fhat = numpy.zeros(Nk, dtype=numpy.complex)
    fac = -2.0j*numpy.pi/N
    for n in range(N):
        for k in range(Nk):
            fhat[k] += f[n]*numpy.exp(fac*k*n)
    return fhat

def rdft_back(fhat, N=None):
    Nk = len(fhat)
    if N is None:
        N = 2*(Nk-1)
    f = numpy.zeros(N, dtype=numpy.complex)
    fac = 2.0j*numpy.pi/N
    for n in range(N):
        f[n] = fhat[0] 
        for k in range(1,Nk):
            f[n] += 2.0*(fhat[k]*numpy.exp(fac*k*n)).real
    f /= N
    return f

def g(x):
    return x**-12 - x**-6

x = numpy.linspace(1, 5, N)

fhat = rdft_forw(g(x))

rec5 = rdft_back(fhat[:5],N)
rec10 = rdft_back(fhat[:10], N)
rec50 = rdft_back(fhat[:50],N)

plt.figure()
plt.plot(x, g(x),
         x, rec5,
         x, rec10,
         x, rec50,
)
plt.legend(('g(x)', '5', '10', '50'))

plt.show()