52 lines
1.5 KiB
Python
52 lines
1.5 KiB
Python
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 |