User:John Randall/FourierTransformAndPolynomialMultiplication

From J Wiki
Jump to: navigation, search

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 can be represented as a sequence of complex numbers. The most familiar is the coefficient representation, where represents . In signal processing terms, this corresponds to the time domain.

A less familiar representation is the point-value form. We choose distinct points , and record the of values . This is again a sequence of 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 is the coefficient form of a polynomial, and is point value form, then the Fourier transform takes to , and its inverse takes to .

For Fourier transforms, the points at which to evaluate are , where is a primitive th root of unity (so and any th root of unity is a power of ). There is some choice here: we use the engineering convention that . 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 is then .

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 . The inverse Fourier transform of is given by .

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