Fast Fourier Transform

1
Fast Fourier Transform
? ? ???????????
2
Fast Fourier Transform

• The Fast Fourier Transform (FFT) is a very
  efficient algorithm for performing a discrete
  Fourier transform
• A.C.Clairaut,cosine-only finite Fourier
  series(1754),J.L.Lagrange,sine-only
  series(1762),Bernoulli,a series of sine and
  cosine.FFT principle first used by Gauss in 1805?
• FFT algorithm published by Cooley Tukey in 1965
• In 1969, the 2048 point analysis of a seismic
  trace took 13 ½ hours. Using the FFT, the same
  task on the same machine took 2.4 seconds!

3
Fourier ??
?????
Jean Baptiste Joseph Fourier (1768 - 1830)
???
??????????,????? ,exp(-i2pkx)
4
? g(x) ????
5
delta ??
TopHat ??
6
cosine ??
sine ??
7
????
????
a
-a
8
energy theorem, Rayleighs theorem
The zero frequency point
???
??
9
???Fourier??
10
???????
??????? (Continuous Fourier Transform)
??????? (Discrete Fourier Transform)
???????(Continuous Fourier Transform)
where
???????(Discrete Fourier Transform)
For u0,1,2,,N-1
For x0,1,2,,N-1
11
??Fourier??(Continuous Fourier Transform)
???
DFT IDFT
12
?????DFT
13
omega exp(-2pii/n) j 0n-1 k j' F
omega.(kj) an easier,and quicker, way to
generate F is F fft(eye(n))
14
??Fourier??(DFT) (????????MATLAB????)

• Requires N2 complex multiplies and N(N-1) complex
  additions

WNe-i 2p/N
??? ???
15
?????N/2?DFT??
16
N/2 DFT
x0,2,4,6
X07
N/2 DFT
x1,3,5,7

• Cross feed of Gk and Hk in flow diagram is
  called a butterfly, due to the shape

or simplify
-1
17
(No Transcript)
18
?N8,
bit reversal (0, 1, 2, 3, 4, 5, 6, 7 is
reordered to 0, 4, 2, 6, 1, 5, 3, 7)
Decimal Binary Binary Decimal
0 000 000 0
1 001 100 4
2 010 010 2
3 011 110 6
4 100 001 1
5 101 101 5
6 110 011 3
7 111 111 7
19
??WN/2 -1, Xk0 ? Xk1 ???? N/2,
Diagrammatically (butterfly),
There are N/2 butterflies for this stage of the
FFT, and each butterfly requires one
multiplication
The splitting of Xk into two half-size DFTs
can be repeated on Xk0 and Xk1 themselves,
20
The FFT for eight data values d1,d2,,d7 proceeds
as follows.








?? P55, C.W.Ueberhuber,Numerical Computation
2,Springer 1995.
21

• Xk00 is the N/4-point DFT of x0, x4,, xN-4,
• Xk01 is the N/4-point DFT of x2, x6,, xN-2,
• Xk10 is the N/4-point DFT of x1, x5,, xN-3,
• Xk11 is the N/4-point DFT of x3, x7,, xN-1.

22
?? c0 2 4 6 1 3 5 7
????????G.H.Golub,C.F.Van Loan??????
23
(No Transcript)
24
(No Transcript)
25
? n2t
function yfft(x,n) if n1 yx else mn/2
we-i2pn yTfft(x(02n),m) yBfft(x(12n),m)
d1,w,,wm-1T zd.yB yyTz yBz end
????????????,?? ??????
?? ??????(bit reversal permutation).
26
fftgui(y)??4?plots real(y), imag(y), real(fft(y)
), imag(fft(y))
print -deps FftGui.eps print depsc2 FftGui.eps
????????Moler??Numerical Computing with Matlab?
27
y01,y1yn-10,
28
y00,y11,y2yn-10,
29
FFT is the sum of two sinusoids
30
The Nyquist point
31
Some symmetries in the FFT. Ignoring the first
point in each plot, the real part is symmetric
about the Nyquist point and the imaginary part is
antisymmetric about the Nyquist point. More
precisely, ?y????n????,Yfft(y),? real(Y0)
?yj imag(Y0) 0 real(Yj) real(Yn-j),
j1,,n/2 imag(Yj) -imag(Yn-j),j1,,n/2
32
the sampling rate. Fs 32768 t
01/Fs0.25 the button in position (k,j) for
k14 for j13 y1
sin(2pifr(k)t) y2 sin(2pifc(j)t)
y (y1 y2)/2 input('Press any
key)') sound(y,Fs) end end
697 770 852 941
1209 1336 1477
33
(No Transcript)
34
For centuries, people have noted that the face of
the sun is not constant or uniform in appearance,
but that darker regions appear at random
locations on a cyclical basis. In1848,Rudolf
Wolfer proposed a rule that combined the number
and size of these sunspots into a single index.
load sunspot.dat t sunspot(,1)'
wolfer sunspot(,2)' n length(wolfer)
c polyfit(t,wolfer,1) trend
polyval(c,t) plot(t,wolfer trend,'-',
t,wolfer,'k.')
35
(No Transcript)
36
Now subtract off the linear trend and take the
FFT.
y wolfer - trend Y fft(y) Fs
1 Sample rate f (0n/2)Fs/n pow
abs(Y(1n/21)) plot(f f,0pow
pow,'c-', f,pow,'b.', ...
'linewidth',2, 'markersize',16)
37
(No Transcript)
38
(No Transcript)
39
plot(fft(eye(17))) axis square
40
Chebyshev Polynomial
????
??????
??????
???zgt1
???Chebyshev???
41
shg,hold on fplot('x',-1,1) fplot('2x2-1',-1,
1) fplot('4x3-3x',-1,1) fplot('8x4-8x21
',-1,1) fplot('16x5-20x35x',-1,1)
42
We do not make things, We make things better.