Lecture 16 Notes: Fast Convolution, Low Pass Filter Approximations, Integral Images (US 6,457,032)

6.801/6.866: Machine Vision, Lecture 16 1 Lecture 16: Fast Convolution, Low Pass Filter Approximations, Integral Images, (US 6,457,032) 1.1 Sampling and Aliasing Sampling is a ubiquitous operation for machine vision and general signal processing. Recall that PatMax, for instance, uses sampling in its methodology. PatMax also performs an interesting operation: low-pass filtering before sampling. Why is this performed? To answer this question, let us revisit the Nyquist Sampling Theorem. 1.1.1 Nyquist Sampling Theorem Before we dive into the mathematics behind this theorem, let us first build some intuition surrounding this theory. • If we can sample a signal at a high enough frequency, we can recover the signal exactly through reconstruction. • How is this reconstruction performed? We will convolve samples from the signal with sinc functions, and then superimpose these convolved results with one another. • It is hard to sample from a signal with infinite support. • What frequency do we need for this? Intuitively, to pick out how fast the signal needs to be moving, we certainly need to sample as quickly as the signal’s fastest-varying component itself. But do we need to sample even faster? It turns out the answer is yes. As we will see below: fmax < fsample =⇒ fsample > 2fmax 2 I.e. we will need to sample at more than twice the frequency of the highest-varying component of the signal. Let us look at this graphically. What happens if we sample at the frequency of the signal? 1 Cosine Function, cos(x) 2 True Function Interpolated Function f (x) 1 0 −1 −2 0 2 4 6 8 10 12 14 x Figure 1: Sampling only once per period provides us with a constant interpolated function, from which we cannot recover the original. Therefore, we must sample at a higher frequency. Note that this holds at points not on the peaks as well: Cosine Function, cos(x) 2 True Function Interpolated Function f (x) 1 0 −1 −2 0 2 4 6 8 10 12 14 x Figure 2: Sampling only once per period provides us with a constant interpolated function, from which we cannot recover the original. Therefore, we must sample at a higher frequency. What if we sample at twice the frequency? I.e. peaks and troughs: 2 Cosine Function, cos(x) 2 True Function Interpolated Function f (x) 1 0 −1 −2 0 2 4 6 8 10 12 14 x Figure 3: Sampling at twice the rate as the highest-varying component almost gets us there! This is known as the Nyquist Rate. It turns out we need to sample at frequencies that are strictly greater than this frequency to guarantee no aliasing - we will see why in the example below. Is this good enough? As it turns out, the inequality for Nyquist’s Sampling Theorem is there for a reason: we need to sample at greater than twice the frequency of the original signal in order to uniquely recover it: Cosine Function, cos(x) 2 True Function Interpolated Function f (x) 1 0 −1 −2 0 2 4 6 8 10 12 14 x Figure 4: It turns out we need to sample at frequencies that are strictly greater than this frequency to guarantee no aliasing we will see why in the example below. Therefore, any rate above 2 times the highest-varying frequency component of the signal will be sufficient to completely avoid aliasing. As a review, let us next discuss aliasing. 1.1.2 Aliasing Aliasing occurs when higher frequencies become indistinguishable from lower frequencies, and as a result they add interference and artifacts to the signal that are caused by sampling at too low of a frequency. Suppose we have a signal given by: f (x) = cos(2πf0 x) And suppose we sample this signal with frequency given by fs = 1δ . Then our sampled signal is given by: 1 f0 sk = cos(2πf0 k ) = cos(2π k) ∀ k ∈ {1, 2, ...} δ fs 3 Now let us consider what happens when we add multiples of 2π to this:  f  0 sk−2π = cos 2π k − 2πk fs  f   0 = cos 2π −1 k fs  f − f   0 s = cos 2π k fs  f − f   s 0 = cos 2π k , since cos(x) = cos(−x) ∀ x ∈ R fs Another way to put this - you cannot distinguish multiples of base frequencies with the base frequencies themselves if you sample at too low a frequency, i.e. below the Nyquist Rate. 1.1.3 How Can We Mitigate Aliasing? There are several strategies to mitigate aliasing effects: • (Note) Anti-aliasing measures must be taken before sampling. After sampling occurs, information is “lost”, so to speak, and the original signal cannot be recovered. • High frequency noise suppression with (approximate) low-pass filtering. As it turns out, exact lowpass filtering is impossible due to convergence properties of Fourier Series at cutoff frequencies (a phenomenon known as the Gibbs Phenomenon [1]. We will touch more on these anti-aliasing techniques and strategies throughout this lecture. 1.2 Integral Image Next, we will shift gears somewhat to discuss the notion of an integral image, and the critical role this technique plays in improving computational efficiency in image processing and machine vision. We will disucss this concept in both 1D and 2D. 1.2.1 Integral Images in 1D Block averaging is a common operation in computer vision in which we take the average over a set of values across an entire vector (1D) or matrix (2D), such as an image. This involves summing and then dividing by the total number of elements, which can become prohibitively computationally-expensive, particularly if this operation is being called many times and the averages that we need to compute are very large. Is there a more computationally-efficient way to do this? It turns out this computationally-simpler solution is through integral images. An integral image is essentially the sum of values from the first value to the ith value, i.e if gi defines the ith value in 1D, then: Δ Gi = i X gk ∀ i ∈ {1, · · · , K} k=1 Why is this useful? Well, rather than compute averages (normalized sums) by adding up all the pixels and then dividing, we simply need to perform a single subtraction between the integral image values (followed by a division by the number of elements we are averaging). For instance, if we wanted to calculate the average of values between i and j, then: j ḡ[i,j] = 1 X 1 (Gj − Gi ) gk = j−i j−i k=i This greatly reduces the amortized amount of computation, because these sums only need to be computed once, when we calculate the initial values for the integral image. 1.2.2 Integral Images in 2D Now, we can extend the notion of an integral image to 2D! Note further that integral “images” can extend beyond images! E.g. these can be done with gradients, Jacobians, Hessians, etc. One example in particular is calculating a Histogram of Gradients (HoG), which is quite useful for feature matching algorithms such as Scale-Invariant Feature Transform (SIFT). These approaches also map nicely onto GPUs, enabling for even faster computation times. 4 Let us now see how block averaging looks in 2D - in the diagram below, we can obtain a block average for a group of pixels in the 2D range (i, j) in x and (k, l) in y using the following formula: j ḡ([i,j],[k,l]) l XX 1 = gx,y (j − i)(l − k) x=i y=k But can we implement this more efficiently? We can use integral images again: Gi,j = j i X X gk,l k=1 l=1 Referencing the figure below, this becomes: Figure 5: Block averaging using integral images in 2D. As pointed out above, block averaging also extends beyond pixels! This can be computed for other measures such as gradients (e.g. Histogram of Gradients). Using the integral image values, the block average in the 2D range (i, j) in x and (k, l) in y becomes: ḡ([i,j],[k,l]) =   1 (Gj,l + Gi,k ) − (Gi,l + Gj,k ) (j − i)(l − k) Some comments about this: • This analysis can be extended to higher dimensions as well! Though the integral image will take longer to compute, and the equations for computing these block averages become less intuitive, this approach generalizes to arbitrary dimensions. • As we saw in the one-dimensional case, here we can also observe that after computing the integral image (a one-time operation that can be amortized), the computational cost for averaging each of these 2D blocks becomes independent of the size of the block being averaged. This stands in stark contrast to the naive implementation - here, the computational cost scales quadratically with the size of the block being averaged (or linearly in each dimension, if we take rectangular block averages). • Why is this relevant? Recall that block averaging implements approximate lowpass filtering, which can be used as a frequency suppression mechanism to avoid aliasing when filtering. • In other domains outside of image processing, the integral image is known as the “Summed-area Table” [2]. Since we intend to use this for approximate lowpass filtering, let us now change topics toward fourier analysis of this averaging mechanism to see how efficacious it is. 1.3 Fourier Analysis of Block Averaging Let us consider a one-dimensional “filter” that implements an approximate lowpass filtering for mechanism through block δ 1 averaging. Let us consider a function such that it takes 0 value outside of the range (− −δ 2 , 2 ), and value δ inside this range: ( 1 δ x ∈ (− −δ 2 , 2) h(x) = δ 0 o.w. Visually: 5 Block averaging h(x) for δ = 2 h(x) 0.6 0.4 0.2 0 −2 −1 0 x 1 2 Figure 6: Example h(x) for δ = 2. Let’s see what this Fourier Transform looks like. Recall that the Fourier Transform (up to a constant scale factor, which varies by domain) is given by: Z ∞ F (jω) = f (x)e−jωx dx −∞ Where jω corresponds to complex frequency. Substituting our expression into this transform: Z ∞ H(jω) = h(x)e−jωx dx −∞ Z δ 2 = − δ2 = = 1 1 −jωx x= δ2 [e ]x=− δ 2 δ jω e− jωδ 2 sin = 1 −jωx e dx δ −e −jωδ  δω 2 jωδ 2  (Sinc function) δω 2 Where in the last equality statement we use the identity given by: sin(x) = ejx − e−jx −2j Graphically, this sinc function appears (for δ = 2): 6  sin Sinc function H(jω) = δω 2  δω 2 H(jω) 1 0.5 0 −8 −6 −4 −2 0 2 4 6 8 jω Figure 7: Example H(jω) for δ = 2. This is the Fourier Transform of our block averaging “filter”. Although sinc functions in the frequency domain help to attenuate higher frequencies, they do not make the best lowpass filters. This is the case because: • Higher frequencies are not completely attenuated. • The first zero is not reached quickly enough. The first zero is given by: π π ω0 δ = =⇒ ω0 = 2 2 δ Intuitively, the best lowpass filters perfectly preserve all frequencies up to the cutoff frequencies, and perfectly attenuate everything outside of the passband. Visually:   sin Sinc function H(jω) = Sinc Filter Ideal lowpass filter 1 H(jω) δω 2 δω 2 0.5 0 −15 −10 −5 0 5 10 15 jω Figure 8: Frequency response comparison between our block averaging filter and an ideal lowpass filter. We also note that the “boxcar” function and the sinc function are Fourier Transform pairs! Although sinc functions in the frequency domain help to attenuate higher frequencies, they do not make the best lowpass filters. This is the case because: • Higher frequencies are not completely attenuated. 7 • The first zero is not reached quickly enough. The first zero is given by: ω0 δ π π = =⇒ ω0 = 2 2 δ Where else might we see this? It turns out cameras perform block average filtering becuase pixels have finite width over which to detect incident photons. But is this a sufficient approximate lowpass filtering technique? Unfortunately, oftentimes it is not. We will see below that we can improve with repeated block averaging. 1.4 Repeated Block Averaging One way we can improve our ability to attenuate higher frequencies - repeated block averaging! If our “boxcar” filter given above is given by b(x) (note that this was h(x) above), then our previous result y(x) can be written as: y1 (x) = f (x) ⊗ b(x) f (x) y1 (x) b(x) What happens if we add another filter? Then, we simply add another element to our convolution: y2 (x) = (f (x) ⊗ b(x)) ⊗ b(x) = y1 (x) ⊗ b(x) f (x) b(x) y1 (x) y2 (x) b(x) Adding this second filter is equivalent to convolving our signal with the convolution of two “boxcar” filters, which is a triangular filter: Triangular filter for δ = 2. 1 h(x) 0.8 0.6 0.4 0.2 0 −3 −2 −1 0 x 1 2 3 Figure 9: Example of a triangular filter resulting from the convolution of two “boxcar” filters. Additionally, note that since convolution is associative, for the “two-stage” approximate lowpass filtering approach above, we do not need to convolve our input f (x) with two “boxcar” filters - rather, we can convolve it directly with our trinagular filter b2 (x) = b(x) ⊗ b(x): y2 (x) = (f (x) ⊗ b(x)) ⊗ b(x) = f (x) ⊗ (b(x) ⊗ b(x)) = f (x) ⊗ b2 (x) Let us now take a brief aside to list out how discontinuities affect Fourier Transforms in the frequency domains: F • Delta Function: δ(x) ← →1 Intuition: Convolving a function with a delta function does not affect the transform, since this convolution simply produces the function. 8 F • Unit Step Function: u(x) ← → 1 jω Intuition: Convolving a function with a step function produces a degree of averaging, reducing the high frequency components and therefore weighting them less heavily in the transform domain. F • Ramp Function: r(x) ← → − ω12 Intuition: Convolving a function with a ramp function produces a degree of averaging, reducing the high frequency F d → jωF (jω) components and therefore weighting them less heavily in the transform domain. Derivative: dx f (x) ← Intuition: Since taking derivatives will increase the sharpness of our functions, and perhaps even create discontinuities, a derivative in the spatial domain corresponds to multiplying by jω in the frequency domain. As we can see from above, the more “averaging” effects we have, the more the high-frequency components of the signal will be filtered out. Conversely, when we take derivatives and create discontinuities in our spatial domain signal, this increases high frequency components of the signal because it introduces more variation. To understand how we can use repeated block averaging in the Fourier domain, please recall the following special properties of Fourier Transforms: 1. Convolution in the spatial domain corresponds to multiplication in the frequency domain, i.e. for all f (x), g(x), h(x) with corresponding Fourier Transforms F (jω), G(jω), H(jω), we have: F h(x) = f (x) ⊗ g(x) ← → H(jω) = F (jω)G(jω) 2. Multiplication in the spatial domain corresponds to convolution in the frequency domain, i.e. f (x), g(x), h(x) with corresponding Fourier Transforms F (jω), G(jω), H(jω), we have: for all F h(x) = f (x)g(x) ← → H(jω) = F (jω) ⊗ G(jω) For block averaging, we can use the first of these properties to understand what is happening in the frequency domain: F → Y (jω) = F (jω)(B(jω)2 ) y2 (x) = f (x) ⊗ (b(x) ⊗ b(x)) ← This operation is equivalent to a sinc2 function in the spatial domain:  Sinc function Squared H 2 (jω) = sin δω 2  δω 2 1 H 2 (jω) 0.8 0.6 0.4 0.2 0 −8 −6 −4 −2 0 2 4 6 8 jω Figure 10: Example H 2 (jω) for δ = 2. This is the Fourier Transform of our block averaging “filter” convolved with itself in the spatial domain. . 9  This is not perfect, but it is an improvement. In fact, the frequencies with this filter drop off with magnitude ω1 )2 . What happens if we continue to repeat this process with more block averaging filters? It turns out that for N “boxcar” filters that we  N use, the magnitude will drop off as ω1 . Note too, that we do not want to go “too far” in this direction, because this repeated block averaging process will also begin to attenuate frequencies in the passband of the signal. 1.4.1 Warping Effects and Numerical Fourier Transforms: FFT and DFT Two main types of numerical transforms we briefly discuss are the Discrete Fourier Transform (DFT) and Fast Fourier Transform (FFT). The FFT is an extension of the DFT that relies on using a “divide and conquer” approach to reduce the computational runtime from f (N ) ∈ O(N 2 ) to f (N ) ∈ O(N log N ) [3]. Mathematically, the DFT is given as a transform that transforms a sequence of N complex numbers {xn }N n=1 into another th sequence of N complex numbers {Xk }N [4]. The transform for the K value of this output sequence is given in closed-form k=1 as: Xk = N X 2π xn e−j N kn i=1 And the inverse transform for the nth value of this input sequence is given as: xn = N X 2π Xk ej N kn k=1 One aspect of these transforms to be especially mindful of is that they introduce a wrapping effect, since transform values are spread out over 2π intervals. This means that the waveforms produced by these transforms, in both the spatial (if we take the inverse transform) and frequency domains may be repeated - this repeating can introduce undesirable discontinuities, such as those seen in the graph below: Repeated Function x2 8 f (x) 6 4 2 0 0 2 4 x 6 8 Figure 11: Example of a repeated waveform that we encounter when looking at DFTs and FFTs. Fun fact: It used to be thought that natural images had a power spectrum (power in the frequency domain) that falls off as It turns out that this was actually caused by warping effects introduced by discrete transforms. 1 ω. This begs the question - how can we mitigate these warping effects? Some methods include: • Apodizing: This corresponds to multiplying your signal by a waveform, e.g. Hamming’s Window, which takes the form akin to a Gaussian, or an inverted cosine. • Mirroring: Another method to mitigate these warping effects is through waveform mirroring - this ensures continuity at points where discontinuties occurred: 10 Repeated Function x2 8 f (x) 6 4 2 0 0 2 4 x 6 8 Figure 12: Example of a mirrored waveform that we can use to counter and mitigate the discontinuity effects of warping from transforms such as the DFT and FFT. With this approach, the power spectrum of these signals falls off at 1 ω2 , rather than 1 ω. • Infinitely Wide Signal: Finally, a less practical, but conceptual helpful method is simply to take an “infinitely wide signal”. Let us now switch gears to talk more about the unit impulse and convolution. 1.5 Impulses and Convolution In this section, we will review impulse functions, and discuss how they relate to many properties with convolution. We will begin by reviewing delta functions and their properties. 1.5.1 Properties of Delta Functions Recall that the Dirac delta function δ(x), or impulse function, is defined according to the two properties: R∞ 1. Unit Area: −∞ δ(x)dx = 1 R∞ 2. Sifting Property: f (x0 ) = −∞ f (x)δ(x − x0 )dx Another way to conceptualize delta functions is through probabilistic distributions. We can use Gaussians (one of the only distributions to have a Fourier Transform). Recall that the (zero-mean) Gaussian Fourier Transform pair is given by: x ω2 σ2 1 F √ e− 2σ2 ← → e− 2 σ 2π An impulse can be conceptualized as the limit in which the variance of this Gaussian distribution σ 2 goes to 0, which corresponds to a Fourier Transform of 1 for all frequencies (which is the Fourier Transform of a delta function). Another way to consider impulses is that they are the limit of “boxcar” functions as their width goes to zero. Let us next generalize from a single impulse function to combinations of these functions. 1.5.2 Combinations of Impulses When are combinations of impulses helpful? it turns out that one particular combination can be used for approximating the derivative using our prior work on finite differences: h(x) = 1    δ(x + ) − δ(x − ) for some  > 0  2 2 11 Correlating (*note that this is not convolution - if we were to use convolution, this derivative would be flipped) this combination of impulse “filter” with an arbitrary function f (x), we compute a first-order approximation of the derivative: Z ∞ f 0 (x) ≈ f (x)h(x)dx −∞ Z ∞  1    = δ(x + ) − δ(x − ) dx 2 2 −∞  Therefore, combinations of impulses can be used to represent the same behavior as the “computational molecules” we identified before. It turns out that there is a close connection between linear, shift-invariant operators and derivative operators. With impulses motivated, let us now formally review convolution. 1.5.3 Convolution Review Recall from the example above that the convolution operation is simply the result of the correlation operation flipped. Mathematically, the convolution of functions f (x) and h(x) is given by the following commutative operation: g(x) = f (x) ⊗ h(x) Z ∞ = f (ξ)h(x − ξ)dξ −∞ = h(x) ⊗ g(x) Z ∞ = h(ξ)f (x − ξ)dξ −∞ 1.5.4 Analog Filtering with Birefringent Lenses Why filter in the analog domain? Returning to our original problem of mitigating aliasing artifacts through high-frequency suppression, unfortunately, if we wait to digitally filter out high frequencies after the image has already been taken by a camera or sensor, then we have already caused aliasing. One way to achieve this analog filtering is through Birefringent Lenses. Here, we essentially take two “shifted” images by convolving the image with a symmetric combination of offset delta functions, given mathematically by: h(x) = 1  1  δ(x + ) + δ(x − ) for some  > 0 2 2 2 2 Let us look at the Fourier Transform of this filter, noting the following Fourier Transform pair: F δ(x − x0 ) ← → e−jωx0 With this we can then express the Fourier Transform of this filter as: Z ∞ 1 F (jω) = √ h(x)e−jωx dx 2π −∞  jω 1  − jω = e 2 +e 2 2   ω = cos 2 With this framework, the first zero to appear here occurs at ω0 = high-frequency noise suppression. π . A few notes about these filters, and how they relate to • When these birefringent lenses are cascaded with a block averaging filter, this results in a combined filtering scheme in which the zeros of the frequency responses of these filters cancel out most of the high-frequency noise. • In the 2D case, we will have 2 birefringent filters, one for the x-direction and one for the y-direction. Physically, these are rotated 90 degrees off from one another, just as they are for a 2D cartesian coordinate system. 12 • High-performance lowpass filtering requires a large support (see definition of this below if needed) - the computational costs grow linearly with the size of the support in 1D, and quadratically with the size of the support in 2D. The support of a function is defined as the set where f (·) is nonzero [5]: supp(f ) = {x : f (x) 6= 0, x ∈ R} • Therefore, one way to reduce the computatonal costs of a filtering system is to reduce the size/cardinality of the support |supp(f )| - in some sense to encourage sparsity. Fortunately, this does not necessarily mean looking over a narrower range, but instead just considering less points overall. 1.5.5 Derivatives and Integrals as Convolution Operators and FT Pairs As we have seen before, convolving a function with a unit step function results in integrating the given function. Let us verify this below: Z ∞ f (x) ⊗ u(x) = f (ξ)u(x − ξ)dξ −∞ Z ∞ = f (x − ξ)u(ξ)dξ −∞ Z ∞ = f (ξ)dξ 0 Therefore, we can represent integral and derivative operators as Fourier Transform pairs too, denoted S for integration and D for derivative: F • S← → 1 jω F → jω • D← Niote that we can verify this by showing that convolving these filter operators corresponds to multiplying these transforms in frequency space, which results in no effect when cascaded together:  1  F F (f (x) ⊗ D) ⊗ S = f (x) ⊗ (D ⊗ S) ← → F (jω) jω = F (jω) ← → f (x) jω f (x) f (x) R∞ S D 0 f (ξ)dξ d dx f (x) D S f (x) f (x) Can we extend this to higher-order derivatives? It turns out we can. One example is the convolution of two derivative operators, which becomes:   F → H(jω) = D(jω)2 = (jω)2 = −ω 2 (Recall that j 2 = −1) h(x) = δ(x + ) − 2δ(x) + δ(x − ) = D ⊗ D ← 2 2 In general, this holds. Note that the number of integral operators S must be equal to the number of derivative operators D, e.g. for K order:      ⊗K ⊗K i=1 S ⊗ i=1 D ⊗ f (x) 1.5.6 Interpolation and Convolution A few notes about interpolation methods: • These methods work well computationally when there is sparsity in the operators we with - for instance, cubic or linear interpolation. Photoshop and other photo-editing software frameworks iuse this interpolation techniques. Performance deteriorates when we switch from cubic interpolation to Nearest Neighbor interpolation. • The bicubic spline is a popular interpolation technique. It was invented by IBM and was used in early satellite image processing. However, it requires many neighboring pixels, as it is composed of 7 points to interpolate over, so it is less computationally-efficient. 13 • Recall that one key element of computational efficiency we pursue is to use integral images for block averaging, which is much more efficient than computing naive sums, especially if (1) This block averaging procedure is repeated many times (the amortized cost of computing the integral image is lessened) and (2) This process is used in higher dimensions. • Linear interpolation can be conceptualized as connecting points together using straight lines between points. This corresponds to piecewise-linear segments, or, convolution with a triangle filter, which is simply the convolution of two “boxcar filters”: f (x) = f (1)x + f (0)(1 − x) Unfortunately, one “not-so-great” property of convolving with triangular filters for interpolation is that the noise in the interpolated result varies depending on how far away we are from the sampled noise. • Nearest Neighbor techniques can also be viewed through a convolutional lens - since this method produces piecewiseconstant interpolation, this is equivalent to convolving our sampled points with a “boxcar” filter! 1.5.7 Rotationally-Symmetric Lowpass Filter in 2D Where u and v are the spatial frequencies of x and y, i.e. the 2D Fourier Transform and Inverse Fourier Transform then take the forms of: Z ∞Z ∞ 1 F (u, v) = √ f (x, y)e−j(ux+vy) dxdy 2π −∞ −∞ Z ∞Z ∞ 1 f (x, y) = √ F (u, v)ej(ux+vy) dudv 2π −∞ −∞ The inverse transform of this can be thought of as a sinc function in polar coordinates: f (ρ, θ) = B 2 J1 (ρB) 2π ρB A few notes about this inverse transform function: • This is the point spread function of a microscope. • J1 (·) is a 1st-order Bessel function. • This relates to our defocusing problem that we encountered before. • In the case of defocusing, we can use the “symmetry” property of the Fourier Transform to deduce that if we have a circular point spread function resulting from defocusing of the lens, then we will have a Bessel function in the frequency/Fourier domain. • Though a pointspread function is a “pillbox” in the ideal case, in practice this is not perfect due to artifacts such as lens aberrations. 1.6 References 1. Gibbs Phenomenon, https://en.wikipedia.org/wiki/Gibbs phenomenon ____________________________________________ 2. Summed-area Table, https://en.wikipedia.org/wiki/Summed-area table ____________________________________________ 3. Fast Fourier Transform, https://en.wikipedia.org/wiki/Fast Fourier transform _______________________________________________ 4. Discrete Fourier Transform, https://en.wikipedia.org/wiki/Discrete Fourier transform __________________________________________________ 5. Support, https://en.wikipedia.org/wiki/Support (mathematics) ______________________________________________ 14