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 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