# 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. ## Cutting a piece out Squiddy was read above as instructions to the inverse transform. The applet below puts that reading to the test. Pick a window on either graph with its two handles. Everything outside the window is set to zero, and the other panel shows what the piece that is left transforms into. Call the kept piece $g$ and its transform $\hat g$, so $g \leftrightarrow \hat g$ is always a genuine transform pair. <div class="applet" data-applet="rect-sinc-window"></div> **Cutting the rect.** Whatever survives of the rect is a narrower rect. A piece on $[c, d]$ has width $w = d - c$ and centre $m = (c + d)/2$, and the same integral as before gives $ \hat g(\omega) = \int_c^d A\,e^{-i\omega x}\,dx = A w\,\operatorname{sinc}\!\left(\frac{\omega w}{2}\right) e^{-i\omega m}. $ The width and the position do separate jobs in that formula. The width $w$ sets the sinc, with its first zero at $2\pi/w$, so a thinner piece has a wider transform. The position $m$ only multiplies by the phase $e^{-i\omega m}$, which has modulus one. Sliding a piece along $x$ leaves $|\hat g|$ alone and trades the real and imaginary parts back and forth, with period $2\pi/|m|$ in $\omega$. This is the shift theorem, and it is why an off-centre piece has a complex transform even though the transform of the whole rect is real. **Cutting the sinc.** Keeping $\hat f$ only on $[\alpha, \beta]$ and inverting gives $ g(x) = \frac{1}{2\pi}\int_\alpha^\beta 2AL\,\frac{\sin(\omega L)}{\omega L}\,e^{i\omega x}\,d\omega . $ For a symmetric window $[-\Omega, \Omega]$ the odd part of the integrand cancels and $g$ is real. Writing $\sin(\omega L)\cos(\omega x) = \tfrac12\bigl[\sin(\omega(L+x)) + \sin(\omega(L-x))\bigr]$ turns what is left into the sine integral $\operatorname{Si}(z) = \int_0^z \frac{\sin t}{t}\,dt$, $ g(x) = \frac{A}{\pi}\Bigl[\operatorname{Si}\bigl(\Omega(L+x)\bigr) + \operatorname{Si}\bigl(\Omega(L-x)\bigr)\Bigr]. $ This is the rect rebuilt from every frequency up to $\Omega$ and none beyond. Parseval says how good the rebuild is. Since $\int |f|^2\,dx = \frac{1}{2\pi}\int |\hat f|^2\,d\omega$, the energy $g$ misses is exactly the energy in the part of $\hat f$ that was thrown away, $ \frac{\|g - f\|^2}{\|f\|^2} = 1 - \frac{1}{2\pi\,(2A^2L)}\int_{\text{kept}} |\hat f(\omega)|^2\,d\omega , $ and the applet reports the share that was kept. One more fact does a lot of work here. Because $\int g\,dx = \hat g(0)$, the area under $g$ is $2AL$ when the window contains $\omega = 0$ and zero when it does not. Things worth trying with it. - **Only the main lobe** (the default). Squiddy's head without the arms is a smooth bump with no edges, peaking at $\frac{2}{\pi}\operatorname{Si}(\pi) \approx 1.179A$. It carries 90.3% of the energy, so the error is $\sqrt{0.097} \approx 0.31$ of $\|f\|$. The sharp corners of the rect live in the arms. - **Grow the window.** Press *grow the window* and watch $[-\Omega, \Omega]$ open from nothing out to $\Omega = 30$. The rebuilt pulse gets steeper edges and the ripples crowd toward $x = \pm L$, but the overshoot next to each jump never goes away. It settles near 9% of the jump. That is the Gibbs phenomenon from [[Fourier Series]], now seen as a statement about truncating a transform. - **No main lobe (mirrored).** Keep $|\omega| \ge \pi/L$ only. With $\omega = 0$ gone, $\int g\,dx = 0$, and what survives is the edges of the rect, two spikes at $x = \pm L$ with ringing around them and nothing in the middle. That is a high-pass filter, and it is how edge detection in images works. Only 8.7% of the energy lives out there. - **One side lobe, not mirrored.** A window that is not symmetric about $0$ gives a complex $g$, because a real $g$ needs $\hat g(-\omega) = \overline{\hat g(\omega)}$ and a lopsided window breaks that. Tick *mirror about 0* and the partner lobe comes back and $g$ is real again. - **Cut the rect, left half.** One piece of width $L$ centred at $-L/2$. $|\hat g|$ is a sinc with half the height and its first zero at $2\pi/L$, twice as far out as for the whole rect. Drag the band sideways and only the phase changes. - **Two slivers (mirrored).** Two narrow pieces of width $w$ at $\pm m$ give $\hat g(\omega) = 2Aw\,\operatorname{sinc}(\omega w/2)\cos(\omega m)$. The narrow width makes a wide sinc envelope and the separation stripes it with $\cos(\omega m)$. This is the two-slit interference pattern, the same one the *two rects* seed showed above. - **Wider than the rect.** A window that reaches past $\pm L$ changes nothing, because $f$ was already zero there. ## 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-sinc-window`, keep one window of $f$ or $\hat f$, zero the rest, and transform the piece that is left > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > L, A = 1.0, 1.0 > > F = lambda w: 2 * A * L * np.sinc(w * L / np.pi) # f-hat = 2AL sin(wL)/(wL); np.sinc(z) = sin(pi z)/(pi z) > > x = np.linspace(-6, 6, 601) > > wa, wb = -np.pi / L, np.pi / L # cut the sinc: keep the main lobe, zero the rest ... > > w = np.linspace(wa, wb, 2001) > > g = np.trapezoid(F(w) * np.exp(1j * np.outer(x, w)), w, axis=1) / (2 * np.pi) # ... and transform back > > c, d = -L, 0.0 # cut the rect: keep the left half, a rect of width d - c at m > > W = np.linspace(-30, 30, 1201); m = (c + d) / 2 > > ghat = A * (d - c) * np.sinc(W * (d - c) / 2 / np.pi) * np.exp(-1j * W * m) > > plt.subplot(1, 2, 1); plt.plot(x, g.real, x, g.imag, "--") > > plt.subplot(1, 2, 2); plt.plot(W, abs(ghat), W, ghat.real, "--", W, ghat.imag, ":"); plt.show() > > print("energy kept:", np.trapezoid(F(w)**2, w) / (2 * np.pi) / (2 * A**2 * L), " peak of g:", g.real.max()) # 0.903, 1.179 > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > L = 1; A = 1; > > 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, 601)'; > > wa = -pi / L; wb = pi / L; % cut the sinc: keep the main lobe, zero the rest ... > > w = linspace(wa, wb, 2001); > > g = trapz(w, F(w) .* exp(1i * x * w), 2) / (2 * pi); % ... and transform back > > c = -L; d = 0; m = (c + d) / 2; % cut the rect: keep the left half, a rect of width d - c at m > > W = linspace(-30, 30, 1201); > > ghat = A * (d - c) * sinc_(W * (d - c) / 2) .* exp(-1i * W * m); > > subplot(1, 2, 1); plot(x, real(g), x, imag(g), '--') > > subplot(1, 2, 2); plot(W, abs(ghat), W, real(ghat), '--', W, imag(ghat), ':') > > fprintf('energy kept: %.4f, peak of g: %.4f\n', trapz(w, F(w).^2) / (2 * pi) / (2 * A^2 * L), max(real(g))) % 0.903, 1.179 > > ``` > > > [!info]- R > > ```r > > L <- 1; A <- 1 > > sinc <- function(z) ifelse(z == 0, 1, sin(z) / z); F <- function(w) 2 * A * L * sinc(w * L) > > trap <- function(y, t) sum((y[-1] + y[-length(y)]) / 2 * diff(t)) > > x <- seq(-6, 6, length.out = 601) > > w <- seq(-pi / L, pi / L, length.out = 2001) # cut the sinc: keep the main lobe, zero the rest ... > > g <- sapply(x, function(xx) trap(F(w) * exp(1i * w * xx), w)) / (2 * pi) # ... and transform back > > cc <- -L; d <- 0; m <- (cc + d) / 2 # cut the rect: keep the left half, a rect of width d - cc at m > > W <- seq(-30, 30, length.out = 1201); ghat <- A * (d - cc) * sinc(W * (d - cc) / 2) * exp(-1i * W * m) > > par(mfrow = c(1, 2)); matplot(x, cbind(Re(g), Im(g)), type = "l", lty = 1:2) > > matplot(W, cbind(Mod(ghat), Re(ghat), Im(ghat)), type = "l", lty = 1:3) > > cat("energy kept:", trap(F(w)^2, w) / (2 * pi) / (2 * A^2 * L), " peak of g:", max(Re(g)), "\n") # 0.903, 1.179 > > ``` > > > [!info]- Mathematica > > ```mathematica > > L = 1; A = 1; F[w_] := 2 A L Sinc[w L]; > > {wa, wb} = {-Pi/L, Pi/L}; (* cut the sinc: keep the main lobe, zero the rest ... *) > > g[x_?NumericQ] := NIntegrate[F[w] Exp[I w x], {w, wa, wb}]/(2 Pi); (* ... and transform back *) > > reG[x_] := A/(2 Pi) (SinIntegral[wb (L + x)] - SinIntegral[wa (L + x)] + SinIntegral[wb (L - x)] - SinIntegral[wa (L - x)]); (* Re g, exact *) > > {c, d} = {-L, 0}; (* cut the rect: keep the left half *) > > ghat[w_] := A (d - c) Sinc[w (d - c)/2] Exp[-I w (c + d)/2]; (* = Integrate[A Exp[-I w x], {x, c, d}] *) > > Print["energy kept: ", NIntegrate[F[w]^2, {w, wa, wb}]/(2 Pi)/(2 A^2 L), ", peak of g: ", N[reG[0]]]; (* 0.903, 1.179 *) > > GraphicsRow[{Plot[{Re[g[x]], Im[g[x]]}, {x, -6, 6}, PlotStyle -> {Automatic, Dashed}], > > Plot[Evaluate@{Abs[ghat[w]], Re[ghat[w]], Im[ghat[w]]}, {w, -30, 30}, PlotRange -> All, PlotStyle -> {Thick, Dashed, Dotted}]}] > > ``` > > [!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`.