# The rect function and its Fourier transform [[Fourier Series]] built periodic functions out of sines and cosines at the harmonics $n\omega_0$ of a fundamental. A function that never repeats has no fundamental, and the sum over harmonics becomes an integral over every angular frequency $\omega$: $ \hat f(\omega) = \int_{-\infty}^{\infty} f(x)\,e^{-i\omega x}\,dx, \qquad f(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty} \hat f(\omega)\,e^{i\omega x}\,d\omega . $ This is the **Fourier transform pair** in the convention used in lecture: angular frequency $\omega$ (radians per unit of $x$), no prefactor on the forward transform, and the $1/2\pi$ on the inverse. If $x$ is time in seconds, $\omega$ is in radians per second and $\omega/2\pi$ is the frequency in hertz. The one function whose transform everyone should know by heart is the rectangular pulse of half-width $L$ and height $A$, $ f(x) = \begin{cases} A, & |x| < L \\ 0, & |x| > L \end{cases} \qquad\Longrightarrow\qquad \hat f(\omega) = \int_{-L}^{L} A\,e^{-i\omega x}\,dx = A\,\frac{e^{-i\omega L} - e^{i\omega L}}{-i\omega} = \frac{2A\sin(\omega L)}{\omega}, $ so that, with $\operatorname{sinc}(z) = \sin(z)/z$, $ f(x) = A \text{ on } [-L, L] \quad\longleftrightarrow\quad \hat f(\omega) = 2AL\,\operatorname{sinc}(\omega L). $ Three things are visible in the formula before the picture is drawn. $\hat f(0) = 2AL$ is the area under the rect, which is what the $\omega = 0$ integral always computes. The zeros of $\hat f$ sit at $\omega = n\pi/L$, so a **narrower** rect (smaller $L$) has a **wider** sinc: the width $2L$ in $x$ times the distance $\pi/L$ to the first zero in $\omega$ is $2\pi$ whatever $L$ is, and that is the uncertainty principle in its most elementary form. And changing $A$ only rescales $\hat f$; it moves none of the zeros. Two more things about the graph. The formula forbids $\omega = 0$, but the limit exists ($\sin(\omega L)/(\omega L) \to 1$, by L'Hôpital), so the point is a *hole*, not a singularity, and we spackle it in with $\hat f(0) = 2AL$; the applet draws a small ring there. And the graph itself, in keeping with a course tradition, is **Squiddy**: the peak $2AL$ at the centre, oscillating arms decaying to either side, roots marching along at $\pi/L, 2\pi/L, 3\pi/L, \dots$, eyes included. Read Squiddy as instructions to the inverse transform: the pulse is *almost* a constant, so take a lot of the $\omega \approx 0$ non-wave; it is not completely constant, so add corrective waves in decreasing amounts, dipping negative in the first side lobes; and at the roots take none at all. Those are the frequencies the pulse simply does not need. The case worth singling out is **unit area**, $A = 1/(2L)$. Then $\hat f(\omega) = \operatorname{sinc}(\omega L)$ and $\hat f(0) = 1$ no matter what $L$ is. Send $L \to 0$: the pulse gets taller and narrower with its area pinned at one, which is the step-function construction of the **delta function**, and its transform flattens out to $\hat f(\omega) \equiv 1$, every frequency in equal measure. Send $L \to \infty$ instead and $f$ spreads into a vanishing constant while $\hat f$ collapses toward a spike at $\omega = 0$. The two limits are the two ends of the uncertainty principle, and the applet's *unit area* toggle lets you slide between them. ## Seeing it, and hearing it <div class="applet" data-applet="rect-transform"></div> For the sound, take $x$ to be time in **milliseconds**, so that $\omega$ is in **kiloradians per second** ($\omega/2\pi$ in kHz) and everything on the axes is in the audible range. Three buttons: - **pulse train.** A single rect has no pitch. Repeat it every $T$ ms and it does: a buzz at $f_0 = 1000/T$ Hz, angular fundamental $\omega_0 = 2\pi/T$. Tick *repeat every T* and the applet draws the repeats and, on the right, the Fourier coefficients of the train as stems at $\omega = n\omega_0$. This is the bridge between the last page and this one. The complex Fourier coefficients of a $T$-periodic function are $c_n = \frac{1}{T}\int_{-T/2}^{T/2} f(x)\,e^{-i n\omega_0 x}\,dx = \frac{1}{T}\,\hat f(n\omega_0),$ *the transform of one period, sampled at the harmonics.* The stems sit on Squiddy. As $T$ grows the stems crowd together (spacing $2\pi/T$) and trace out the whole curve; as $T \to \infty$ the series becomes the integral, and that limit *is* the Fourier transform. What you hear is the sinc envelope: the harmonics at $n f_0$ have exactly the amplitudes the stems show. - **one pulse.** The rect itself, $2L$ ms wide, played four times. It is a click, and its spectrum is the sinc, with nothing to hear at $\omega = n\pi/L$ krad/s. A wide pulse ($L = 3$) is a dull tick, because its spectrum is squeezed toward zero frequency; a narrow pulse ($L = 0.1$) is a sharp one, because its sinc is thirty times wider. - **the sinc as a sound.** Duality: with this convention, if $f \leftrightarrow \hat f$ then $\hat f(x) \leftrightarrow 2\pi f(-\omega)$, so a sinc in time has a **rectangular** spectrum. The button plays $\operatorname{sinc}(2\pi L x)$, whose spectrum is a perfectly flat band from $0$ out to $\omega = 2\pi L$ krad/s ($L$ kHz) and silence beyond. At $L = 0.1$ that is a low thump; at $L = 3$ it is a bright click. Same pair of functions, roles swapped. Things worth trying: - Set $L = 1$ and $T = 4$, duty cycle $\tfrac12$. The stems at $\omega = n\pi/2$ with $n = 2, 4, 6, \dots$ land exactly on the zeros $n\pi/L$ of the sinc: every even harmonic vanishes, and what you hear is the hollow, clarinet-like square wave of the previous page. Now set $L = 0.5$: every fourth harmonic is missing instead, and the buzz is brighter. The *zeros of the transform of one period* decide which harmonics the periodic signal is allowed to have. - Tick **unit area** and drag $L$ down toward $0.1$. The pulse shoots up, $\hat f(0)$ stays at $1$, and the sinc spreads until the first zero, at $\omega = 10\pi$, has left the window and the curve is nearly flat: that is $\delta(x) \leftrightarrow 1$ arriving. With the pulse train on, every stem has nearly the same height, so the click has every harmonic in equal strength and sounds as bright as a click can. Drag $L$ up to $3$ instead: the pulse fills half the window and $\hat f$ is a narrow spike with zeros every $\pi/3$. - Drag $T$ from $1$ to $6$ with $L$ fixed. The envelope does not move; only the sampling of it does. The pitch drops (from $1$ kHz to $167$ Hz) and the timbre stays put, which is the audible version of "$\hat f$ depends on the shape of one pulse, not on how often it repeats." - Tick **|f̂| as brightness**. The negative arms reflect up and a strip under the axis renders $|\hat f|$ as a grey level: blazing at the centre, dimming outward, dark points at the roots. That strip is one arm of the cross you will see in the two-dimensional pictures below, and it is how a spectrum is usually *displayed*, since brightness has no sign. - Compare the click and the thump at the same $L$. One is a rect in time with a sinc spectrum, the other a sinc in time with a rect spectrum; the applet's two panels are the same picture either way, with the axes relabelled. - Tick **draw my own f(x)** and draw a shape in the left panel (or start from a triangle, a Gaussian, two rects, a ramp). The transform is now computed numerically, $\hat f(\omega) \approx \sum_j f(x_j)\,e^{-i\omega x_j}\,\Delta x$, and unless your $f$ is even it is complex: the solid curve is $|\hat f|$ and the dashed one is $\operatorname{Re} \hat f$. A triangle is the rect convolved with itself, so its transform is $\operatorname{sinc}^2$: same zeros, no negative lobes, $1/\omega^2$ decay instead of $1/\omega$. A Gaussian's transform is a Gaussian ($e^{-2x^2} \leftrightarrow \sqrt{\pi/2}\,e^{-\omega^2/8}$). Two rects a distance $d$ apart give the single rect's sinc *modulated* by $\cos(\omega d/2)$, the same interference pattern as two slits. The pulse train and the click work on whatever you drew, so you can hear a shape as a timbre. ## Two dimensions: a bar and its cross For a function of two variables the transform integrates over both, $ \hat f(\omega_x, \omega_y) = \iint f(x, y)\,e^{-i(\omega_x x + \omega_y y)}\,dx\,dy, \qquad f(x, y) = \frac{1}{(2\pi)^2}\iint \hat f(\omega_x, \omega_y)\,e^{i(\omega_x x + \omega_y y)}\,d\omega_x\,d\omega_y, $ and when $f$ is a product of a function of $x$ and a function of $y$, so is $\hat f$. A rectangle of full width $a$ and full height $b$, $f = h$ for $|x| < a/2$, $|y| < b/2$, therefore has $ \hat f(\omega_x, \omega_y) = a\,b\,h\;\operatorname{sinc}\!\left(\tfrac{a\omega_x}{2}\right)\operatorname{sinc}\!\left(\tfrac{b\omega_y}{2}\right), $ a product of two sincs: a bright cross of ridges along the axes with dark lines every $2\pi/a$ in $\omega_x$ and every $2\pi/b$ in $\omega_y$. Rotate the coordinates by $\theta$, $ u = x\cos\theta + y\sin\theta, \qquad v = -x\sin\theta + y\cos\theta, $ and put a bar of (full) length $L$ along $u$ and (full) width $W$ along $v$, $f = h$ for $|u| < L/2$, $|v| < W/2$. Rotating $f$ rotates $\hat f$ by the same angle, so $ \hat f(\omega_x, \omega_y) = L\,W\,h\;\operatorname{sinc}\!\left(\tfrac{L\omega_u}{2}\right)\operatorname{sinc}\!\left(\tfrac{W\omega_v}{2}\right), \qquad \omega_u = \omega_x\cos\theta + \omega_y\sin\theta,\quad \omega_v = -\omega_x\sin\theta + \omega_y\cos\theta . $ (In this section $L$ and $W$ are the bar's full length and width, as in the course notebook, so the zeros sit every $2\pi/L$ and $2\pi/W$; in the 1-D section above $L$ was a half-width, which is why its zeros were at $\pi/L$.) <div class="applet" data-applet="rect-transform-2d"></div> The spectrum is drawn as $\log(1 + |\hat f|)$, because the central peak $\hat f(0,0) = LWh$ would otherwise swamp the side lobes. Things worth trying: - Start from the **square** preset and stretch $L$. The bar gets longer in $x$ and the pattern gets *narrower* in $\omega_x$: the long side of the bar is the short side of its transform, and vice versa. That is the 1-D uncertainty principle, once per axis. - Turn $\theta$. The whole pattern turns with it, so the bright streak in $\hat f$ always runs *perpendicular* to the bar. Rotation in space is rotation in frequency. - Set $L = W$ and any $\theta$: the square's cross, rotated. The pattern is not a picture of the object; it is a picture of which spatial frequencies the object contains, and in which directions. ## Discrete images Everything above is the continuous transform. A picture is an $n \times n$ array $f(m, n)$, and its transform is the two-dimensional **discrete Fourier transform** $ F(p, q) = \sum_{m=0}^{n-1}\sum_{m'=0}^{n-1} f(m, m')\,e^{-2\pi i (pm + qm')/n}, $ computed by the FFT ($\texttt{np.fft.fft2}$ in Python, $\texttt{Fourier}$ with $\texttt{FourierParameters} \to \{1, -1\}$ in Mathematica). The frequencies $p, q$ are integers, in cycles per image (the angular frequency of the $(p, q)$ term is $2\pi p/n$ radians per pixel across and $2\pi q/n$ down); the FFT returns them in the order $0, 1, \dots, n/2, -n/2+1, \dots, -1$, and the customary $\texttt{fftshift}$ rolls the array by $n/2$ so that $(0, 0)$ sits at the centre and the picture reads like the continuous one. <div class="applet" data-applet="image-dft"></div> The four test patterns are the ones in the course's `fourier_image.py`, and any picture from your own computer works too (it is converted to grayscale and resampled to $256 \times 256$ in the browser; nothing is uploaded anywhere). Things worth trying: - **square** and **rotated bar** reproduce the continuous pictures above, now with a pixel grid. The bar at $46 \times 6$ px and $45°$ is `test_images/bar_rotated.png`. - **disk.** Not separable, so no cross: the transform of a disk is the **Airy pattern**, $2J_1(r)/r$ in the radial coordinate, concentric rings. This is the diffraction pattern of a circular aperture, which is the physical reason telescopes and microscopes have a resolution limit. - **stripes.** A single cosine $\tfrac12 + \tfrac12\cos\bigl(2\pi(k_x j + k_y i)/n\bigr)$ is *one* pair of dots at $\pm(k_x, k_y)$, plus the spike at the origin from the mean. Tick *subtract the mean first* to remove the spike. Tighter stripes (more cycles) push the dots outward; tilting the stripes tilts the dot pair by the same angle. This is the whole idea of a spectrum in one picture: the position of a dot is a frequency and a direction, its brightness is an amplitude. - **draw a shape.** Paint with the mouse or a finger (brush size, paint or erase, or start from the current test image and add to it) and watch the spectrum follow. A single stroke in any direction makes a streak perpendicular to it; a blob makes rings; a row of dots makes a grating, which is a pair of dots in the spectrum. Draw a letter and see why each straight edge of it shows up as its own streak. - **a letter.** Type a letter (or a short word), choose a face, slant it. Every straight stroke is a streak in the spectrum, perpendicular to the stroke, so an **A** — a tripod — gives a three-way star that survives changes of face and slant, a **D** plants the clean stripe of its straight back, and an **S** has almost nothing straight in it and no streaks at all. Each letter carries a spectral fingerprint that varies less across handwriting than the letter itself does; that dual view, the letters *and* their transforms, is extra data a classifier already owns. - **come back.** The second row of panels is the way home. $F(p, q)$ is complex; the spectrum picture shows only its magnitude, and the panel on the lower left shows the part it hides, the **phase** $\arg F(p, q)$, which looks like noise and is anything but. Press **all of F** and the inverse transform returns the image to round-off, $10^{-16}$: the transform loses nothing. Press **magnitude only** and what comes back is a centred, symmetric blob no matter what you started with, because with every phase set to zero every wave is a cosine centred on the origin: the magnitudes say *how much* of each frequency, the phases say *where*. Press **phase only**, with every magnitude set to one, and the outline of the picture is still there: the edges live in the phase. This is why feeding the spectrum picture back into the transform does not give the picture back; it gives you the transform of $|F|$, which is (up to a flip) the autocorrelation of the image. Finally **low frequencies only** keeps the coefficients inside a radius you choose and throws the rest away: a percent of the coefficients holds most of the energy, and the picture comes back blurred but recognisable. That is compression, and the reason JPEG lives in $\hat f$-land. - **your own file.** Edges in the picture become streaks in the spectrum perpendicular to the edges; a picture with a strong texture shows the texture's frequency as a pair of dots; a blurry picture has a spectrum that dies quickly away from the centre. Look at a photograph of a brick wall, or a striped shirt. ## Where this goes The pulse-train applet is the whole chain the Day 16 notes describe, in one picture: **Fourier series** (stems at $n/T$) $\to$ **Fourier transform** ($T \to \infty$, the stems become the curve) $\to$ **sampling** (an $n$-point grid can only hold frequencies up to $n/2$ cycles per image, which is why the image spectra above stop where they do) $\to$ the **FFT** as the tool that computes it. When the input is a time series rather than a picture, the same $\texttt{fft}$ call and the same $\log(1 + |F|)$ display are how we find the annual and semi-annual cycles in the sea level record, and everything else that is hiding in it. ## The applets in code Each applet reduced to the few lines that compute what it shows, in the four languages of the course. Fully working versions of these — with the checks against closed forms, the figures, and the sound written to WAV files — are kept in the `fourier-codes` repository (`python/`, `matlab/`, `r/`, `mathematica/`; one file per applet, named after it). > [!example]- `rect-transform` — $f = A$ on $[-L, L] \leftrightarrow \hat f(\omega) = 2AL\,\mathrm{sinc}(\omega L)$, the pulse-train stems at $n\omega_0$, and the train as a tone > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > L, A, T = 1.0, 1.0, 4.0 # half-width, height, train period; unit area: A = 1 / (2 * L) > > F = lambda w: 2 * A * L * np.sinc(w * L / np.pi) # 2AL sinc(wL); np.sinc(z) = sin(pi z)/(pi z), so sinc(wL) = np.sinc(wL/pi) > > x = np.linspace(-6, 6, 1201); w = np.linspace(-30, 30, 1201); n = np.arange(-int(30 * T / (2 * np.pi)), int(30 * T / (2 * np.pi)) + 1) > > plt.subplot(1, 2, 1); plt.plot(x, A * (abs(x) < L)) > > plt.subplot(1, 2, 2); plt.plot(w, F(w), "r"); plt.stem(2 * np.pi * n / T, F(2 * np.pi * n / T)); plt.show() # stems: T c_n = F(n w0), w0 = 2 pi / T > > f0 = 1000 / T; ts = np.arange(44100 * 1.5) / 44100 # hear the train: x in ms, so f0 = 1000/T Hz > > nn = np.arange(1, int(22050 / f0) + 1) > > y = (2 * F(2 * np.pi * nn / T) / T) @ np.cos(2 * np.pi * f0 * np.outer(nn, ts)) > > from scipy.io import wavfile; wavfile.write("train.wav", 44100, (0.5 * y / abs(y).max() * 32767).astype(np.int16)) > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > L = 1; A = 1; T = 4; % half-width, height, train period; unit area: A = 1/(2*L) > > sinc_ = @(z) sin(z + (z == 0)) ./ (z + (z == 0)); % sin(z)/z with sinc(0) = 1, no toolbox needed > > F = @(w) 2 * A * L * sinc_(w * L); > > x = linspace(-6, 6, 1201); w = linspace(-30, 30, 1201); n = -floor(30*T/(2*pi)):floor(30*T/(2*pi)); > > subplot(1,2,1); plot(x, A * (abs(x) < L)) > > subplot(1,2,2); plot(w, F(w), 'r'); hold on; stem(2*pi*n/T, F(2*pi*n/T)) % stems: T c_n = F(n w0), w0 = 2 pi/T > > f0 = 1000/T; ts = (0:44100*1.5-1)/44100; nn = (1:floor(22050/f0))'; % hear the train: x in ms, so f0 = 1000/T Hz > > y = (2*F(2*pi*nn/T)/T)' * cos(2*pi*f0*nn*ts); sound(0.5*y/max(abs(y)), 44100) > > ``` > > > [!info]- R > > ```r > > L <- 1; A <- 1; T <- 4 # half-width, height, train period; unit area: A <- 1 / (2 * L) > > sinc <- function(z) ifelse(z == 0, 1, sin(z) / z); F <- function(w) 2 * A * L * sinc(w * L) > > x <- seq(-6, 6, length.out = 1201); w <- seq(-30, 30, length.out = 1201); n <- -floor(30 * T / (2 * pi)):floor(30 * T / (2 * pi)) > > par(mfrow = c(1, 2)); plot(x, A * (abs(x) < L), type = "l") > > plot(w, F(w), type = "l", col = "red"); segments(2 * pi * n / T, 0, 2 * pi * n / T, F(2 * pi * n / T)) # stems: T c_n = F(n w0) > > f0 <- 1000 / T; ts <- (0:(44100 * 1.5 - 1)) / 44100; nn <- 1:floor(22050 / f0) # hear the train: f0 = 1000/T Hz > > y <- colSums((2 * F(2 * pi * nn / T) / T) * cos(outer(nn, 2 * pi * f0 * ts))) > > tuneR::play(tuneR::Wave(round(32767 * 0.5 * y / max(abs(y))), samp.rate = 44100, bit = 16)) > > ``` > > > [!info]- Mathematica > > ```mathematica > > L = 1; A = 1; T = 4; F[w_] := 2 A L Sinc[w L]; (* unit area: A = 1/(2 L). F = FourierTransform[A UnitBox[x/(2 L)], x, w, FourierParameters -> {1, -1}] *) > > GraphicsRow[{Plot[A UnitBox[x/(2 L)], {x, -6, 6}, Exclusions -> None], > > Show[Plot[F[w], {w, -30, 30}, PlotStyle -> Red], > > ListPlot[Table[{2 Pi n/T, F[2 Pi n/T]}, {n, -Floor[30 T/(2 Pi)], Floor[30 T/(2 Pi)]}], Filling -> Axis]]}] (* stems: T c_n = F(n w0) *) > > f0 = 1000/T; (* hear the train: x in ms, so f0 = 1000/T Hz *) > > EmitSound @ Play[Evaluate @ Sum[2 F[2 Pi n/T]/T Cos[2 Pi f0 n t], {n, Floor[22050/f0]}], {t, 0, 1.5}] > > ``` > > [!example]- `rect-transform-2d` — the rotated bar and $\log(1 + |\hat f(\omega_x, \omega_y)|)$ from the closed form > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > L, W, th, h = 2.0, 0.4, np.deg2rad(30), 1.0 # full length, full width, angle, height > > x = np.linspace(-3, 3, 301); X, Y = np.meshgrid(x, x) > > U, V = X*np.cos(th) + Y*np.sin(th), -X*np.sin(th) + Y*np.cos(th) > > f = h * ((abs(U) <= L/2) & (abs(V) <= W/2)) > > w = np.linspace(-40, 40, 301); WX, WY = np.meshgrid(w, w) > > WU, WV = WX*np.cos(th) + WY*np.sin(th), -WX*np.sin(th) + WY*np.cos(th) > > F = L * W * h * np.sinc(L * WU / (2 * np.pi)) * np.sinc(W * WV / (2 * np.pi)) # L W h sinc(L wu/2) sinc(W wv/2); np.sinc(z) = sin(pi z)/(pi z) > > plt.subplot(1, 2, 1); plt.imshow(f, extent=[-3, 3, -3, 3], origin="lower", cmap="gray", vmin=0, vmax=2) > > plt.subplot(1, 2, 2); plt.imshow(np.log1p(abs(F)), extent=[-40, 40, -40, 40], origin="lower", cmap="gray"); plt.show() > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > L = 2; W = 0.4; th = 30*pi/180; h = 1; % full length, full width, angle, height > > sinc_ = @(z) sin(z + (z == 0)) ./ (z + (z == 0)); > > x = linspace(-3, 3, 301); [X, Y] = meshgrid(x, x); U = X*cos(th) + Y*sin(th); V = -X*sin(th) + Y*cos(th); > > f = h * (abs(U) <= L/2 & abs(V) <= W/2); > > w = linspace(-40, 40, 301); [WX, WY] = meshgrid(w, w); WU = WX*cos(th) + WY*sin(th); WV = -WX*sin(th) + WY*cos(th); > > F = L * W * h * sinc_(L*WU/2) .* sinc_(W*WV/2); > > colormap gray; subplot(1,2,1); imagesc(x, x, f, [0 2]); axis xy square > > subplot(1,2,2); imagesc(w, w, log1p(abs(F))); axis xy square > > ``` > > > [!info]- R > > ```r > > L <- 2; W <- 0.4; th <- 30 * pi / 180; h <- 1 # full length, full width, angle, height > > sinc <- function(z) ifelse(z == 0, 1, sin(z) / z) > > x <- seq(-3, 3, length.out = 301); w <- seq(-40, 40, length.out = 301) > > f <- outer(x, x, function(x, y) h * (abs(x*cos(th) + y*sin(th)) <= L/2 & abs(-x*sin(th) + y*cos(th)) <= W/2)) > > F <- outer(w, w, function(wx, wy) L * W * h * sinc(L*(wx*cos(th) + wy*sin(th))/2) * sinc(W*(-wx*sin(th) + wy*cos(th))/2)) > > par(mfrow = c(1, 2)); image(x, x, f, zlim = c(0, 2), col = grey.colors(256, 0, 1), asp = 1) > > image(w, w, log1p(abs(F)), col = grey.colors(256, 0, 1), asp = 1) > > ``` > > > [!info]- Mathematica > > ```mathematica > > L = 2; W = 0.4; th = 30 Degree; h = 1; (* full length, full width, angle, height *) > > f[x_, y_] := h UnitBox[(x Cos[th] + y Sin[th])/L] UnitBox[(-x Sin[th] + y Cos[th])/W]; > > F[wx_, wy_] := L W h Sinc[L (wx Cos[th] + wy Sin[th])/2] Sinc[W (-wx Sin[th] + wy Cos[th])/2]; > > GraphicsRow[{DensityPlot[f[x, y], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 201, MaxRecursion -> 0, Exclusions -> None, > > ColorFunction -> GrayLevel, ColorFunctionScaling -> False, PlotRange -> {0, 2}], > > DensityPlot[Log[1 + Abs[F[wx, wy]]], {wx, -40, 40}, {wy, -40, 40}, PlotPoints -> 201, MaxRecursion -> 0, ColorFunction -> GrayLevel]}] > > ``` > > [!example]- `image-dft` — an $n \times n$ image and its centred 2-D DFT > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > n = 256; c = (n - 1) / 2; j, i = np.meshgrid(np.arange(n), np.arange(n)); x, y = j - c, c - i > > th = np.deg2rad(45); u, v = x*np.cos(th) + y*np.sin(th), -x*np.sin(th) + y*np.cos(th) > > f = ((abs(u) <= 23) & (abs(v) <= 3)).astype(float) # a 46 x 6 px bar at 45 degrees > > # or any picture: from PIL import Image; f = np.asarray(Image.open("photo.png").convert("L"), float) / 255 > > F = np.fft.fftshift(np.fft.fft2(f)) # the 2-D DFT, (0, 0) moved to the centre > > back = np.fft.ifft2(np.fft.ifftshift(F)).real # the way home: all of F (not |F|) recovers f to round-off > > print(np.abs(back - f).max()) > > plt.subplot(1, 2, 1); plt.imshow(f, cmap="gray") > > plt.subplot(1, 2, 2); plt.imshow(np.log1p(abs(F)), cmap="gray"); plt.show() > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > n = 256; c = (n-1)/2; [j, i] = meshgrid(0:n-1, 0:n-1); x = j - c; y = c - i; > > th = 45*pi/180; u = x*cos(th) + y*sin(th); v = -x*sin(th) + y*cos(th); > > f = double(abs(u) <= 23 & abs(v) <= 3); % a 46 x 6 px bar at 45 degrees > > % or any picture: f = double(rgb2gray(imread('photo.png'))) / 255; > > F = fftshift(fft2(f)); % the 2-D DFT, (0, 0) moved to the centre > > back = real(ifft2(ifftshift(F))); max(abs(back(:) - f(:))) % the way home: all of F (not |F|) recovers f to round-off > > colormap gray; subplot(1,2,1); imagesc(f); axis image off > > subplot(1,2,2); imagesc(log1p(abs(F))); axis image off > > ``` > > > [!info]- R > > ```r > > n <- 256; cc <- (n - 1) / 2; j <- matrix(0:(n - 1), n, n, byrow = TRUE); i <- matrix(0:(n - 1), n, n) > > x <- j - cc; y <- cc - i; th <- 45 * pi / 180 > > u <- x * cos(th) + y * sin(th); v <- -x * sin(th) + y * cos(th) > > f <- (abs(u) <= 23 & abs(v) <= 3) * 1 # a 46 x 6 px bar at 45 degrees (or: f <- png::readPNG("photo.png")) > > fftshift <- function(A) { h <- nrow(A) / 2; A[c((h + 1):nrow(A), 1:h), c((h + 1):ncol(A), 1:h)] } > > F <- fftshift(fft(f)) # fft() on a matrix is the 2-D DFT > > back <- Re(fft(fftshift(F), inverse = TRUE)) / n^2; max(abs(back - f)) # the way home (fftshift undoes itself for even n) > > par(mfrow = c(1, 2)); image(t(f[n:1, ]), col = grey.colors(256, 0, 1), asp = 1, axes = FALSE) > > image(t(log1p(Mod(F))[n:1, ]), col = grey.colors(256, 0, 1), asp = 1, axes = FALSE) > > ``` > > > [!info]- Mathematica > > ```mathematica > > n = 256; L = 46; W = 6; th = 45 Degree; > > f = Table[With[{x = j - (n + 1)/2, y = (n + 1)/2 - i}, > > Boole[Abs[x Cos[th] + y Sin[th]] <= L/2 && Abs[-x Sin[th] + y Cos[th]] <= W/2]], {i, n}, {j, n}]; > > (* or any picture: f = ImageData[ColorConvert[Import["photo.png"], "Grayscale"]] *) > > F = RotateRight[Fourier[N[f], FourierParameters -> {1, -1}], {n/2, n/2}]; (* numpy's fft2, then fftshift *) > > back = Re @ InverseFourier[RotateLeft[F, {n/2, n/2}], FourierParameters -> {1, -1}]; Max[Abs[back - f]] (* the way home *) > > GraphicsRow[{ArrayPlot[f], ArrayPlot[Log[1 + Abs[F]]]}] > > ``` > > [!quote] Attribution > The see-it-and-hear-it approach carried through this page continues from [[Fourier Series]], and with it the debt to Jez Swanson's *[An Interactive Introduction to Fourier Transforms](https://www.jezzamon.com/fourier/)* (source on [GitHub](https://github.com/Jezzamonn/fourier), MIT License), where hearing the harmonics while watching them was the original idea. The transform applets follow the course notebook `rect_fourier_transform.nb`.