Filters in Image Processing Fourier Transform in one dimension David Svoboda email: svoboda@fi.muni.cz Centre for Biomedical Image Analysis Faculty of Informatics, Masaryk University, Brno, CZ September 20, 2019 David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 1 / 58 Outline 1 Discrete Fourier Transform 2 Continuous Fourier Transform 3 Fourier Transform Properties Visualization Common Properties Discrete specific 4 Fast Fourier Transform David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 2 / 58 Fourier Transform Jean Baptiste Joseph Fourier (1768–1830) David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 4 / 58 Fourier Transform Basis functions Common request: the basis should be orthonormal, i.e. ϕk · ϕl = δk,l ϕk(m) = 1 √ N e 2πimk N = 1 √ N cos 2πmk N + i sin 2πmk N the basis functions for N = 16: 0 10 −1 0 1 freq=0 0 10 −1 0 1 freq=1 0 10 −1 0 1 freq=2 0 10 −1 0 1 freq=3 0 10 −1 0 1 freq=4 0 10 −1 0 1 freq=5 0 10 −1 0 1 freq=6 0 10 −1 0 1 freq=7 0 10 −1 0 1 freq=8 0 10 −1 0 1 freq=9 0 10 −1 0 1 freq=10 0 10 −1 0 1 freq=11 0 10 −1 0 1 freq=12 0 10 −1 0 1 freq=13 0 10 −1 0 1 freq=14 0 10 −1 0 1 freq=15 0 5 10 15 −1 0 1 freq=0 0 5 10 15 −1 0 1 freq=1 0 5 10 15 −1 0 1 freq=2 0 5 10 15 −1 0 1 freq=3 0 5 10 15 −1 0 1 freq=4 0 5 10 15 −1 0 1 freq=5 0 5 10 15 −1 0 1 freq=6 0 5 10 15 −1 0 1 freq=7 0 5 10 15 −1 0 1 freq=8 0 5 10 15 −1 0 1 freq=9 0 5 10 15 −1 0 1 freq=10 0 5 10 15 −1 0 1 freq=11 0 5 10 15 −1 0 1 freq=12 0 5 10 15 −1 0 1 freq=13 0 5 10 15 −1 0 1 freq=14 0 5 10 15 −1 0 1 freq=15 David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 5 / 58 Fourier Transform Basis functions Properties periodical: ϕk+N (m) = 1 √ N e 2πim(k+N) N = 1 √ N e 2πimk N = ϕk (m) symmetrical: ϕN−k (m) = 1 √ N e 2πim(N−k) N = 1 √ N e− 2πimk N = 1 √ N e 2πi(−m)k N = ϕk (−m) orthonormal: ϕk (m) · ϕl (m) = N−1 m=0 1 √ N e 2πmki N 1 √ N e 2πmli N = 1 N N−1 m=0 e 2πm(k−l)i N = δk,l David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 6 / 58 Fourier Transform Definition Given 1D discrete function f of N samples and a basis (ϕk(m), k = {0, . . . , N − 1}), let us define: forward 1D discrete Fourier transform: F(k) ≡ f · ϕk = 1 √ N N−1 m=0 f (m)e−2πimk N inverse 1D discrete Fourier transform: f (m) ≡ F · ϕm = 1 √ N N−1 k=0 F(k)e 2πimk N f · ϕk = 1 √ N N−1 m=0 f (m)e 2πimk N = 1 √ N N−1 m=0 f (m)e−2πimk N David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 7 / 58 Fourier Transform Matrix notation If ϕk(m) = 1 √ N e −2πimk N = 1 √ N e−2πi N mk = 1 √ N ψmk then A = 1 √ N        1 1 1 . . . 1 1 ψ1 ψ2 . . . ψN−1 1 ψ2 ψ2·2 . . . ψ(N−1)2 1 ... ... ... ... 1 ψN−1 ψ2(N−1) . . . ψ(N−1)(N−1)        and F = Af ⇒ f = A−1 F = A T F Notice: ψ = e−2πi N is called the Nth root of unity. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 8 / 58 Fourier Transform The obvious question When we apply FT, we usually say ”let us decompose our signal into the sine waves . . . ” Why do we use another (so complicated) basis? Basis function is a ”sine wave” we avoid complex numbers more intuitive if basis function is a simple ”sine wave” sine waves without phase shift do not generate the whole space possible basis function: sin (km − α) α hidden in the sine function spoils the linearity; matrix multiplication cannot be used Basis function is ϕk (m) we have to use complex numbers ϕk (m) functions generate the whole space (form basis) this basis is orthonormal transform is linear, i.e. realized via matrix multiplication David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 9 / 58 Fourier Transform The meaning of Fourier coefficients If you perform inverse FT f (m) = F · ϕm = N−1 k=0 F(k)ϕk(m) = 1 √ N N−1 k=0 F(k)ei 2πkm N you literally compose the original signal f by combining the individual basis functions. Each basis function ϕk(m) = ei 2πmk N defines only the frequency. Fourier coefficient F(k) = |zk| (cos αk + i sin αk) = |zk|eiαk modifies the corresponding basis function ϕk(m) by scaling it and shifting it. F(k)ϕm(k) = |zk|eiαk ei 2πmk N = |zk|eiαk +i 2πmk N = |zk| cos 2πk N m + αk + i sin 2πk N m + αk David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 10 / 58 Fourier Transform The meaning of Fourier coefficients The number of samples is equal, i.e. |f | = |F| = |ϕk| −4 −3 −2 −1 0 1 2 3 4 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 m f(m) Original function −4 −3 −2 −1 0 1 2 3 4 −8 −6 −4 −2 0 2 4 6 8 k F(k) Real part David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 11 / 58 Fourier Transform The meaning of Fourier coefficients F(0) matches the lowest frequency in the signal f and corresponds to the “mean” of f : F(0) ≡ 1 √ N N−1 m=0 f (m)e−2πim0 N = 1 √ N N−1 m=0 f (m) F(0) is usually called DC term DC . . . “direct current” (current of zero frequency) F(1) . . . F(N − 1) are called AC terms AC . . . “alternating current” F(N−1 2 ) matches the highest frequency in the signal f Exercise: Why does the component F(N−1 2 ) correspond to the highest frequency? David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 12 / 58 Fourier Transform Continuous version We again deal with a projection the basis functions are: ϕω(x) = e2πixω = cos 2πxω + i sin 2πxω the projection of the function f onto the basis function ϕω(x) is the inner product as well as in the discrete case: f · ϕω = b a f (x)ϕω(x)dx David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 14 / 58 Fourier Transform Definition Given 1D integrable function f and a basis (ϕω, ω ∈ R), let us define: forward 1D continuous Fourier transform F(ω) ≡ ∞ −∞ f (x)e−2πixω dx inverse 1D continuous Fourier transform f (x) ≡ ∞ −∞ F(ω)e2πixω dω Notice: The sign flip explanation is beyond the scope of this course. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 15 / 58 Fourier spectrum visualization Centering Analysis of zero frequency Notice: Verify in FTutor1D (https://cbia.fi.muni.cz/software/ftutor1d.html). David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 17 / 58 Fourier spectrum visualization Real and Imaginary part Analysis of time shift David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 18 / 58 Fourier spectrum visualization Magnitude and phase Analysis of time shift David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 19 / 58 FT – Properties & Theorems Oddness and evenness 1 Each function f (x) is sum of its odd and even part: E(x) = 1 2 [f (x) + f (−x)] O(x) = 1 2 [f (x) − f (−x)] f (x) = E(x) + O(x) 2 sin (x) is an odd function cos (x) is an even function 3 Any FT basis function is composed of sine and cosine waves Corollary: FT of even function misses imaginary part (sine waves) FT of odd function misses real part (cosine waves) David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 20 / 58 FT – Properties & Theorems Oddness −4 −3 −2 −1 0 1 2 3 4 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 m f(m) Original function −4 −3 −2 −1 0 1 2 3 4 −8 −6 −4 −2 0 2 4 6 8 k F(k) Imaginary part original function (odd) its transform pair David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 21 / 58 FT – Properties & Theorems Evenness −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 m f(m) Original function −4 −3 −2 −1 0 1 2 3 4 −2 0 2 4 6 8 10 12 k F(k) Real part original function (even) its transform pair David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 22 / 58 FT – Properties & Theorems Transform pairs original function transform pair even even odd odd real and even real and even real and odd imaginary and odd real complex and hermitian Hermitian . . . conjugated symmetry, i.e. F(−x) = F(x) Exercise: How does change the FT if the original function is “real and even”? Prove your assertion. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 23 / 58 FT – Properties & Theorems Transform pairs If F(ω) = FT [f (x)] (ω) we usually write: f (x) ⊃ F(ω) and call f (x) and F(ω) the Fourier transform pair. Keep in mind that: not all the continuous functions are integrable, i.e. suitable for FT. For example sin x, cos x, f (x) = 1, any discrete function is suitable for FT, since it is finite, if f is even then |F(ω)| = |Re(F(ω))|, i.e. imaginary part is null, it is common to say that f (x) belongs to time domain. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 24 / 58 FT – Properties & Theorems Transform pairs Statement: Gaussian ⊃ Gaussian Proof: FT e−ax2 (ω) = ∞ −∞ e−ax2 e−2πixω dx = ∞ −∞ e−ax2 (cos 2πxω − i sin 2πxω) dx ... = π a e − π2 a ω2 David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 25 / 58 FT – Properties & Theorems Transform pairs Gaussian ⊃ Gaussian David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 26 / 58 FT – Properties & Theorems “sinc” function and its importance sinc(x) = sin x x David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 27 / 58 FT – Properties & Theorems “Π” function and its importance −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Π(x) =    1 |x| < 1 2 (1 2 |x| = 1 2) 0 |x| > 1 2 David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 28 / 58 FT – Properties & Theorems Transform pairs Statement: Rectangle Π function ⊃ Sinc Proof: FT [Π(ax)] (ω) = ∞ −∞ Π(ax)e−2πixω dx = 1 2a − 1 2a 1 · (cos 2πxω − i sin 2πxω) dx = · · · = 1 a sinc π a ω David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 29 / 58 FT – Properties & Theorems Transform pairs Rectangle Π function ⊃ Sinc David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 30 / 58 FT theorems Scaling Statement: f (ax) ⊃ 1 |a| F ω a Proof: FT [f (ax)] (ω) = ∞ −∞ f (ax)e−2πixω dx = 1 |a| ∞ −∞ f (ax)e−2πi(ax)(ω/a) d(ax) = 1 |a| F ω a Notice: Stretch in time domain corresponds to shrinkage in Fourier domain and vice versa. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 31 / 58 FT theorems Scaling Obstacles in discrete domain Continuous domain: scaling Discrete domain: upsampling/downsampling (loss of information) f (am) ≈ downsamplea−times(f (m)) Sampling may introduce new frequencies David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 32 / 58 FT theorems Scaling Fine sampling Rough sampling David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 33 / 58 FT theorems Scaling/Reciprocity David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 34 / 58 FT theorems Scaling/Reciprocity David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 35 / 58 FT theorems Scaling/Reciprocity David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 36 / 58 FT theorems Scaling/Reciprocity David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 37 / 58 FT theorems Shift Statement: f (x − a) ⊃ e−2πiaω F(ω) Proof: FT [f (x − a)] (ω) = ∞ −∞ f (x − a)e−2πixω dx = ∞ −∞ f (x − a)e−2πi(x−a)ω e−2πiaω d(x − a) = e−2πiaω F(ω) Notice: Shift affects only phase. The higher the frequency ω is the more the corresponding cosine wave is affected. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 38 / 58 FT theorems Shift David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 39 / 58 FT theorems Shift David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 40 / 58 FT theorems Linearity Statement: αf (x) + βg(x) ⊃ αF(ω) + βG(ω) Proof: FT [αf (x) + βg(x)] (ω) = ∞ −∞ [αf (x) + βg(x)] e−2πixω dx = α ∞ −∞ f (x)e−2πixω dx + β ∞ −∞ g(x)e−2πixω dx = αF(ω) + βG(ω) David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 41 / 58 FT theorems Linearity David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 42 / 58 FT theorems Convolution theorem Statement: f (x) ∗ g(x) ⊃ F(ω)G(ω) f (x)g(x) ⊃ F(ω) ∗ G(ω) Proof: FT [f (x) ∗ g(x)] (ω) = ∞ −∞   ∞ −∞ f (x )g(x − x )dx   e−2πixω dx = ∞ −∞ f (x )   ∞ −∞ g(x − x )e−2πixω dx   dx = ∞ −∞ f (x )e−2πix ω G(ω)dx = F(ω)G(ω) David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 43 / 58 FT theorems Convolution theorem The meaning: The convolution in time domain corresponds to point-wise multiplication in the Fourier domain, and vice versa. time domain domain Fourier = = f F G f*g FT FT IFT g FT(f)FT(g) faster? slower? * . David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 44 / 58 FT theorems Rayleigh’s energy theorem (Parseval’s theorem) Statement: ∞ −∞ |f (x)|2 dx = ∞ −∞ |F(ω)|2 dω Proof: ∞ −∞ |f (x)|2dx = ∞ −∞ f (x)f (x)dx = ∞ −∞ f (x)f (x)e−2πixω dx ω = 0 = F(ω ) ∗ F(−ω ) ω = 0 = ∞ −∞ F(ω )F(ω − ω )dω ω = 0 = ∞ −∞ F(ω)F(ω)dω David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 45 / 58 FT theorems Rayleigh’s energy theorem (Parseval’s theorem) Rayleigh’s energy theorem – the integral of the square of a function is equal to the integral of the square of its transform. Notice: |F(k)|2 is called a power spectrum. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 46 / 58 DFT theorems Convolution theorem Let f and g be 1D discrete periodic signals of length N, then: f ∗ g = IDFT [DFT(f ) · DFT(g)] √ N Proof: f (m) ∗ g(m) = N−1 k=0 f (k)g(m − k) = N−1 k=0 1 √ N N−1 n=0 F(n)e2πikn/N 1 √ N N−1 l=0 G(l)e2πi(m−k)l/N = 1 N N−1 n=0 F(n) N−1 l=0 G(l) N−1 k=0 e2πikn/N e2πi(m−k)l/N = 1 N N−1 n=0 F(n) N−1 l=0 G(l)e2πiml/N N−1 k=0 e2πik(n−l)/N = 1 N N−1 n=0 F(n) N−1 l=0 G(l)e2πiml/N δ(n − l)N = N−1 n=0 [F(n) · G(n)] e2πimn/N David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 47 / 58 DFT theorems Stretch If f (m) is a 1D function of length N, p ∈ N, and stretchp{f } = {g}, where g(n) = f (n/p) n = 0, p, 2p, . . . , (N − 1)p 0 otherwise then G(k) = 1 √ p    F(k) k = 0, . . . , N − 1 F(k − N) k = N, . . . , 2N − 1 ... ... F(k − (p − 1)N) k = (p − 1)N, . . . , pN − 1 Notice: Stretch by a factor p in the time domain results in p-fold repetition of F(k) in the frequency domain. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 48 / 58 DFT theorems Stretch An example of stretch for p = 2 0 1 2 3 4 5 6 7 0 2 4 6 8 10 m f(m) time domain 1 2 3 4 5 6 7 8 0 2 4 6 8 10 u |F(u)| Fourier (frequency) domain 0 5 10 15 0 2 4 6 8 10 m g(m) time domain 0 5 10 15 0 2 4 6 8 10 u |G(u)| Fourier (frequency) domain David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 49 / 58 DFT theorems Periodicity of DFT Statement: F(k + N) = F(k) Proof: F(k + N) = 1 √ N N−1 m=0 f (m)e− 2πim(k+N) N = 1 √ N N−1 m=0 f (m)e−2πimk N e−2πimN N = 1 √ N N−1 m=0 f (m)e−2πimk N = F(k) David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 50 / 58 DFT theorems Symmetry of real DFT Statement: F(−k) = F(k) Proof: F(−k) = 1 √ N N−1 m=0 f (m)e− 2πim(−k) N = 1 √ N N−1 m=0 f (m)e 2πimk N = F(k) iff f ∈ RN David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 51 / 58 Fast Fourier Transform Idea: N-point signal (N = 2m, m ∈ N) is decomposed into two N/2-point signals: one with all odd samples one with all even samples Example: 1 input signal: {0 1 2 3 4 5 6 7 } ⊃ ? 2 2× simpler DFT: {0 2 4 6} ⊃ {A B C D} 2× simpler DFT: {1 3 5 7} ⊃ {P Q R S} 3 stretch: {0 0 2 0 4 0 6 0} ⊃ 1√ 2 {A B C D A B C D} 4 stretch: {1 0 3 0 5 0 7 0 } ⊃ 1√ 2 {P Q R S P Q R S} shift: {0 1 0 3 0 5 0 7} ⊃ 1√ 2 {P ψQ ψ2R ψ3S ψ4P ψ5Q ψ6R ψ7S} 5 linearity: {0 1 2 3 4 5 6 7} = {0 0 2 0 4 0 6 0} + {0 1 0 3 0 5 0 7} ⇓ {0 1 2 3 4 5 6 7} ⊃ 1√ 2 {(A + P) (B + ψQ) (C + ψ2R) . . . } David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 53 / 58 Fast Fourier Transform Derivation: F(k) = 1 √ N N−1 m=0 f (m)e−2πimk N = 1 √ N N/2−1 m=0 f (2m)e− 2πi(2m)k N + 1 √ N N/2−1 m=0 f (2m + 1)e− 2πi(2m+1)k N = 1 √ N N/2−1 m=0 f (2m)e −2πimk N/2 + e−2πik N 1 √ N N/2−1 m=0 f (2m + 1)e −2πimk N/2 = Fe (k) + ψk Fo (k) Notice: ψ = e−2πi N David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 54 / 58 Fast Fourier Transform Idea: While it is possible repeat the division. F(k) → Fe (k), Fo (k) → Fee (k), Feo (k), Foe (k), Foo (k) → Feee (k), Feeo (k), Feoe (k), Feoo (k), Foee (k), Foeo (k), Fooe (k), Fooo (k) → . . . After log2(n) divisions we have Feeeoe...eoeooee(k) which is just one point long signal in Fourier domain. You should know that: {X} ⊃ {X} Exercise: What is the complexity of FFT? David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 55 / 58 Fast Fourier Transform One more illustration 150 1 2 3 4 5 6 7 10 11 12 13 1498 0 2 4 14121086 1 3 5 7 9 11 13 15 0 4 8 12 2 6 10 14 1 5 9 13 3 7 1511 0 8 4 12 2 10 6 14 1 9 5 13 11 7 153 0 8 4 12 2 10 6 14 1 9 5 13 3 11 15716 signals (1 point each) 2 signals (8 points each) 8 signals (2 points each) 4 signals (4 points each) 1 signal (16 points) David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 56 / 58 Bibliography Bracewell, R. N., Fourier transform and its applications / 2nd ed. New York: McGraw-Hill, 474 pages, ISBN 0070070156 Gonzalez, R. C., Woods, R. E., Digital image processing / 2nd ed., Upper Saddle River: Prentice Hall, 2002, pages 793, ISBN 0201180758 Veit, J., Integr´aln´ı transformace, Praha, 1983 Smith, Steven W. Digital signal processing: A practical guide for engineers and scientists; Amsterdam: Newnes, 2003, 650 pages; on-line version: http://www.dspguide.com/ David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 57 / 58 You should know the answers . . . Express the discrete Fourier transform as a matrix multiplication. Derive this matrix. How many Fourier basis functions do we need if we transform the signal of length N into the frequency domain? Formulate the forward discrete Fourier transform. Explain all the variables and constants. What does DC and AC terms mean? Explain the meaning of one particular Fourier coefficient in inverse Fourier transform. What is the product of the projection of an even function into a sin wave? Why are the wide functions in time domain transformed into their narrow counterparts in frequency domain and vice versa? Derive FFT for a signal of length 3m, m ∈ N. David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 58 / 58