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 $n$ can be represented as a sequence of $n$ complex numbers. The most familiar is the coefficient representation, where $a_{0},\dots ,a_{n-1}$ represents $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 $n$ distinct points $w_{0},\dots ,w_{n-1}$, and record the of values $p(w_{0}),\dots ,p(w_{n-1})$. This is again a sequence of $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 $c$ is the coefficient form of a polynomial, and $v$ is point value form, then the Fourier transform takes $c$ to $v$, and its inverse takes $v$ to $c$.

For Fourier transforms, the points $w_{0},\dots ,w_{n-1}$ at which to evaluate are $\omega ^{0},\omega ^{1},\dots ,\omega ^{n-1}$, where $\omega$ is a primitive $n$th root of unity (so $\omega ^{n}=1$ and any $n$th root of unity is a power of $\omega$). There is some choice here: we use the engineering convention that $\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 $p$ is then $[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 ${\bar {\omega }}=e^{2\pi i/n}$. The inverse Fourier transform of $q$ is given by $(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