Applications

1. Value of a polynomial
   Define a polynomial P(x) = $\sum_{{j=0}}^{{N-1}}$c_jx^j by the vector of its coefficients c : = [c₀, c₁,..c_N-1], where zeroes may be added so that N is a power of 2.
   • Compute the values of P(x) at
     x = a_k = $$\omega_{N}^{{-k}}$$ = exp($${\frac{{-2ik\pi}}{{N}}}$$),    k = 0..N - 1
     This is just the discrete Fourier transform of c since
     P(a_k) = $$\sum_{{j=0}}^{{N-1}}$$c_j($$\omega_{N}^{{-k}}$$)^j = F_N(c)_k
     Input, for example :
     P(x):=x+x^2; w:=i
     Here the coefficients of P are [0,1,1,0], N = 4 and $\omega$ = exp(2i$\pi$/4) = i.
     Input :
     fft([0,1,1,0])
     Output :
     [2,-1-i,0,-1+i]
     hence
     • P(0)=2,
     • P(-i)=P(w^-1)=-1-i,
     • P(-1)=P(w^-2)=0,
     • P(i)=P(w^-3)=-1+i.

   • Compute the values of P(x) at
     x = b_k = $$\omega_{N}^{{k}}$$ = exp($${\frac{{2ik\pi}}{{N}}}$$),    k = 0..N - 1
     This is N times the inverse fourier transform of c since
     P(a_k) = $$\sum_{{j=0}}^{{N-1}}$$c_j($$\omega_{N}^{{k}}$$)^j = NF_N^-1(c)_k
     Input, for example :
     P(x):=x+x^2 and w:=i
     Hence, the coefficients of P are [0,1,1,0], N = 4 and $\omega$ = exp(2i$\pi$/4) = i.
     Input :
     4*ifft([0,1,1,0])
     Output :
     [2,-1+i,0,-1-i]
     hence :
     • P(0)=2,
     • P(i)=P(w^1)=-1+i,
     • P(-1)=P(w^2)=0,
     • P(-i)=P(w^3)=-1-i.
     We find of course the same values as above...

2. Trigonometric interpolation
   Let f be periodic function of period 2$\pi$, assume that f (2k$\pi$/N) = f_k for k = 0..(N - 1). Find a trigonometric polynomial p that interpolates f at x_k = 2k$\pi$/N, that is find p_j, j = 0..N - 1 such that
   p(x) = $$\sum_{{j=0}}^{{N-1}}$$p_jexp(ijx),    p(x_k) = f_k
   Replacing x_k by it's value in p(x) we get:
   $$\sum_{{j=0}}^{{N-1}}$$p_jexp(i$${\frac{{j2k\pi}}{{N}}}$$) = f_k
   In other words, (f_k) is the inverse DFT of (p_k), hence
   (p_k) = $${\frac{{1}}{{N}}}$$F_N( (f_k) )
   If the function f is real, p_-k = $\overline{p}_{k}^{}$, hence depending whether N is even or odd:
   

   $$\Re \sum_{{k=0}}^{{\frac{N}{2}-1}} \Re {\frac{{N}}{{2}}} {\frac{{Nx}}{{2}}}$$ p(x) =

   $$\Re \sum_{{k=0}}^{{\frac{N-1}{2}}}$$ p(x) =

   
   

3. Fourier series
   Let f be a periodic function of period 2$\pi$, such that
   f (x_k) = y_k,    x_k = $${\frac{{2k\pi}}{{N}}}$$, k = 0..N - 1
   Suppose that the Fourier serie of f converges to f (this will be the case if for example f is continuous). If N is large, a good approximation of f will be given by:
   $$\sum_{{-\frac{N}{2} \leq n<\frac{N}{2}}}^{}$$c_nexp(inx)
   Hence we want a numeric approximation of
   c_n = $${\frac{{1}}{{2\pi}}}$$$$\int_{0}^{{2\pi}}$$f (t)exp(- int)dt
   The numeric value of the integral $\int_{0}^{{2\pi}}$f (t)exp(- int)dt may be computed by the trapezoidal rule (note that the Romberg algorithm would not work here, because the Euler Mac Laurin developpement has it's coefficients equal to zero, since the integrated function is periodic, hence all it's derivatives have the same value at 0 and at 2$\pi$). If $\tilde{{c_n}}$ is the numeric value of c_n obtained by the trapezoidal rule, then
   $$\tilde{{c_n}}$$ = $${\frac{{1}}{{2\pi}}}$$$${\frac{{2\pi}}{{N}}}$$$$\sum_{{k=0}}^{{N-1}}$$y_kexp(- 2i$${\frac{{nk\pi}}{{N}}}$$),     - $${\frac{{N}}{{2}}}$$ $$\leq$$ n < $${\frac{{N}}{{2}}}$$
   Indeed, since x_k = 2k$\pi$/N and f (x_k) = y_k:
   

   $${\frac{{nk\pi}}{{N}}}$$ f (x_k)exp(- inx_k) =

   $$\pi {\frac{{nN\pi}}{{N}}}$$ = y₀ = y_N

   
   Hence :
   [$$\tilde{{c}}_{0}^{}$$,..$$\tilde{{c}}_{{\frac{N}{2}-1}}^{}$$,$$\tilde{{c}}_{{\frac{-N}{2}}}^{}$$,..c_-1] = $${\frac{{1}}{{N}}}$$F_N([y₀, y₁...y_(N-1)])
   since
   • if n $\geq$ 0, $\tilde{{c}}_{n}^{}$ = y_n
   • if n < 0 $\tilde{{c}}_{n}^{}$ = y_n+N
   • $\omega_{N}^{}$ = exp(${\frac{{2i\pi}}{{N}}}$), then $\omega_{N}^{n}$ = $\omega_{N}^{{n+N}}$

   Properties

   • The coefficients of the trigonometric polynomial that interpolates f at x = 2k$\pi$/N are
     p_n = $$\tilde{{c}}_{n}^{}$$,     - $${\frac{{N}}{{2}}}$$ $$\leq$$ n < $${\frac{{N}}{{2}}}$$

   • If f is a trigonometric polynomial P of degree m $\leq$ ${\frac{{N}}{{2}}}$, then
     f (t) = P(t) = $$\sum_{{k=-m}}^{{m-1}}$$c_kexp(2ik$$\pi$$t)
     the trigonometric polynomial that interpolate f = P is P, the numeric approximation of the coefficients are in fact exact ( $\tilde{{c}}_{n}^{}$ = c_n).
   • More generally, we can compute $\tilde{{c}}_{n}^{}$ - c_n.
     Suppose that f is equal to it's Fourier serie, i.e. that :
     
     f (t) = $$\sum_{{m=-\infty}}^{{+\infty}}$$c_mexp(2i$$\pi$$mt),    $$\sum_{{m=-\infty}}^{{+\infty}}$$| c_m| < $$\infty$$
     Then :
     f (x_k) = f ($${\frac{{2k\pi}}{{N}}}$$) = y_k = $$\sum_{{m=-\infty}}^{{+\infty}}$$c_m$$\omega_{N}^{{km}}$$,    $$\tilde{{c_n}}$$ = $${\frac{{1}}{{N}}}$$$$\sum_{{k=0}}^{{N-1}}$$y_k$$\omega_{N}^{{-kn}}$$
     Replace y_k by it's value in $\tilde{{c_n}}$:
     $$\tilde{{c_n}}$$ = $${\frac{{1}}{{N}}}$$$$\sum_{{k=0}}^{{N-1}}$$$$\sum_{{m=-\infty}}^{{+\infty}}$$c_m$$\omega_{N}^{{km}}$$$$\omega_{N}^{{-kn}}$$
     If m $\neq$ n(mod N), $\omega_{N}^{{m-n}}$ is a N-th root of unity different from 1, hence:
     $$\omega_{N}^{{(m-n)N}}$$ = 1,    $$\sum_{{k=0}}^{{N-1}}$$$$\omega_{N}^{{(m-n)k}}$$ = 0
     Therefore, if m - n is a multiple of N ( m = n + l ^. N) then $\sum_{{k=0}}^{{N-1}}$$\omega_{N}^{{k(m-n)}}$ = N, otherwise $\sum_{{k=0}}^{{N-1}}$$\omega_{N}^{{k(m-n)}}$ = 0. By reversing the two sums, we get
     

     $$\tilde{{c_n}} {\frac{{1}}{{N}}} \sum_{{m=-\infty}}^{{+\infty}} \sum_{{k=0}}^{{N-1}} \omega_{N}^{{k(m-n)}}$$ =

     $$\sum_{{l=-\infty}}^{{+\infty}}$$ =

     = ...c_{n-2 ^. N} + c_n-N + c_n + c_n+N + c_{n+2 ^. N} + .....

     
     Conclusion: if | n| < N/2, $\tilde{{c_n}}$ - c_n is a sum of c_j of large indices (at least N/2 in absolute value), hence is small (depending on the rate of convergence of the Fourier series).
   Example, input
   f(t):=cos(t)+cos(2*t)
   x:=f(2*k*pi/8)$(k=0..7)
   Output :
   $\tt x=\{2,(\sqrt 2)/2,-1,-((\sqrt2)/2),0,-((\sqrt2)/2),-1,(\sqrt2)/2\}$
   fft(x)=[0.0,4.0,4.0,0.0,0.0,0.0,4.0,4.0]
   After a division by N = 8, we get
   c₀ = 0, c₁ = 4.0/8, c₂ = 4.0/2, c₃ = 0.0,
   c_-4 = 0, c_-3 = 0, c_-2 = 4.0/8, = c_-1 = 4.0/8
   Hence b_k = 0 and a_k = c_-k + c_k is equal to 1 if k = 1, 2 and 0 otherwise.

4. Convolution Product
   If P(x) = $\sum_{{j=0}}^{{n-1}}$a_jx^j and Q(x) = $\sum_{{j=0}}^{{m-1}}$b_jx^j are given by the vector of their coefficients a = [a₀, a₁,..a_n-1] and b = [b₀, b₁,..b_m-1], we may compute the product of these two polynomials using the DFT. The product of polynomials is the convolution product of the periodic sequence of their coefficients if the period is greater or equal to (n + m). Therefore we complete a (resp b) with m + p (resp n + p) zeros, where p is choosen such that N = n + m + p is a power of 2. If a = [a₀, a₁,..a_n-1, 0..0] and b = [b₀, b₁,..b_m-1, 0..0], then:
   P(x)Q(x) = $$\sum_{{j=0}}^{{n+m-1}}$$(a*b)_jx^j
   We compute F_N(a), F_N(b), then ab = F_N^-1(F_N(a) ^. F_N(b)) using the properties
   NF_N(x ^. y) = F_N(x)*F_N(y),    F_N(x*y) = F_N(x) ^. F_N(y)