Convolutions and the Fast Fourier Transform

In mathematics and, in particular, functional analysis, convolution is a mathematical operation on two functions f and g, producing a third function. Using the fast Fourier transform can implement discrete convolutions in \( O(nlogn) \) time.

 

Preliminaries

  • Complex number and representation:

A complex number can be represented as \( c=a+bi \) where \( a,b\in \mathbb{R} \)

Euler’s formula shows that:

$$ e^{ix}=cos(x)+isin(x) $$

So, complex number can be also represented as:

$$ c=a+bi=re^{i\theta}\quad where\quad r=c\bar{c}=\sqrt{a^2+b^2}\quad \theta=tan^{-1}(\frac{b}{a}) $$

  • Root of unity:

An \( n^{th} \) root of unity, where n is a positive integer (i.e. \( n \) = 1, 2, 3, …), is a number \( z \) satisfying the equation:

$$ z^n=1 $$

\( z \) is assumed to be a complex number, let \( z=re^{i\theta} \) and substituted into the equation:

$$ r^ne^{in\theta}=r^n(cos(n\theta)+isin(n\theta))=1 $$

Hence:\( r=1 \) and \( \theta=\frac{2\pi}{n}k \)  where  \( k=0,1,2,…,n-1 \). So, let \( w_{n}=e^{\frac{2\pi}{n}} \) and \( z=w_{n}^{0},w_{n}^{1},w_{n}^{2},…,w_{n}^{n-1} \)

  • Property

Periodicity: \( w_{n}^{k}=w_{n}^{k+n} \) and \( w_{n}^{k}=-w_{n}^{k+\frac{n}{2}} \)

$$ w_{n}^{k+\frac{n}{2}}=w_{n}^{k}\times w_{n}^{\frac{n}{2}}=-w_{n}^{k}\quad (if\ n\ is \ even) $$

Summation: for any \( z=w_n^{i} \), \( 0 < i < n \)

$$ \sum_{i=0}^{n-1}{z^i}=0 $$

Since the equation \( z^n-1=0 \) can be factorized as:

$$ (z-1)(\sum_{i=0}^{n-1}{z^i})-1=\sum_{i=1}^{n}{z^i}-\sum_{i=1}^{n}{z^i}+z^n-1=z^n-1=0 $$

And \( z\ne 1 \) since \( z=w_n^{i} \), \( 0 < i < n \) and \( z \) is the root of \( z^n=1 \)

\( z-1\ne 1 \) So: \( \displaystyle \sum_{i=0}^{n-1}{z^i}=0 \)

 

Convolution definition

The convolution of two measurable function \( f \) and \( g \) is denoted as \( f\ast g \). It is defined as the integral of the product of the two functions where one is reversed and shifted. As such, it is a sort of integral transform:

$$ (f\ast g)(n)=\int_{-\infty}^{\infty} {f(\tau)g(n-\tau)d\tau} $$

Or:

$$ (f\ast g)(n)=\int_{-\infty}^{\infty} {f(n-\tau)g(\tau)d\tau} $$

As for the symbol t is used above, it can represent the different quantity according to the context such as representing  the time domain.

Discrete convolution

For two finite discrete series, the discrete convolution of \( f \) and \( g \) is given by:

$$ (f*g)[n]=\sum_{i=0}^{N-1}{f[i]g[n-i]} $$

When the elements of the sequences are the corresponding coefficients of two polynomials, then the coefficients of the ordinary product of the two polynomials are the convolution of the original two sequences called as the Cauchy product.

$$ \left( \sum_{i=0}^{n}{a_i} \right)\cdot \left( \sum_{j=0}^{n}{b_j} \right) =\sum_{k=0}^{n}{c_k} \quad where \quad c_k=\sum_{l=0}^{k}{a_{l}b_{k-l}} $$

In many cases, discrete convolutions can be converted to circular convolutions so that fast transforms with a convolution property can be used to implement the computation. The most common fast convolution algorithms use fast Fourier transform (FFT) algorithms via the circular convolution theorem.

 

Fast Fourier Transform  (FFT)

polynomial-multiplication

Calculating the discrete convolution of two digit sequences \( a\ast b \)is equivalent to multiply two polynomials which its coefficients are the elements of the digit sequences as mentioned above. The result is the coefficients sequences of the polynomial by multiplying two corresponding polynomials.

As for the multiplication of two polynomials, naive algorithm take \( O(n^2) \) time to executing by multiplying one polynomial with each items of another polynomial and adding \( n \) polynomials. However a polynomial of \( n \) degree only contain \( n+1 \) coefficients and the polynomial can be determined by solving \( (n+1)\times(n+1) \) linear equations with given \( n+1 \) points satisfying the polynomial.

To determined the multiplication \( c=a\times b \) of two polynomials \( a \) and \( c \) of degree \( n \), need to find at least \( 2n+1 \) points satisfying polynomial \( c \):

$$ \{x_i,y_i|c(x_i)=y_i\}\quad where \quad y_i=a(x_i)\times b(x_i)\quad since\quad c=a\times b $$

Now, the polynomials multiplication \( c=a\times b \) can be divided into three step(see the figure above), the technique known as the fast Fourier transform:

  1. For two polynomials \( a \) and \( b \) of degree \( n-1 \), transform them into the point-value representation(calculating the set contains \( 2n \) points \( \{x_i,y_i|a(x_i)=y_i\} \) and \( \{x_i,y_i|b(x_i)=y_i\} \)).
  2. Multiplying the point-value representation of two polynomials to construct the point-value representation of \( c \) (\( \{x_i,y_i|c(x_i)=y_i\}\;where\ y_i=a(x_i)\times b(x_i) \)).
  3. Interpolation, revert the polynomials \( c \) from the point-value representation.

Both step 1 and step 3 take \( O(nlogn) \) time and step 2  takes \( O(n) \) time.

The key for convert the polynomial to point-value representation and revert in \( O(nlogn) \) is substituting selected value into the polynomial. We will assume the degree \( n \) of polynomials \( a \) and \( b \) is a power of 2, if not, padding the higher degree items with 0 coefficient.

 

  • Convert to point-value representation

The \( 2n \) values to substitute into the polynomial are the root of unity mention in the preliminaries: \( w_{2n}^{0},w_{2n}^{1},w_{2n}^{2},…,w_{2n}^{2n-1} \).

For a polynomial:

\begin{split}f(x)=&\sum_{i=0}^{2n-1}{a_ix^i}\\=&\sum_{i=2k}^{2n-1}{a_ix^i}+\sum_{i=2k+1}^{2n-1}{a_ix^i}\\=&\sum_{k}^{n-1}{a_{2k}x^{2k}}+\sum_{k}^{n-1}{a_{2k+1}x^{2k+1}}\\=&\sum_{k}^{n-1}{a_{2k}{(x^2)}^k}+x\sum_{k}^{n-1}{a_{2k}{(x^2)}^k}\end{split}

So, we can decompose a polynomial to two sub-polynomial

define:

$$ f_{even}(x)=\sum_{k}^{n-1}{a_{2k}x^k} $$

$$ f_{odd}(x)=\sum_{k}^{n-1}{a_{2k+1}x^k} $$

And then:

$$ f(x)=f_{even}(x^2)+xf_{odd}(x^2) $$

Substitute \( w_{2n}^{0},w_{2n}^{1},w_{2n}^{2},…,w_{2n}^{2n-1} \) into the polynomial \( f \):

$$ f(w_{2n}^{i})=f_{even}(w_{2n}^{2i})+w_{2n}^{i}f_{odd}(w_{2n}^{2i}) $$

Notice: \( w_{2n}^{2i}=w_{n}^{i} \)

$$ f(w_{2n}^{i})=f_{even}(w_{n}^{i})+w_{2n}^{i}f_{odd}(w_{n}^{i}) $$

Both \( f_{even} \) and \( f_{odd} \) are polynomials contains \( n \) items, and the sequences of \( w_{2n}^{i} \) reduces to \( w_{n}^{i} \), which means the points need to calculate halved.

Restrict \( i<\frac{n}{2} \) so that we can reuse \( f_{even}(w_{n}^{i}) \) and \( f_{odd}(w_{n}^{i}) \) to calculate \( f(w_{2n}^{i}) \) and \( f(w_{2n}^{i+\frac{n}{2}}) \)

Because of the periodicity of the root of unity \( w_{n}^{k}=-w_{n}^{k+\frac{n}{2}} \)

$$ f(w_{2n}^{i})=f_{even}(w_{n}^{i})+w_{2n}^{i}f_{odd}(w_{n}^{i}) $$

$$ f(w_{2n}^{i+\frac{n}{2}})=f_{even}(w_{n}^{i})-w_{2n}^{i}f_{odd}(w_{n}^{i}) $$

For each recursive step: \( T(n)=2T(\frac{n}{2})+O(n) \)

Total time complexity: \( O(nlogn) \)

 

  • Interpolation, revert from point-value representation

Consider a polynomial \( \displaystyle c(x)=\sum_{i=0}^{2n-1}{c_ix^i} \) that we want to revert from its point-value representation \( \{x_i,y_i|c(x_i)=y_i\} \). Define a new polynomial \( d \):

$$ d(x)=\sum_{i=0}^{2n-1}{c(w_{2n}^i)x^i} $$

And then, substitute \( w_{2n}^{0},w_{2n}^{1},w_{2n}^{2},…,w_{2n}^{2n-1} \) into \( d \):

\begin{split}d(w_{2n}^{j})&=\sum_{i=0}^{2n-1}{c(w_{2n}^i)(w_{2n}^{j})^i}\\&=\sum_{i=0}^{2n-1}{(\sum_{k=0}^{2n-1}{c_kw_{2n}^{ik}})(w_{2n}^{j})^i}\\&=\sum_{k=0}^{2n-1}{c_k\sum_{i=0}^{2n-1}{w_{2n}^{ik}}w_{2n}^{ij}}\\&=\sum_{k=0}^{2n-1}{c_k\sum_{i=0}^{2n-1}{w_{2n}^{ik+ij}}}\\\end{split}

Because of summation of root of unity \( w_n^i \)

$$ \sum_{k=0}^{n-1}{(w_n^i)^k}=0\quad (for\ any\ w_n^i\ne1) $$

$$ \sum_{k=0}^{n-1}{(w_n^i)^k}=n\quad (for\ w_n^i=1) $$

Hence:

\begin{split}d(w_{2n}^{j})&=\sum_{k=0}^{2n-1}{c_k\sum_{i=0}^{2n-1}{w_{2n}^{i(k+j)}}}\\&=2nc_k \quad (for\ j+k=2n,\ then\ w_{2n}^{i(k+j)}=1)\end{split}

if \( j+k=2n \), then \( d(w_{2n}^{j})=2nc_k \),

So: \( \displaystyle c_k=\frac{d(w_{2n}^{2n-k})}{2n} \)

The inverse transform process of solving \( c_k \) is similar to step 1, converting the polynomial \( \displaystyle d(x)=\sum_{i=0}^{2n-1}{c(w_{2n}^i)x^i} \) to point-value representation at points \( w_{2n}^{0},w_{2n}^{1},w_{2n}^{2},…,w_{2n}^{2n-1} \). Therefore, this process can be also done in \( O(nlogn) \) time.

 

Implement

typedef std::vector<double> real_list;
typedef std::vector<std::complex<double>> complex_list;

complex_list epsilon;

inline
  size_t align_up(size_t n)
{
  size_t i=1;
  for (; i < n; i = i << 1);
  return i;
}

void init_epsilon(size_t n)
//initiate the list of root of unity
{
  if (epsilon.size() != n)
  {
    epsilon.clear();
    for (size_t i = 0; i != n; i++)
    {
      epsilon.emplace_back(cos(2.0 * M_PI * i / n), sin(2.0 * M_PI * i / n));
    }
  }
}

void __FFT(complex_list &coeff, complex_list &result, size_t n, size_t init_pos, size_t coeff_init_pos, size_t step)
{
  if (n == 1)
  {
    result[init_pos] = coeff[coeff_init_pos];
  }
  else
  {
    size_t m = n >> 1;
    __FFT(coeff, result, m, init_pos, coeff_init_pos, step << 1);
    __FFT(coeff, result, m, init_pos + m, coeff_init_pos + step, step << 1);
    complex_list __temp;
    __temp.resize(n);
    for (size_t i = 0; i < m; i++)
    {
      __temp[i] = result[init_pos + i] + epsilon[i*step] * result[init_pos + i + m];
      __temp[i + m] = result[init_pos + i] - epsilon[i*step] * result[init_pos + i + m];
    }
    for (size_t i = 0; i < n; i++)
    {
      result[init_pos + i] = __temp[i];
    }
  }
}

void FFT(complex_list &coeff, complex_list &result)
{
  size_t n = coeff.size();
  size_t align_n = align_up(n);
  for (; n < align_n; n++)
  {
    coeff.emplace_back(0, 0);
  }
  //Peddling
  init_epsilon(n);
  result.clear();
  result.resize(n);
  __FFT(coeff, result, n, 0, 0, 1);
}

void convolutions(complex_list &coeff_A, complex_list &coeff_B, complex_list &result_coeff)
{
  complex_list __FFT_A, __FFT_B, __FFT_C, __C;
  size_t n = max(coeff_A.size(), coeff_B.size());
  {
    size_t align_n = align_up(n) * 2;
    for (; n < align_n; n++)
    {
      coeff_A.emplace_back(0, 0);
      coeff_B.emplace_back(0, 0);
    }
    //Peddling
  }
  FFT(coeff_A, __FFT_A);
  FFT(coeff_B, __FFT_B);
  __C.resize(n);
  for (size_t i = 0; i < n; i++)
  {
    __C[i] = __FFT_A[i] * __FFT_B[i];
  }
  FFT(__C, __FFT_C);
  result_coeff.clear();
  result_coeff.resize(n);
  for (size_t i = 1; i < n; i++)
  {
    result_coeff[i-1] = __FFT_C[i] / static_cast<double>(n);
  }
  result_coeff[n - 1] = __FFT_C[0] / static_cast<double>(n);
  //Notice: w^0=w^n so the first coeff appears at last
}

 

Leave a Reply

Your email address will not be published. Required fields are marked *