iconEuler Home

Fast Fourier Transformation

This notebook introduce the Fast Fourier Transform (FFT). We will use
it to analyze a WAV file.

Mathematically, the FFT evaluates a polynomial in all complex n-th
roots of unity simultanuously. If n is a power of 2 (or has at least
not many prime large prime factors) the algorithm is rather quick. To
demonstrate this, we set xi equal to all 8-th roots of unity.
>reset; xi=exp(I*2*Pi/8); z=xi^(0:7);  ...
   plot2d(re(z),im(z),points=1):

15 - Introduction to Fast Fourier Transform

These are the complex values

15 - Introduction to Fast Fourier Transform

Now we generate a complex polynomial of degree 7. We take a random
polynomial.
>p=random(1,8)+1i*random(1,8);
We evaluate this polynomial at all points of z. Note
that a polynomial is represented in EULER in a vector
starting with the lowest coefficient.
>longestformat; evalpoly(z,p)
Column 1:
            3.61513775646999+4.027250420449434i 
Column 2:
       -0.06067706070439005+0.8797662696327804i 
Column 3:
         0.1325483481709888-0.1451728061492038i 
Column 4:
     -0.0005680160112071508-0.9346386176109833i 
Column 5:
           1.422171838062166-1.867849911920374i 
Column 6:
         0.3965191106395581+0.1482400937388588i 
Column 7:
        -0.4481910174229353-0.0126947295288784i 
Column 8:
         0.1863898275842544+0.1666181435935085i 
This is exactly the same as FFT of p.
>fft(p)
Column 1:
            3.61513775646999+4.027250420449434i 
Column 2:
       -0.06067706070438983+0.8797662696327807i 
Column 3:
         0.1325483481709888-0.1451728061492037i 
Column 4:
     -0.0005680160112072064-0.9346386176109834i 
Column 5:
           1.422171838062166-1.867849911920373i 
Column 6:
         0.3965191106395586+0.1482400937388583i 
Column 7:
       -0.4481910174229351-0.01269472952887946i 
Column 8:
         0.1863898275842536+0.1666181435935069i 
The inverse oparation interpolates values with a polynomial
in the roots of unity. We can go back and forth.
>fft(ifft(fft(p)))
Column 1:
            3.61513775646999+4.027250420449434i 
Column 2:
       -0.06067706070438991+0.8797662696327807i 
Column 3:
         0.1325483481709884-0.1451728061492038i 
Column 4:
     -0.0005680160112070398-0.9346386176109837i 
Column 5:
           1.422171838062166-1.867849911920373i 
Column 6:
         0.3965191106395585+0.1482400937388582i 
Column 7:
       -0.4481910174229351-0.01269472952887912i 
Column 8:
         0.1863898275842537+0.1666181435935071i 
The inverse FFT is almost the same as the original FFT.
>8*conj(ifft(conj(p)))
Column 1:
            3.61513775646999+4.027250420449434i 
Column 2:
       -0.06067706070438983+0.8797662696327807i 
Column 3:
         0.1325483481709888-0.1451728061492037i 
Column 4:
     -0.0005680160112072064-0.9346386176109834i 
Column 5:
           1.422171838062166-1.867849911920373i 
Column 6:
         0.3965191106395586+0.1482400937388583i 
Column 7:
       -0.4481910174229351-0.01269472952887946i 
Column 8:
         0.1863898275842536+0.1666181435935069i 
FFT can be used to evaluate a trigonometric series at equidistant
points. Let us take a cos series and evaluate it directly.
>t=0:Pi/128:2*Pi-Pi/128; s=1+2*cos(t)+cos(2*t); plot2d(t,s):

15 - Introduction to Fast Fourier Transform

Then we evaluate it with FFT. This is very much faster.
>a=[1,2,1]|zeros(1,253); totalmax(abs(re(fft(a))-s))
 4.973799150320701e-014 

Sound

IFFT (and FFT) can thus be used to determine frequencies in a
signal. The following signal consists of the frequencies 20 and
55.
>t=0:Pi/128:2*Pi-Pi/128; s=2*cos(20*t)+cos(55*t);
We can see the frequencies with IFFT.
>f=abs(ifft(s)); plot2d(0:127,f[1:128]):

15 - Introduction to Fast Fourier Transform

Adding randomly distributed noise does still shows the frequencies
clearly.
>s=s+normal(size(s)); f=abs(ifft(s)); plot2d(0:127,f[1:128]):

15 - Introduction to Fast Fourier Transform

Sound in Euler

At first we create 2 seconds of sound parameters.
>t=soundsec(2);
Then we generate the sound.
>s=sin(440*t)+sin(220*t)/2+sin(880*t)/2;
We can easily make a Fourier analysis of it.
>analyzesound(s):

15 - Introduction to Fast Fourier Transform

If you have a soundcard and speakers installed, you can hear the
sound file. Else you will get an error in the following command.
>playwave(s);
Now we add a lot of noise and listen to it.
>sa=s+normal(size(s)); playwave(sa);
We can also make a picture of the sound frequences versus the
length of the sound.
>mapsound(sa):

15 - Introduction to Fast Fourier Transform

Fourier Series

Finally we demonstrate Gibb's effect. Summing up sin(nt)/n for
odd n approaches the rectangle wave but with an error at 0 and
pi.
>t=linspace(0,pi,500)';
To compute the sum sin(k*t)/k of k=1,...,100, we use the matrix
language of Euler.
>s=sum(sin(t*(1:100))/(1:100))*2;
The plot shows that the Fourier sum converges to pi-x.
>plot2d(t',s'); plot2d("pi-x",color=3,add=1):

15 - Introduction to Fast Fourier Transform

But there is an overshoot close to 0, which is called Gibb's phenomenon.
Let us measure this overshoot for a very large n.
>t=linspace(0,pi,5000)'; max(sum(sin(t*(1:1000))/(1:1000))'*2-(pi-t'))
     0.5622809267763276 
Using a Riemann sum, one can see that the limit of this value
for n to infinity is the following.
>integrate("sinc(x)",0,pi)*2-pi
     0.5622814503751399 
Let us plot some the truncated series for n=1,3,...,61 using a
3D plot.
>n=1:2:61;
>t=linspace(0,pi,500)';
>S=cumsum(sin(t*n)/n);
>plot3d(n/10,t,S,hue=1,zoom=6):

15 - Introduction to Fast Fourier Transform

Euler Home