import numpy as np


def dft_coeff(x):
    """
    Evaluate c_n = 1/N*sum_{k=0}^{N-1}x_k exp(-i*2*pi*nk/N)
    :param x: array of samples of function
    """
    ns = np.size(x)
    c = np.zeros(ns, dtype=np.complex64)  # zero coefficients before summing
    for n in range(ns):
        for k in range(ns):
            arg = -2*np.pi*1j*n*k/ns
            c[n] = c[n]+x[k]*np.exp(arg)
        c[n] = c[n]/ns
    return c


def dft_synthesis(c, nc=0):
    """
    Evaluate the Fourier series as described above
    :param c: complex DFT coefficients
    :param nc: number of conjugated pairs of coefficients to be used
               (nc <= N/2-1 for even N or nc <= (N-1)/2 for odd N)
               (default: 0 indicating all coefficients)
    """
    ns = np.size(c)

    # initialize x_k with c_0
    x = np.ones(ns) * np.real(c[0])

    # distinguish even and odd N, we later sum up to nmax
    if ns % 2 == 0:
        nmax = ns // 2 - 1
    else:
        nmax = (ns - 1) // 2

    # if N is even and nc == 0, add coefficient c[N/2]*exp(i k pi)
    # the exponential is 1 for even k and -1 for odd k
    if ns % 2 == 0 and nc == 0:
        x = x + np.real(c[ns // 2]) * np.where(np.arange(0, ns) % 2 == 0, 1, -1)

    # check input nc, reset to nmax if greater
    if nc == 0: nc = nmax
    if nc > nmax: nc = nmax

    # do synthesis
    for n in range(1, nc + 1):
        for k in range(ns):
            arg = +2 * np.pi * 1j * n * k / ns
            x[k] = x[k] + 2 * np.real(c[n] * np.exp(arg))
    return x