# User:John Randall/FourierTransformAndPolynomialMultiplication

While the discrete Fourier transform is best known from digital signal processing, it has a clear interpretation as a way of transforming between different representations of polynomials, and providing another way to do polynomial multiplication.

Complex polynomials of degree less than ${\displaystyle n}$ can be represented as a sequence of ${\displaystyle n}$ complex numbers. The most familiar is the coefficient representation, where ${\displaystyle a_{0},\dots ,a_{n-1}}$ represents ${\displaystyle p(z)=a_{0}+a_{1}z+\dots +a_{n-1}z^{n-1}}$. In signal processing terms, this corresponds to the time domain.

A less familiar representation is the point-value form. We choose ${\displaystyle n}$ distinct points ${\displaystyle w_{0},\dots ,w_{n-1}}$, and record the of values ${\displaystyle p(w_{0}),\dots ,p(w_{n-1})}$. This is again a sequence of ${\displaystyle n}$ complex numbers. It is not completely obvious that this contains the same information, but Lagrange interpolation shows us that it does. For a particular choice of points at which we evaluate, this corresponds to the frequency domain of signal processing.

If ${\displaystyle c}$ is the coefficient form of a polynomial, and ${\displaystyle v}$ is point value form, then the Fourier transform takes ${\displaystyle c}$ to ${\displaystyle v}$, and its inverse takes ${\displaystyle v}$ to ${\displaystyle c}$.

For Fourier transforms, the points ${\displaystyle w_{0},\dots ,w_{n-1}}$ at which to evaluate are ${\displaystyle \omega ^{0},\omega ^{1},\dots ,\omega ^{n-1}}$, where ${\displaystyle \omega }$ is a primitive ${\displaystyle n}$th root of unity (so ${\displaystyle \omega ^{n}=1}$ and any ${\displaystyle n}$th root of unity is a power of ${\displaystyle \omega }$). There is some choice here: we use the engineering convention that ${\displaystyle \omega =e^{-2\pi i/n}}$. The roots of unity below will be conjugated for the Fourier transform.

NB. Roots of unity
rou=:[: r. 2p1 * i. % ]


The Fourier transform of a polynomial ${\displaystyle p}$ is then ${\displaystyle [p(\omega ^{i})]}$.

NB. Fourier transform
F=:(] p. (+ @: rou @: #))"1


It would appear that the inverse Fourier transform (corresponding to interpolation) is more difficult than evaluation. The choice of roots of unity for the point-value form makes it expressible (up to a scale factor) as another evaluation, at powers of ${\displaystyle {\bar {\omega }}=e^{2\pi i/n}}$. The inverse Fourier transform of ${\displaystyle q}$ is given by ${\displaystyle (1/n)[q({\bar {\omega }}^{i})]}$.

NB. Inverse Fourier transform
require 'numeric'
FI =:(clean @: (# %~ ] p. (rou @: #)))"1

NB. Declaration of inverse
FT=:F :. FI


Here is a simple example of using filtering in the frequency domain. We create a 100-point signal c and transform it to v. The filter multiplies frequency 2 by 4, multiplies frequency 8 by 2, and cuts all others. This can be done by simple multiplication in the frequency domain. The inverse Fourier transform gives the modified signal.

require 'trig'
norm=:%: @: (+/ .*)~
n=:100
x=:2p1*(i.n)%n

NB. Original signal
c=:(2*cos 2*x)+(5*cos 6*x)+(4*sin 8*x)
NB. Transformed signal
v=:clean F c
NB. Filter
filter=:101 {. 0 0 4 0 0 0 0 0 2
filter=:}: (+|.)filter
NB. Filtered signal
f=:FI filter*v
NB. Expected result
c2=:(8*cos 2*x)+(8*sin 8*x)
NB. Check
norm f-c2
1.10942e_12


Returning to the interpretation of signals as polynomials, there is a direct connection between simple multiplication in the frequency domain and polynomial multiplication in the time domain. If  ppr is polynomial multiplication, and  p and  q are polynomials, then roughly speaking  (FT p ppr q)-:(FT p) * (FT q) . The only complication is that  p and  q must be padded with zeros to accommodate the size of the product.

NB. Auxiliary verbs for adding and removing zeros
pad=:,. 0 $~$
trim=:] #~ 0 +./\.@:~: ]

NB. Polynomial product using Fourier transform
fpr=:trim @: (*/&.: FT) @: pad @: ,:

NB. Regular polynomial product
ppr=: +//.@(*/)

NB. Comparison
a=:?. 20 # 10
b=:?. 100 # 20
a (norm @: (fpr-ppr)) b
1.19929e_10


We can also lift the filter to the time domain and use it as a convolution kernel, a fixed polynomial which multiplies the signal. Here the polynomial multiplication has to be interpreted mod the length of the signal.

NB. Convolution kernel
k=:FI  filter
NB. Convolution filtering
cf=:+/ (-#c)]\ k ppr c
norm cf-c2
5.01941e_13