% last modified 9 Feb. 1995

\include{psfig}
\documentstyle[times,12pt,twoside]{article}

\setlength{\textwidth}{6in}
\setlength{\textheight}{8.5in}
\setlength{\topmargin}{-.5in}	% top of page to top of running head minus 1in
% \setlength{\headsep}{.5in}	% bottom of running head to top of text
\setlength{\oddsidemargin}{.5in}	% left margin minus 1in on odd pages
\setlength{\evensidemargin}{0in}	% left margin minus 1in on even pages
\setlength{\parskip}{2ex}		% vertical space between paragraphs

% convolution symbol
\def\conv{\mathbin{\hbox{$\,*$\kern-1.5ex$\odot$}}}
%
% definite integrals and sums
\def\intinf{\int_{-\infty}^{+\infty}}
%
% matrix with rounded brackets
\def\rmatrix #1{\left\lgroup\matrix{ #1 }\right\rgroup}

\def\U{\hskip-.6em\relax&}
\def\ur{\cr\noalign{\vskip 1.5ex}}
\def\ud{\cdot}

% figure with caption and label
\def\fig #1 {       % #1=psfig args
    \begin{figure}[h]
        \centerline{\psfig{figure=#1}}
    \end{figure}
}


% NASTY STUFF: ==============================
\catcode`@=11	% make at signs act like letters, temporarily

% put back eqalign, which was removed by latex; it's better than eqnarray
%
\def\eqalign#1{\null\,\vcenter{\openup\jot\m@th
  \ialign{\strut\hfil$\displaystyle{##}$&$\displaystyle{{}##}$\hfil
      \crcr#1\crcr}}\,}

\catcode`@=12	% at signs are no longer letters
% END NASTY STUFF ==============================


\begin{document}
\noindent
Notes 7, Computer Graphics 2, 15-463

\vspace{.5in}
\begin{center}
{\noindent \bf \LARGE Discrete Fourier Transforms and the \\ \vskip 1ex
Fast Fourier Transform (FFT) Algorithm}

Paul Heckbert
\end{center}

We start in the continuous world; then we get discrete.

\section*{Definition of the Fourier Transform}

The Fourier transform (FT) of the function $f(x)$ is the function $F(\omega)$,
where:
$$
F(\omega) = \int_{-\infty}^{\infty} f(x) e^{-i \omega x} \,dx
$$
and the inverse Fourier transform is
$$
f(x) = {1 \over 2 \pi} \int_{-\infty}^{\infty}
    F(\omega) e^{i \omega x} \,d\omega
$$
Recall that $i=\sqrt{-1}$ and
$e^{i\theta}=\cos\theta + i\sin\theta$.

Think of it as a transformation into a different set of basis functions.
The Fourier transform uses complex exponentials (sinusoids)
of various frequencies as its basis functions.
(Other transforms, such as Z, Laplace, Cosine, Wavelet, and Hartley,
use different basis functions).

A Fourier transform pair is often written $f(x) \leftrightarrow F(\omega)$,
or ${\cal F}(f(x)) = F(\omega)$ where $\cal F$ is the Fourier transform
operator.

If $f(x)$ is thought of as a signal (i.e. input data) then we call $F(\omega)$
the signal's {\it spectrum}.
If $f$ is thought of as the {\it impulse response}
of a filter (which operates on input data to
produce output data) then we call $F$ the filter's {\it frequency response}.
(Occasionally the line between what's signal and what's filter becomes blurry).

\section*{Example of a Fourier Transform}

In signal processing terms,
the ideal low pass filter with cutoff frequency $\omega_c$ has a
frequency response that is a box:
$$
F(\omega) = \cases{
    1 & $|\omega| \le \omega_c$ \cr
    0 & $|\omega| > \omega_c$ \cr
}
$$
What is its impulse response?

We know that the impulse response is the inverse Fourier transform of
the frequency response, so
taking off our signal processing hat and putting on our mathematics hat,
all we need to do is evaluate:
$$
f(x) = {1 \over 2 \pi} \int_{-\infty}^{\infty}
    F(\omega) e^{i \omega x} \,d\omega
$$
for this particular $F(\omega)$:
$$
\eqalign{
f(x) &= {1 \over 2 \pi} \int_{-\omega_c}^{\omega_c} e^{i \omega x} \,d\omega \cr
     &= {1 \over 2 \pi} {e^{i \omega x} \over i x}
	\bigg|_{\omega=-\omega_c}^{\omega_c} \cr
     &= {1 \over \pi x} {e^{i \omega_c x} - e^{- i \omega_c x} \over 2 i} \cr
     &= {\sin \omega_c x \over \pi x}
	\qquad \qquad \hbox{since }
	\sin \theta = {e^{i \theta} - e^{-i \theta} \over 2 i} \cr
     &= {\omega_c \over \pi} {\rm sinc}({ \omega_c \over \pi } x) \cr
}
$$
where $\rm sinc (x) = \sin (\pi x)/(\pi x)$.
For antialiasing with unit-spaced samples, you want the cutoff frequency to
equal the Nyquist frequency, so $\omega_c = \pi$.

\section*{Fourier Transform Properties}

Rather than write ``the Fourier transform of an $X$ function is a $Y$
function'', we write the shorthand: $X \leftrightarrow Y$.
If $z$ is a complex number and $z=x+iy$ where $x$ and $y$ are
its real and imaginary parts,
then the complex conjugate of $z$ is $z^*=x-iy$.
A function $f(u)$ is {\it even} if $f(u)=f(-u)$,
it is {\it odd} if $f(u)=-f(-u)$,
it is conjugate symmetric if $f(u)=f^*(-u)$,
and it is conjugate antisymmetric if $f(u)=-f^*(-u)$.

\begin{tabular}{l}
    discrete $\leftrightarrow$ periodic \\
    periodic $\leftrightarrow$ discrete \\
    discrete, periodic $\leftrightarrow$ discrete, periodic \\
\hline
    real $\leftrightarrow$ conjugate symmetric \\
    imaginary $\leftrightarrow$ conjugate antisymmetric \\
\hline
    box $\leftrightarrow$ sinc \\
    sinc $\leftrightarrow$ box \\
    Gaussian $\leftrightarrow$ Gaussian \\
    impulse $\leftrightarrow$ constant \\
    impulse train $\leftrightarrow$ impulse train \\
\end{tabular}

{\it (can you prove the above?)}

When a signal is scaled up spatially, its spectrum is scaled down in
frequency, and vice versa: $f(ax) \leftrightarrow F(\omega/a)$
for any real, nonzero $a$.

\section*{Convolution Theorem}

The Fourier transform of a convolution of two signals is
the product of their Fourier transforms: $f \conv g \leftrightarrow FG$.
Recall that the convolution of signals $f$ and $g$ is
$(f \conv g)(x) = \intinf f(t) g(x-t) \,dt$.
So $\intinf f(t) g(x-t) \,dt \leftrightarrow F(\omega)G(\omega)$.

The Fourier transform of a product of two signals is
the convolution of their Fourier transforms:
$fg \leftrightarrow F \conv G / 2 \pi$.

\section*{Delta Functions}

The (Dirac) delta function $\delta(x)$ is defined such that
$\delta(x)=0$ for all $x \ne 0$, $\intinf \delta(t)\,dt = 1$, and
for any $f(x)$:
$$
(f \conv \delta)(x) = \intinf f(t) \delta(x-t) \,dt = f(x)
$$
The latter is called the {\it sifting property} of delta functions.
Because convolution with a delta is linear shift-invariant filtering,
translating the delta by $a$ will translate the output by $a$:
$$
\Bigl( f(x) \conv \delta(x-a) \Bigr)(x) = f(x-a)
$$

\section*{Discrete Fourier Transform (DFT)}

When a signal is discrete and periodic, we don't need the continuous Fourier
transform.
Instead we use the discrete Fourier transform, or DFT.
Suppose our signal is $a_n$ for $n=0 \ldots N-1$,
and $a_n=a_{n+jN}$ for all $n$ and $j$.
The spectrum of $a$ is:
\begin{equation}
A_k = \sum_{n=0}^{N-1} W_N^{kn} a_n
\label{dft.eq}
\end{equation}
where
$$
W_N = e^{-i {2 \pi \over N}}
$$
and $W_N^k$ for $k=0 \ldots N-1$ are called the {\it Nth roots of unity}.
They're called this because, in complex arithmetic,
$(W_N^k)^N = 1$ for all $k$.
They're vertices of a regular polygon inscribed in the unit
circle of the complex plane, with one vertex at $(1,0)$.

The sequence $A_k$ is the discrete Fourier transform of the sequence $a_n$.
Each is a sequence of $N$ complex numbers.


\subsection*{Two-point DFT (N=2)}

$W_2 = e^{-i \pi} = -1$,
and
$$
A_k = \sum_{n=0}^1 (-1)^{kn} a_n
= (-1)^{k \cdot 0} a_0 + (-1)^{k \cdot 1} a_1
= a_0 + (-1)^k a_1
$$
so
$$
\eqalign{
    A_0 &= a_0 + a_1 \cr
    A_1 &= a_0 - a_1 \cr
}
$$

\subsection*{Four-point DFT (N=4)}

$W_4 = e^{-i \pi/2} = -i$,
and
$$
A_k = \sum_{n=0}^3 (-i)^{kn} a_n
= a_0 + (-i)^k a_1 + (-i)^{2k} a_2 + (-i)^{3k} a_3
= a_0 + (-i)^k a_1 + (-1)^k a_2 + i^k a_3
$$
so
$$
\eqalign{
    A_0 &= a_0 + a_1 + a_2 + a_3 \cr
    A_1 &= a_0 - i a_1 - a_2 + i a_3 \cr
    A_2 &= a_0 - a_1 + a_2 - a_3 \cr
    A_3 &= a_0 + i a_1 - a_2 - i a_3 \cr
}
$$
This can also be written as a matrix multiply:
$$
\rmatrix{ A_0 \ur A_1 \ur A_2 \ur A_3 }
=
\rmatrix{
    1 &  1 &  1 &  1 \cr
    1 & -i & -1 &  i \cr
    1 & -1 &  1 & -1 \cr
    1 &  i & -1 & -i \cr
}
\rmatrix{ a_0 \ur a_1 \ur a_2 \ur a_3 }
$$
More on this later.

To compute $A$ quickly, we can pre-compute common subexpressions:
$$
\eqalign{
    A_0 &= (a_0 + a_2) + (a_1 + a_3) \cr
    A_1 &= (a_0 - a_2) - i (a_1 - a_3) \cr
    A_2 &= (a_0 + a_2) - (a_1 + a_3) \cr
    A_3 &= (a_0 - a_2) + i (a_1 - a_3) \cr
}
$$
This saves a lot of adds.
(Note that each add and multiply here is a complex (not real)
operation.)

If we use the following diagram for a complex
multiply and add:

\fig{mulad.eps,width=1.8in}

then we can diagram the 4-point DFT like so:

\fig{4pt.eps,width=2.5in}

If we carry on to $N=8$, $N=16$, and other power-of-two
discrete Fourier transforms, we get...

\section*{The Fast Fourier Transform (FFT) Algorithm}

The FFT is a fast algorithm for computing the DFT.
If we take the 2-point DFT and 4-point DFT and generalize them
to 8-point, 16-point, ..., $2^r$-point, we get the FFT algorithm.

To compute the DFT of an $N$-point sequence using equation (\ref{dft.eq})
would take $O(N^2)$ multiplies and adds.
The FFT algorithm computes the DFT using $O(N \log N)$ multiplies and adds.

There are many variants of the FFT algorithm.
We'll discuss one of them, the ``decimation-in-time'' FFT
algorithm for sequences whose length is
a power of two ($N=2^r$ for some integer $r$).

Below is a diagram of an 8-point FFT,
where $W=W_8=e^{-i \pi / 4} = (1-i)/ \sqrt{2}$:

\fig{8pt.eps,width=3in}

\section*{Butterflies and Bit-Reversal}

The FFT algorithm decomposes the DFT into $\log_2 N$ stages,
each of which consists of $N/2$ {\it butterfly} computations.
Each butterfly takes two complex numbers $p$ and $q$ and
computes from them two other numbers, $p+\alpha q$ and $p-\alpha q$,
where $\alpha$ is a complex number.
Below is a diagram of a butterfly operation:

\fig{butterfly.eps,width=1.5in}

In the diagram of the 8-point FFT above,
note that the inputs aren't in normal order:
$a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7$,
they're in the bizarre order:
$a_0, a_4, a_2, a_6, a_1, a_5, a_3, a_7$.
Why this sequence?

Below is a table of $j$ and the index of the $j$th input sample, $n_j$:
\begin{center}
\begin{tabular}{l|rrrrrrrr}
$j$ & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\
\hline
$n_j$ & 0 & 4 & 2 & 6 & 1 & 5 & 3 & 7 \\
\hline
$j$ base 2 & 000 & 001 & 010 & 011 & 100 & 101 & 110 & 111 \\
$n_j$ base 2 & 000 & 100 & 010 & 110 & 001 & 101 & 011 & 111 \\
\end{tabular}
\end{center}
The pattern is obvious if  $j$ and $n_j$ are written
in binary (last two rows of the table).
Observe that each $n_j$ is the {\it bit-reversal} of $j$.
The sequence is also related to breadth-first traversal of a binary tree.

It turns out that this FFT algorithm is simplest if the input array
is rearranged to be in bit-reversed order.
The re-ordering can be done in one pass through the array $a$:
\begin{verbatim}
    for j = 0 to N-1
        nj = bit_reverse(j)
        if (j<nj) swap a[j] and a[nj]
\end{verbatim}

\section*{FFT Explained Using Matrix Factorization}

The 8-point DFT can be written as a matrix product, where
we let $W=W_8=e^{-i \pi / 4} = (1-i)/ \sqrt{2}$:
{\scriptsize
$$
\rmatrix{
    A_0 \ur A_1 \ur A_2 \ur A_3 \ur A_4 \ur A_5 \ur A_6 \ur A_7 \ur
}
=
\rmatrix{
    W^0 \U W^0 \U W^0 \U W^0 \U W^0 \U W^0 \U W^0 \U W^0 \ur
    W^0 \U W^1 \U W^2 \U W^3 \U W^4 \U W^5 \U W^6 \U W^7 \ur
    W^0 \U W^2 \U W^4 \U W^6 \U W^0 \U W^2 \U W^4 \U W^6 \ur
    W^0 \U W^3 \U W^6 \U W^1 \U W^4 \U W^7 \U W^2 \U W^5 \ur
    W^0 \U W^4 \U W^0 \U W^4 \U W^0 \U W^4 \U W^0 \U W^4 \ur
    W^0 \U W^5 \U W^2 \U W^7 \U W^4 \U W^1 \U W^6 \U W^3 \ur
    W^0 \U W^6 \U W^4 \U W^2 \U W^0 \U W^6 \U W^4 \U W^2 \ur
    W^0 \U W^7 \U W^6 \U W^5 \U W^4 \U W^3 \U W^2 \U W^1 \ur
}
\rmatrix{
    a_0 \ur a_1 \ur a_2 \ur a_3 \ur a_4 \ur a_5 \ur a_6 \ur a_7 \ur
}
$$
}
Rearranging so that the input array $a$ is bit-reversed
and factoring the $8 \times 8$ matrix:
{\scriptsize
$$
\eqalign{
\rmatrix{
    A_0 \ur A_1 \ur A_2 \ur A_3 \ur A_4 \ur A_5 \ur A_6 \ur A_7 \ur
}
&=
\rmatrix{
    W^0 \U W^0 \U W^0 \U W^0 \U W^0 \U W^0 \U W^0 \U W^0 \ur
    W^0 \U W^4 \U W^2 \U W^6 \U W^1 \U W^5 \U W^3 \U W^7 \ur
    W^0 \U W^0 \U W^4 \U W^4 \U W^2 \U W^2 \U W^6 \U W^6 \ur
    W^0 \U W^4 \U W^6 \U W^2 \U W^3 \U W^7 \U W^1 \U W^5 \ur
    W^0 \U W^0 \U W^0 \U W^0 \U W^4 \U W^4 \U W^4 \U W^4 \ur
    W^0 \U W^4 \U W^2 \U W^6 \U W^5 \U W^1 \U W^7 \U W^3 \ur
    W^0 \U W^0 \U W^4 \U W^4 \U W^6 \U W^6 \U W^2 \U W^2 \ur
    W^0 \U W^4 \U W^6 \U W^2 \U W^7 \U W^3 \U W^5 \U W^1 \ur
}
\rmatrix{
    a_0 \ur a_4 \ur a_2 \ur a_6 \ur a_1 \ur a_5 \ur a_3 \ur a_7 \ur
}\cr\noalign{\vskip 3ex}
&=
\rmatrix{
      1 \U \ud \U \ud \U \ud \U W^0 \U \ud \U \ud \U \ud \ur
    \ud \U   1 \U \ud \U \ud \U \ud \U W^1 \U \ud \U \ud \ur
    \ud \U \ud \U   1 \U \ud \U \ud \U \ud \U W^2 \U \ud \ur
    \ud \U \ud \U \ud \U   1 \U \ud \U \ud \U \ud \U W^3 \ur
      1 \U \ud \U \ud \U \ud \U W^4 \U \ud \U \ud \U \ud \ur
    \ud \U   1 \U \ud \U \ud \U \ud \U W^5 \U \ud \U \ud \ur
    \ud \U \ud \U   1 \U \ud \U \ud \U \ud \U W^6 \U \ud \ur
    \ud \U \ud \U \ud \U   1 \U \ud \U \ud \U \ud \U W^7 \ur
}
\rmatrix{
      1 \U \ud \U W^0 \U \ud \U \ud \U \ud \U \ud \U \ud \ur
    \ud \U   1 \U \ud \U W^2 \U \ud \U \ud \U \ud \U \ud \ur
      1 \U \ud \U W^4 \U \ud \U \ud \U \ud \U \ud \U \ud \ur
    \ud \U   1 \U \ud \U W^6 \U \ud \U \ud \U \ud \U \ud \ur
    \ud \U \ud \U \ud \U \ud \U   1 \U \ud \U W^0 \U \ud \ur
    \ud \U \ud \U \ud \U \ud \U \ud \U   1 \U \ud \U W^2 \ur
    \ud \U \ud \U \ud \U \ud \U   1 \U \ud \U W^4 \U \ud \ur
    \ud \U \ud \U \ud \U \ud \U \ud \U   1 \U \ud \U W^6 \ur
}
\rmatrix{
      1 \U W^0 \U \ud \U \ud \U \ud \U \ud \U \ud \U \ud \ur
      1 \U W^4 \U \ud \U \ud \U \ud \U \ud \U \ud \U \ud \ur
    \ud \U \ud \U   1 \U W^0 \U \ud \U \ud \U \ud \U \ud \ur
    \ud \U \ud \U   1 \U W^4 \U \ud \U \ud \U \ud \U \ud \ur
    \ud \U \ud \U \ud \U \ud \U   1 \U W^0 \U \ud \U \ud \ur
    \ud \U \ud \U \ud \U \ud \U   1 \U W^4 \U \ud \U \ud \ur
    \ud \U \ud \U \ud \U \ud \U \ud \U \ud \U   1 \U W^0 \ur
    \ud \U \ud \U \ud \U \ud \U \ud \U \ud \U   1 \U W^4 \ur
}
\rmatrix{
    a_0 \ur a_4 \ur a_2 \ur a_6 \ur a_1 \ur a_5 \ur a_3 \ur a_7 \ur
}\cr
}
$$
}
where ``$\ud$'' means 0.

These are sparse matrices (lots of zeros), so multiplying by
the dense (no zeros) matrix on top is more expensive
than multiplying by the three sparse matrices on the bottom.

For $N=2^r$, the factorization would involve $r$ matrices
of size $N \times N$, each with 2 non-zero entries in each
row and column.

\section*{How Much Faster is the FFT?}

To compute the DFT of an $N$-point sequence using the definition,
$$
A_k = \sum_{n=0}^{N-1} W_N^{kn} a_n ,
$$
would require $N^2$ complex multiplies and adds,
which works out to $4N^2$ real multiplies and $4N^2$ real adds
(you can easily check this, using the definition of complex
multiplication).

The basic computational step of the FFT algorithm
is a butterfly.
Each butterfly computes two complex numbers of the form $p+\alpha q$
and $p-\alpha q$, so it requires one complex multiply ($\alpha \cdot q$)
and two complex adds.
This works out to 4 real multiplies and 6 real adds
per butterfly.

There are $N/2$ butterflies per stage, and $\log_2 N$ stages, so
that means about $4 \cdot N/2 \cdot \log_2 N = 2 N \log_2 N$
real multiplies and $3 N \log_2 N$ real adds for an $N$-point FFT.
(There are ways to optimize further, but this is the basic FFT algorithm.)

Cost comparison:

\begin{center}
\begin{tabular}{r|r|r|r|r}
& & {\bf BRUTE FORCE} & {\bf FFT} \\
$N$ & $r=\log_2 N$ & $4N^2$ & $2N \log_2 N$ & speedup \\
\hline
2 & 1 & 16 & 4 & 4 \\
4 & 2 & 64 & 16 & 4 \\
8 & 3 & 256 & 48 & 5 \\
\hline
1,024 & 10 & 4,194,304 & 20,480 & 205 \\
65,536 & 16 & $1.7 \cdot 10^{10}$ & $2.1 \cdot 10^6$ & \~{}$10^4$ \\
\end{tabular}
\end{center}

The FFT algorithm is a LOT faster for big $N$.

There are also FFT algorithms for $N$ not a power of two.
The algorithms are generally fastest when $N$ has many factors,
however.

An excellent book on the FFT is:
E. Oran Brigham,
{\it The Fast Fourier Transform},
Prentice-Hall,
Englewood Cliffs, NJ,
1974.

\section*{Why Would We Want to Compute Fourier Transforms, Anyway?}

The FFT algorithm is used for fast convolution
(linear, shift-invariant filtering).
If $h=f \conv g$ then
convolution of continuous signals involves an integral:\\
$h(x) = \intinf f(t) g(x-t) \,dt$,
but convolution of discrete signals involves a sum:
$h[x] = \sum_{t=-\infty}^{\infty} f[t] g[x-t]$.
We might think of $f$ as the signal and $g$ as the filter.

When working with finite sequences, the definition of convolution simplifies
if we assume that $f$ and $g$ have the same length $N$ and
we regard the signals as being periodic, so that $f$ and $g$ ``wrap around''.
Then we get {\it circular convolution}:
$$
h[x] = \sum_{t=0}^{N-1} f[t] g[x-t \hbox{ mod }N] \hbox{    for } x=0 \ldots N-1
$$

The convolution theorem says that the Fourier transform of the convolution
of two signals is the product of their Fourier transforms:
$f \conv g \leftrightarrow FG$.
The corresponding theorem for discrete signals is that
the DFT of the circular convolution of two signals is the product
of their DFT's.

Computing the convolution with a straightforward algorithm
would require $N^2$ (real) multiplies and adds --
too expensive!

We can do the same computation faster using discrete Fourier transforms.
If we compute the DFT of sequence $f$ and the DFT of sequence $g$,
multiply them point-by-point, and then
compute the inverse DFT, we'll get the same answer.
This is called {\it Fourier Convolution}:

\fig{fourconv.eps}

If we use the FFT algorithm, then the two DFT's and the one inverse
DFT have a total cost of $6N \log_2 N$ real multiplies,
and the multiplication of transforms in the frequency domain
has a negligible cost of $4N$ real multiplies.
The straightforward algorithm, on the other hand,
required $N^2$ real multiplies.

Fourier convolution wins big for large $N$.

Often, circular convolution isn't what you want,
but this algorithm can be modified to do standard ``linear'' convolution
by padding the sequences with zeros appropriately.

\section*{Fourier Transforms of Images}

To compute the Fourier transform of an image, you
\begin{itemize}
\item Compute DFT of each row, in place.
\item Compute DFT of each column, in place.
\end{itemize}

Practical issues:
You probably want to cyclically translate the picture so that pixel
(0,0), which now contains frequency $(\omega_x,\omega_y)=(0,0)$, moves
to the center of the image.
And you probably want to display pixel values proportional to
log(magnitude) of each complex number
(this looks more interesting than just magnitude).
For color images, do the above to each of the three channels
(R, G, and B) independently.

For an $N \times N$ picture, the cost is proportional to $N^2 \log N$.
The straightforward algorithm has a cost proportional to $N^4$.

FFT's are also used for synthesis of fractal textures and to create
images with a given spectrum.

\section*{Fourier Transforms and Arithmetic}

The FFT is also used for fast extended precision arithmetic,
and multiplication of high-degree polynomials,
since they also involve convolution.
If polynomials $p$ and $q$ have the form:
$p(x) = \sum_{n=0}^{N-1} f_n x^n$ and
$q(x) = \sum_{n=0}^{N-1} g_n x^n$ then their product is the polynomial
$$
\eqalign{
r(x) = p(x)q(x)
&=
    \Bigl( \sum_{n=0}^{N-1} f_n x^n \Bigr)
    \Bigl( \sum_{n=0}^{N-1} g_n x^n \Bigr) \cr
&=
    ( f_0 + f_1 x + f_2 x^2 + \cdots )
    ( g_0 + g_1 x + g_2 x^2 + \cdots ) \cr
&=
    f_0 g_0 + (f_0 g_1 + f_1 g_0) x +
    (f_0 g_2 + f_1 g_1 + f_2 g_0) x^2 + \cdots \cr
&=
    \sum_{n=0}^{2N-2} h_n x^n
}
$$
where $h_n = \sum_{j=0}^{N-1} f_j g_{n-j}$, and $h = f \conv g$.
Thus, computing the product of two polynomials involves the convolution
of their coefficient sequences.

Extended precision numbers (numbers with
hundreds or thousands of significant figures) are typically stored
in the form of a polynomial, where $x$ has a fixed value of
perhaps $2^{32}$ or $10^9$.
Multiplication of extended precision numbers involves the
convolution of high-degree polynomials.

When $N$ is small ($<100$, say), then straightforward convolution is
reasonable, but for large $N$, it makes sense to compute
convolutions using Fourier convolution.

\end{document}
