Relatives of the Fourier transform: multidimensional, Hankel, Mellin, z, Abel, Radon, Hilbert, fractional

Relatives of the Fourier transform: multidimensional, Hankel, Mellin, z, Abel, Radon, Hilbert, fractional

• Apply the two-dimensional Fourier theorems (rotation, shear, affine, moments) and the spherically symmetric three-dimensional transform.
• Use the zero-order and n-dimensional Hankel transforms under circular symmetry, and understand Fourier kernels in the broad sense through Hankel's theorem.
• Move between Laplace, Mellin and z, and use the serial product to multiply and divide z polynomials.
• Compute the Abel transform numerically and use the Abel-Fourier-Hankel ring, the projection-slice theorem and modified back-projection in tomography.
• Build the analytic signal with the Hilbert transform, find the envelope and instantaneous frequency, use causality (Hilbert pairs) and know the fractional Fourier transform.

Use the buttons, the dots or the ← → keys to move between slides.

Close and distant relatives

Many linear transforms in common use have a direct connection with either the Fourier or the Laplace transform. The closest is the generalisation to two or more dimensions, and the Hankel transforms of zero and higher order under symmetry. The Mellin transform is illuminated by the Laplace transform and is the tool by which the theory of Fourier kernels is built. The Hilbert transform simplifies when studied by Fourier, and there is an intimate relationship whereby the Abel, Fourier and Hankel transformations applied in succession regenerate the original function. The Radon transform is a generalisation of the Abel transform when circular symmetry is dropped, and the Abel-Fourier-Hankel cycle generalises to the projection-slice theorem, a thinking tool for reconstruction from projections (Bracewell, pp. 329 to 330).

The two-dimensional Fourier transform

F(u,v)=\iint f(x,y)e^{-i2\pi(ux+vy)}dx\,dy and the inverse f=\iint F\,e^{i2\pi(ux+vy)}du\,dv. Stretched membranes, antennas and arrays of antennas, lenses, diffraction gratings and pictures on television screens are two-dimensional systems. A separable function f(x)g(y) has transform F(u)G(v); for example the Gaussian e^{-\pi(x^2+y^2)} transforms into e^{-\pi(u^2+v^2)}. Measured by a sum on a grid of step 0.04: at (u,v)=(0.5,0.3) we get 0.3436, and the formula 0.3436 deviates by 6e-16. The two-dimensional convolution is f**g=\iint f(x',y')g(x-x',y-y')dx'dy' (one function rotated half a turn about the origin, displaced, multiplied and integrated) (pp. 330 to 331).

Product and convolution: two ways to express a function

The two-dimensional rectangle \Pi(x,y) (equal to 1 when |x|,|y|<\tfrac12) has two representations: \Pi(x,y)=\Pi(x)\Pi(y) (a product) and =[\Pi(x)\delta(y)]*[\Pi(y)\delta(x)] (a convolution of two perpendicular bars) (Fig. 13.2, p. 331). This helps recognise a function as a convolution and gives the transform \mathrm{sinc}\,u\,\mathrm{sinc}\,v. Discrete measurement: the two-dimensional convolution of a row a (5 elements) and a column b (4 elements) equals the outer product b\,a^T, deviation 1e-16, with total sum \sum a\sum b = 72.

Table of two-dimensional theorems

The one-dimensional theorems generalise readily (Table 13.1, pp. 332 to 333). The definite integral \iint f=F(0,0); the first moment gives the centre of gravity \bar x=-F_u(0,0)/(i2\pi F(0,0)); the second moment \iint x^2f=-F_{uu}(0,0)/4\pi^2 and \iint(x^2+y^2)f=-[F_{uu}+F_{vv}](0,0)/4\pi^2. For the Gaussian e^{-\pi r^2}: F(0,0) = 1 and \iint r^2f = 0.3183, both ways (direct and from the derivatives of F).

The rotation and shear theorems

Rotation: if f is rotated through an angle \theta on the (x,y)-plane then its transform is rotated through the same angle in the same sense on the (u,v)-plane: f(x\cos\theta-y\sin\theta,\ x\sin\theta+y\cos\theta)\supset F(u\cos\theta-v\sin\theta,\ u\sin\theta+v\cos\theta). Shear: f(x+by,y)\supset F(u,v-bu): if the function is sheared its transform is sheared to the same degree in the perpendicular direction (pp. 332 to 333). Measured with f=e^{-\pi(4x^2+y^2/4)} rotated by \theta=0.6 at (u,v)=(0.3,0.2): the two sides deviate by 4e-16; with a circular Gaussian sheared by b=0.5 at (0.4,0.3): deviation 1e-15.

The affine theorem

A function subjected to an affine transformation of the plane, f(ax+by+c,\ dx+ey+f), has the transform |ae-bd|^{-1}\exp\{i2\pi(ae-bd)^{-1}[(ec-bf)u+(af-cd)v]\}\,F\!\left[\tfrac{eu-dv}{ae-bd},\ \tfrac{-bu+av}{ae-bd}\right] (Bracewell et al., 1993; p. 333). It contains the rotation, shear, similarity and shift theorems as special cases. Measured with an anisotropic Gaussian (a,b,c,d,e,f=1.2,0.3,0.4,-0.2,0.9,-0.3) at (u,v)=(0.35,-0.2): the formula with a + sign in the exponent deviates by 1e-16, while flipping it to - deviates by 0.78: the sign matters.

The three-dimensional Fourier transform

The classical example: X-ray diffraction by crystals: F(u,v,w)=\iiint f\,e^{-i2\pi(ux+vy+wz)}dx\,dy\,dz. Vector notation \mathbf x,\mathbf s is convenient in n dimensions. With spherical symmetry f=k(r): K(s)=4\pi\int_0^\infty k(r)\,\mathrm{sinc}(2sr)\,r^2dr, also reciprocal. Table 13.4: the ball \Pi(r/2) (radius 1) has transform \dfrac{\sin2\pi s-2\pi s\cos2\pi s}{2\pi^2s^3}; its value at 0 is the volume 4.1888 (4\pi/3), at s=1 it is -0.3183 (-1/\pi), matching the numerical integral -0.3183. The three-dimensional Gaussian e^{-\pi r^2} is self-transforming, for example at s=0.5: 0.4559. The separable pairs f(x)g(y)h(z)\supset F(u)G(v)H(w) generate many more pairs (pp. 340 to 343).

Cylindrical and spherical symmetry

Circular symmetry about the z-axis (f independent of angle): H(s,w)=2\pi\iint h(r,z)J_0(2\pi sr)e^{-i2\pi wz}r\,dr\,dz (Hankel in r, Fourier in z). Cylindrical symmetry (f depending only on r): F(u,v,w)=K(s)\delta(w) with K(s)=2\pi\int k(r)J_0(2\pi sr)r\,dr: an infinite cylinder has as transform a sheet in frequency space. Table 13.4: the bar 2\Pi(x,y) has transform \mathrm{sinc}\,u\,\mathrm{sinc}\,v\,\delta(w); the slab \Pi(x) has \mathrm{sinc}\,u\,\delta(v)\delta(w); the cube \Pi(x,y,z) has \mathrm{sinc}\,u\,\mathrm{sinc}\,v\,\mathrm{sinc}\,w. Measured: the unit cube at (0.5,0.25,0.25) equals the product of three sincs 0.516 (two ways: product and a summed measure).

Self-check, part 1

• What does the two-dimensional rotation theorem say?
• Why does the sign in the exponent of the affine theorem matter?
• What is the transform of the unit-radius ball at 0?

Hints: the transform rotates through the same angle in the same sense; it sets the linear phase due to shift; the volume 4\pi/3.

Circular symmetry leads to Hankel

Optical systems are often built from circularly symmetrical components, and waves spreading in two dimensions from a source exhibit symmetry naturally. When f(x,y)=f(r), r^2=x^2+y^2, then F(u,v)=F(q), q^2=u^2+v^2. Change to polar coordinates and integrate over the angle (using J_0(z)=\tfrac1{2\pi}\int_0^{2\pi}e^{iz\cos\phi}d\phi): F(q)=2\pi\int_0^\infty f(r)J_0(2\pi qr)\,r\,dr and the inverse f(r)=2\pi\int_0^\infty F(q)J_0(2\pi qr)\,q\,dq (Bracewell, pp. 335 to 336). Call F(q) the Hankel transform of zero order of f: strictly reciprocal, as with the cos and sin kernels. The 2\pi in parentheses comes from measuring q in whole cycles per unit of r, the 2\pi before the integral from the area element 2\pi r\,dr.

Hankel pairs: disk, exponential, Gaussian

Table 13.2 (p. 338): the disk \Pi(r/2a)\supset a\,J_1(2\pi aq)/q (also called jinc); the exponential e^{-ar}\supset2\pi a\,(a^2+4\pi^2q^2)^{-3/2}; the Gaussian e^{-\pi r^2}\supset e^{-\pi q^2}; the ring impulse \delta(r-a)\supset2\pi aJ_0(2\pi aq); 1/r\supset1/q. Measured for a=1.5 at q=0.4: the disk by numerical integration 0.094 matches the formula 0.094 (jinc); the exponential with a=1.3: integral 0.3605 matches 0.3605; the Gaussian at q=0.5: 0.4559. Applying it twice to the exponential returns f(0.5) = 0.522 with deviation 7e-09: reciprocity.

Theorems for the Hankel transform

Table 13.3 (p. 339): similarity f(ar)\supset a^{-2}F(q/a); addition; circular two-dimensional convolution \supset F(q)G(q); Rayleigh \int|f|^2r\,dr=\int|F|^2q\,dq; definite integral 2\pi\int f\,r\,dr=F(0); second moment 2\pi\int r^2f\,r\,dr=-F''(0)/2\pi (up to the sign convention of the table); equivalent width =F(0)/f(0) or f(0)/F(0). A shift of origin destroys circular symmetry so there is no shift theorem. For the Gaussian e^{-\pi r^2}: 2\pi\int f\,r\,dr = 1, and Rayleigh gives 0.5 and 0.5.

Numerical evaluation of the Hankel transform

Bracewell: quite straightforward if the Bessel function J_0 is available. Alternative: take the Abel transform f_A(x) for x\ge0, make it an even function of double length by supplying values for negative x, and call an FFT; the Abel-Fourier-Hankel ring is explained later (p. 337). Measured for the Gaussian e^{-\pi r^2}: the direct J_0 integral at q=0.5 gives 0.4559, the Abel path followed by a Fourier sum (grid step 0.01) gives 0.4559, and the formula e^{-\pi q^2} gives 0.4559: three ways agree.

Fourier kernels in the broad sense

Two functions f,g related by g(x)=\int f(s)k(s,x)ds with inverse f(x)=\int g(s)k(s,x)ds (the same kernel) make k a Fourier kernel in the broad sense: 2\cos2\pi ax, 2\sin2\pi ax, \mathrm{cas}\,2\pi ax and J_0. Hankel's theorem: g(x)=\int_0^\infty f(s)(xs)^{1/2}J_\nu(xs)\,ds, or G(s)=\int g(x)(xs)^{1/2}J_\nu(xs)\,dx, and g(x)=\int G(s)(xs)^{1/2}J_\nu(xs)\,ds. Splitting off s^{1/2}, x^{1/2}: F(s)=\int_0^\infty x\,f(x)J_\nu(xs)\,dx and f(x)=\int s\,F(s)J_\nu(xs)\,ds (general order \nu; \nu=0 above) (pp. 339 to 340). With \nu=\pm\tfrac12: J_{1/2}(z)=(2/\pi z)^{1/2}\sin z and J_{-1/2}(z)=(2/\pi z)^{1/2}\cos z, so the cosine and sine are included in Hankel's theorem; measured at z=2.3: the two identities deviate by 5e-16.

The n-dimensional Hankel transform

With spherical symmetry in n dimensions the resulting one-dimensional transform is F(q)=2\pi q^{1-n/2}\int_0^\infty f(r)J_{n/2-1}(2\pi qr)\,r^{n/2}dr (p. 343). With n=1,2,3 we recover the Fourier and Hankel transforms and the case of spherical symmetry stated before (recalling J_{1/2} and J_{-1/2} are sine and cosine). Measured for the Gaussian e^{-\pi r^2} at q=0.5: n=1,2,3 all give the same 0.4559 (e^{-\pi q^2}), maximum deviation 3e-16. This is why the Gaussian is the fixed point of the Fourier transform in every number of dimensions.

Fig. 13.4: Hankel pairs pictured

• A ring impulse \delta(r-a) becomes 2\pi aJ_0(2\pi aq): a set of concentric wave rings (diffraction from a circular slit).
• A disk (circular aperture) \Pi(r/2a) becomes the jinc: the Airy rings, with zeros at the roots of J_1; the first zero at 2\pi aq\approx 3.8317 (the value 3.8317).
• The Gaussian transforms into itself with no side rings, an advantage for light beams; 1/r transforms into 1/q.

Self-check, part 2

• Why is the zero-order Hankel transform strictly reciprocal?
• What is the Hankel transform of a disk of radius a?
• Which kernel does Hankel's theorem give for \nu=\tfrac12?

Hints: because it derives from the symmetric 2D Fourier transform; aJ_1(2\pi aq)/q; the sine, since J_{1/2}\propto\sin z/\sqrt z.

The Mellin transform

For complex s: F_M(s)=\int_0^\infty f(x)x^{s-1}dx (Bracewell, p. 343). Put x=e^{-t}: dx=-e^{-t}dt, x^{s-1}=e^{-(s-1)t} so F_M(s)=\int_{-\infty}^\infty f(e^{-t})e^{-st}dt: the Laplace transform of f(e^{-t}); replotting a function of time as a function of e^{-t} compresses all positive time into the range from 1 to 0 (Fig. 13.5). The inverse is f(x)=\tfrac1{2\pi i}\int F_M(s)x^{-s}ds. Measured: for f=e^{-x} we have F_M(s)=\Gamma(s); at s=2.5 this is 1.3293 by both the direct integral and the Laplace integral \int\exp(-e^{-t})e^{-st}dt (deviation 3e-12).

A table of Mellin pairs

(Table 13.5, p. 345.) Measured: \sin x at s=0.5 gives 1.2533 (\Gamma(\tfrac12)\sin\tfrac\pi4) through an oscillatory integral; \tfrac1{1+x} at s=0.3 gives 3.8832 (\pi/\sin0.3\pi); and the box H(1-x) at s=2 gives 0.5.

Mellin as a moment

F_M(s) is the (s-1)th moment of f. Area F_M(1); first moment F_M(2); second moment F_M(3); abscissa of centroid \bar x=F_M(2)/F_M(1); radius of gyration [F_M(3)/F_M(1)]^{1/2}; second moment about the centroid F_M(3)/F_M(1)-[F_M(2)/F_M(1)]^2 (Bracewell, p. 344, normalised form). Example f=\Pi(x-\tfrac12) (equal to 1 on (0,1)): F_M(s)=1/s; centroid 0.5, radius of gyration 0.5774, variance 0.0833 (1/12). The similarity theorem f(ax)\supset a^{-s}F_M(s): with e^{-x}, a=2, s=2.5 we get 0.235.

Mellin and the Fourier-kernel condition

From g(x)=\int f(s)k(xs)ds and f(x)=\int g(s)h(xs)ds, taking the Mellin transform of both gives K_M(s)H_M(1-s)=1, i.e. the solving kernel for integral equations with kernel k(ax). When h=k this is the condition for k(ax) to be a Fourier kernel: K_M(s)K_M(1-s)=1 (Bracewell, pp. 346 to 347). Check with k(x)=2\cos2\pi x: K_M(s)=2(2\pi)^{-s}\Gamma(s)\cos\tfrac12\pi s; at s=0.3: K_M(0.3) = 3.0715 (numerical integral 3.0715), and the product K_M(0.3)K_M(0.7) = 1. This proves through Mellin that the cosine kernel is reciprocal.

The z transform: a polynomial of samples

A signal whose sample values at equispaced times constitute the full information (pulse modulation, sampled-data control systems). With interval 1 and f(t) known at t=0,1,\ldots,n, the polynomial F(z)=f(0)+f(1)z^{-1}+f(2)z^{-2}+\cdots+f(n)z^{-n} is the z transform (p. 347). The impulse string g(t)=\sum f(n)\delta(t-n) is fully equivalent to the sample set. The relation with Laplace: F_L(p)=\sum f(n)e^{-np}, so with z=e^{p} it is exactly the z transform; at p=0.7 the measurement gives 5.2559 equal to F(e^{0.7}) (deviation 9e-16). On the unit circle z=e^{i2\pi k/N} we get the DFT (module 14): F(z_k) deviates by 1e-14 from the FFT.

The serial product is a polynomial product

A system with impulse response \{h(0)\,h(1)\,\ldots\} receiving the signal \{f(0)\,f(1)\,\ldots\} gives the serial product, whose z transform is the product F(z)H(z). Example: the signal \{2\,1\} and the response \{8\,4\,2\,1\}: (2+z^{-1})(8+4z^{-1}+2z^{-2}+z^{-3})=16+16z^{-1}+8z^{-2}+4z^{-3}+z^{-4}, so the output is \{16\,16\,8\,4\,1\} (16 the first element, 1 the last) (p. 348). Conversely dividing the output z transform by the input gives the response: polynomial division returns \{8\,4\,2\,1\} (remainder 0). Factoring polynomials gives theoretical convenience, which is why z is used.

z tables and theorems

Table 13.7 for sequences with formulas (p. 350): \{1,1,1,\ldots\}\to\dfrac1{1-z^{-1}}; \{n\}\to\dfrac{z^{-1}}{(1-z^{-1})^2}; \{n^2\}\to\dfrac{z^{-1}(1+z^{-1})}{(1-z^{-1})^3}; \{a^n\}\to\dfrac1{1-az^{-1}}; \{\cos\alpha n\}\to\dfrac{1-z^{-1}\cos\alpha}{1-2z^{-1}\cos\alpha+z^{-2}}. Measured at z=2 (a 300-term series sum): \sum n z^{-n} = 2, \sum n^2z^{-n} = 6, \sum\cos(0.7n)z^{-n} = 1.2729 matching the formulas. Theorems (Table 13.8): f(n-1)\supset z^{-1}F; nf(n)\supset-zF'(z) (measured deviation 2e-09); f_1*f_2\supset F_1F_2. The z table and the Mellin theorems are nearly identical because Mellin and inverse z both derive from Laplace.

Self-check, part 3

• What is the Mellin transform of e^{-x} at s=2.5?
• What is the z transform of \{2\,1\}*\{8\,4\,2\,1\}?
• What is the Fourier-kernel condition in Mellin terms?

Hints: \Gamma(2.5)=1.3293; the polynomial 16+16z^{-1}+8z^{-2}+4z^{-3}+z^{-4}; K_M(s)K_M(1-s)=1.

The Abel transform

When a circularly symmetric distribution is projected in one dimension (the response of a television camera scanning across a narrow line, or of a microdensitometer whose slit scans over a circular density distribution on a photographic plate) one meets the Abel transform: f_A(x)=2\int_x^\infty\dfrac{f(r)\,r\,dr}{(r^2-x^2)^{1/2}} (Bracewell, p. 351). It is also tied to fractional-order derivatives (heat conduction in solids, signal transmission through cables). Putting \xi=x^2, \rho=r^2 gives the modified Abel transform F_A=K*F with K(\xi)=\xi^{-1/2} for \xi<0 (a convolution form), so Fourier transforming; the solution is f(r)=-\dfrac1\pi\int_r^\infty\dfrac{f_A'(x)\,dx}{(x^2-r^2)^{1/2}} (pp. 352 to 353). Checks: f_A(0)=2\int_0^\infty f\,dr and \int f_A\,dx=2\pi\int f\,r\,dr.

Abel pairs: disk, semi-ellipse, triangle

Table 13.9 (p. 354): the disk \Pi(r/2a)\to2(a^2-x^2)^{1/2}\Pi(x/2a) (a semi-ellipse); Gaussian e^{-r^2/2\sigma^2}\to(2\pi)^{1/2}\sigma\,e^{-x^2/2\sigma^2}; the cone a\Lambda(r/a) \to a function with \cosh^{-1}; the ring impulse \delta(r-a)\to2a(a^2-x^2)^{-1/2}\Pi(x/2a). Measured: the disk with a=1 at x=0.6 gives f_A = 1.6 (2\sqrt{1-0.36}); the Gaussian e^{-\pi r^2} at x=0.5 gives 0.4559 (e^{-\pi x^2}, self-transforming); at 0: f_A(0)=2\int f\,dr for f=1-r (on r<1) is 1, and \int f_A\,dx = 1.0472 (2\pi/6=\pi/3).

Numerical Abel transform by a table of coefficients

After the change of variable, sum products K(\rho)F(\xi-\rho) at discrete steps. The values of K turn out to be the same however fine the interval, save for a normalising factor, so a universal table (Table 13.10) can be set up: the coefficient for the interval [n,n+1] is c_n=\int_n^{n+1}\rho^{-1/2}d\rho=2(\sqrt{n+1}-\sqrt n): c_0 = 2, c_1 = 0.828, c_2 = 0.636, c_3 = 0.536. The book's example F(p)=(10-p)^{1/2} has F_A(\xi)=\tfrac\pi2(10-\xi): at \xi=5 the sum \sum c_nF(5+n+\tfrac12) gives 7.79 (book 7.78) against the exact 7.854 (pp. 353 to 355). The inverse problem is computed from the bottom up with the same table.

The program for f(r)=1-r

An alternative avoiding the singularity as r\to x: define the function f(r) interpolated between samples, then integrate f_A(x)=2\int_0^\infty f(\sqrt{x^2+t^2})\,dt (with r^2=x^2+t^2), avoiding the root in the denominator. For f=1-r on [0,1] the book prints the values f_A(x) at x=0,0.1,\ldots,0.9: 1,\ .9651,\ .8881,\ .7853,\ .6658,\ .5368,\ .4045,\ .2753,\ .1564,\ .0575 (small mean absolute error); recomputed here: f_A(0.1) = 0.9651, f_A(0.5) = 0.5368, f_A(0.9) = 0.0575, agreeing with the printed values within 0.001 (p. 356).

The Abel-Fourier-Hankel ring

Start with an even function f(r): take the Abel transform, then the one-dimensional Fourier transform, then the Hankel transform: we return to the original function (Bracewell, 1956). Example f(r)=\delta(r-a): Abel gives 2a/\sqrt{a^2-x^2}\,\Pi(x/2a), its Fourier transform is 2\pi aJ_0(2\pi as), and Table 13.2 gives the Hankel transform of J_0 as \delta(r-a) (p. 357). Measured with the Gaussian e^{-\pi r^2} at r=0.6: the exact value 0.3227, and carrying out all three steps numerically (Abel integral, Fourier sum, Hankel sum) reproduces it with deviation 2e-07. This means Hankel = Fourier \circ Abel (problem 16 of the chapter).

The Radon transform and the projection-slice theorem

With an unsymmetrical density f(x,y) scanned in direction \theta: g_\theta(R)=\iint f(x,y)\,\delta(R-x\cos\theta-y\sin\theta)\,dx\,dy is the integral along the line x\cos\theta+y\sin\theta=R (the slit at distance R from the origin). This technique revolutionised medical diagnosis (CAT, NMR) and radio astronomy. The projection-slice theorem: the slice of F(u,v) at angle \theta equals the one-dimensional Fourier transform of the projection g_\theta(R) (Bracewell, 1956; pp. 356 to 358). Measured with the elliptical Gaussian f=e^{-\pi(x^2/4+y^2/0.25)} (a=2, b=0.5) at \theta=0.6: the projection g_\theta(0.3) by numerical integration 0.5399 matches the formula 0.5399; and the one-dimensional Fourier transform of the projection at q=0.4 equals the slice F(q\cos\theta,q\sin\theta): deviation 1e-16.

Reconstruction by modified back-projection

The projections give F(u,v) on spokes \theta=\mathrm{const}; the density of points is inversely proportional to radius, so multiply by |q| to correct: the factor \Pi(q/2M)-\Lambda(q/M)=|q|/M for |q|<M. Multiplication in the q domain corresponds to a convolution of the data g_0(R)=g(R)*(2M\,\mathrm{sinc}\,2MR-M\,\mathrm{sinc}^2MR). The procedure: (a) modify each scan by convolution, (b) back-project (spread the modified scan uniformly over the plane perpendicular to the R-axis) and (c) accumulate over the (x,y) plane: modified back-projection (Bracewell and Riddle, 1967; p. 358). Measured: the Fourier transform of the kernel at q=2 equals |q|/M = 0.333 (with M=6), at q=8>M it is 0; reconstructing two Gaussians from 180 directions: at the centre 0.9994 (exact 1), maximum error over four points 0.013.

Example: two Gaussians and their projections

• The object: a central Gaussian (width 0.5, amplitude 1) and a smaller one (width 0.3, amplitude 0.6) at (0.7,-0.5); the projection of each Gaussian is a Gaussian in R with centre shifted to R_c=x_c\cos\theta+y_c\sin\theta.
• Without filtering back-projection gives a blurred image (weight 1/|q|); filtering with the convolution kernel gives a sharp image: the value at the centre 0.9994.
• The remaining error comes from the finite bandwidth M and linear interpolation, not from the number of directions: 90, 180 and 360 directions agree to three digits.

Self-check, part 4

• What is the Abel transform of a unit-radius disk at x=0.6?
• What does the Abel-Fourier-Hankel ring say?
• Why does back-projection need the ramp filter |q|?

Hints: 2\sqrt{1-0.36}=1.6; Abel, then Fourier, then Hankel returns the original function; because the density of points on spokes goes as 1/|q|.

The Hilbert transform

Definition F_{\mathrm{Hi}}(x)=\tfrac1\pi\,\mathrm{P.V.}\!\int\dfrac{f(x')\,dx'}{x'-x}=\left(-\dfrac1{\pi x}\right)*f(x), the divergence at x=x' taken as the Cauchy principal value. The Fourier transform of -1/(\pi x) is i\,\mathrm{sgn}\,s (+i for s>0, -i for s<0): Hilbert transformation is a curious kind of filtering, leaving the amplitudes of spectral components unchanged but shifting their phases by \pi/2 according to the sign of s. Applying it twice reverses the phase of every component so HHf=-f. All cosines become -sines and all sines become cosines; even functions become odd (Bracewell, pp. 359 to 361). Measured on a periodic signal: the Hilbert transform of \cos differs from -\sin by 3e-15, and applying it twice differs from -f by 1e-15. The Hilbert transform of \Pi(x) is \tfrac1\pi\ln|(x-\tfrac12)/(x+\tfrac12)|: at x=1 it is -0.3497, at x=0.2 (a point inside the singular range) it is -0.2697, both matching numerical principal-value integrals.

The analytic signal

For a real f(t), associate the complex function f(t)-iF_{\mathrm{Hi}}(t) whose real part is f: the analytic signal; F_{\mathrm{Hi}} is the quadrature function. Example: the quadrature function of \cos t is -\sin t and the analytic signal is e^{it}: it bears to f the same relation as e^{it} to \cos t. It contains no negative-frequency components: it is obtained by suppressing negative frequencies and doubling: 2H(f)F(f). Measured with a random narrow-band signal: computing f-iF_{\mathrm{Hi}} (filter i\,\mathrm{sgn}\,s) matches \mathrm{scipy.signal.hilbert} with deviation 5e-16, and its spectrum at negative frequencies is below 1e-16 (Bracewell, pp. 361 to 363).

Envelope and instantaneous frequency

If a mean angular frequency \omega can be assigned, the analytic signal is written V(t)e^{i\omega t}; V(t) is a time-varying phasor; |V(t)| is the instantaneous amplitude, or envelope, and the time rate of change of phase is the instantaneous frequency (Bracewell, pp. 362 to 363). Bracewell warns: \omega and V(t) are in general not uniquely defined, so when you see the term envelope be prepared for surprises. Measured: the AM signal (1+0.5\cos3t)\cos50t has envelope 1+0.5\cos3t with deviation 8e-14; the FM signal \cos(50t+2\sin3t) has maximum instantaneous frequency deviation \beta\Omega=6 from the carrier: 6.

Causality and Hilbert pairs

The impulse response I(t) of a physical system is zero for t<0: causal (\"physical realisability\" is a bad term since a system may be impossible to build for other reasons). Split I=E+O into even and odd with O=\mathrm{sgn}\,t\,E so I=(1+\mathrm{sgn}\,t)E, and the transform T(f)=G(f)+iB(f) with B=F_{\mathrm{Hi}}[G]: the real and imaginary parts of a causal transfer function are a Hilbert pair; B(f)=\tfrac{2f}\pi\int_0^\infty\dfrac{G(u)\,du}{u^2-f^2} (principal value) and G(f)=\ldots (Kramers-Kronig, pp. 363 to 364). Example I(t)=e^{-t}H(t): T=\dfrac1{1+i2\pi f}, G=\dfrac1{1+4\pi^2f^2}, B=\dfrac{-2\pi f}{1+4\pi^2f^2}. At f=0.3: B = -0.414 from the formula and -0.414 from a numerical principal-value integral, deviation 1e-16.

Computing the Hilbert transform

It seems simple: discretise, take the DFT, multiply by i\,\mathrm{sgn}\,s, invert. Simpler still: convolve with coefficients of -1/\pi x at integers. With a semispan of N=5 sample spacings, the 11 coefficients are \{.064\ .080\ .106\ .159\ .318\ 0\ -.318\ -.159\ -.106\ -.080\ -.064\} (0.064 at \pm5, 0.318 at \pm1). The concern is that the discarded coefficients have an infinite sum, hoping the two infinities cancel; and the sampling interval must be small against aliasing. The frequency response of the 11 coefficients has the right phase \pm\pi/2 but the amplitude is not 1: only a narrow band near s=0.08 has amplitude near 1 (maximum 1.006), at s=0.1 it is 0.972, at 0.25 it is 0.552, at 0.4 it is 0.18 (pp. 364 to 366). The Fourier way fails similarly due to aliasing.

The Hilbert transform with the Hartley transform

For data f(\tau), \tau=0,\ldots,N-1: take the DHT H(\nu), swap N/2-1 pairs of values, with one element of each pair changing sign: \{H(0),H(1),\ldots,H(7)\}\to\{H(0),-H(7),-H(6),-H(5),H(4),H(3),H(2),H(1)\} (for N=8); then take the inverse DHT to get the Hilbert transform (Pei and Jaw, 1989; Bracewell, pp. 366 to 367). Measured: for a zero-mean signal with no Nyquist component, the DHT swap agrees with the FFT i\,\mathrm{sgn}\,s filter to deviation 1e-15; with an overall sign flip the deviation is 2.21, so the signs must be followed exactly. (The dc and Nyquist components, treated differently by the two methods, are excluded from the test.)

The fractional Fourier transform: definition

Applying \mathcal F twice gives f(-x), four times returns f(x): if \mathcal F^a with integer a has consistency, additivity, commutativity and linearity, can \mathcal F^{1/2} exist that applied twice gives \mathcal F? In two dimensions \mathcal F^2 is a rotation by \pi and \mathcal F^4 by 2\pi. Physics: the field in an intermediate plane between the input plane and the transform plane (Fraunhofer diffraction), or in an optical fibre with quadratically graded index n=n_0(1-r^2/h^2) where the transform of order a appears at distance aL with L=\pi h/2\sqrt{2n_0} (Bracewell, pp. 367 to 369). The definition (Namias, 1980), with \phi=\tfrac12\pi a: \mathcal F^af(\xi)=\sqrt{1-i\cot\phi}\,e^{i\pi\xi^2\cot\phi}\int f(t)\,e^{i\pi t^2\cot\phi}e^{-i2\pi\xi t/\sin\phi}dt (normalised to the familiar Fourier units); at a=1 (\phi=\pi/2) it is the Fourier transform: measured with a Gaussian the deviation is 1e-12.

Additivity and Hermite-Gauss eigenfunctions

Bracewell: the orthogonal functions H_n(\sqrt{2\pi}t)e^{-\pi t^2} (H_n Hermite) are eigenfunctions of the fractional Fourier transform: \mathcal F^ah_n=e^{-in\phi}h_n, \phi=\tfrac12\pi a; and strikingly the quantum harmonic-oscillator wave functions are these same functions, also the modes of a graded-index fibre and of a laser cavity (p. 370). Measured for a=0.7 (discretised on 700 points): the largest deviation over n=0,1,2,3 is 1e-13; the additivity \mathcal F^{0.4}\mathcal F^{0.5}=\mathcal F^{0.9} deviates by 3e-14; the phase of n=1 at a=0.5 equals -\pi/4 = -0.7854. Applications: the Schrödinger equation (Namias), fibre optics, optical information processing, rotation of the Wigner distribution (time-frequency analysis, module 20) and the real fractional Hartley transform (Mendlovic, 1995).

Self-check, part 5

• What does applying the Hilbert transform twice give?
• What is the analytic signal of \cos t?
• Why are the real and imaginary parts of a causal transfer function a Hilbert pair?

Hints: -f; e^{it}; because I=(1+\mathrm{sgn}\,t)E so the spectrum is G plus the Hilbert transform of G.

Takeaways

• Two-dimensional Fourier: rotation and shear keep the structure; the affine theorem contains all cases; three dimensions with spherical symmetry K(s)=4\pi\int k\,\mathrm{sinc}(2sr)r^2dr.
• Hankel F(q)=2\pi\int f J_0(2\pi qr)r\,dr is reciprocal; the disk \to aJ_1(2\pi aq)/q; Hankel's theorem gives J_\nu, with \nu=\pm\tfrac12 the cos and sin.
• Mellin = Laplace after x=e^{-t}, a moment of order s-1; z is the polynomial of samples, the serial product = the polynomial product, z=e^{p} links Laplace.
• Abel \to Fourier \to Hankel returns the function; the projection-slice theorem; modified back-projection filters by |q| then back-projects.
• Hilbert shifts phase by \pm\tfrac\pi2, HH=-1, builds the analytic signal, the envelope, causal pairs; the fractional Fourier transform is additive with Hermite-Gauss eigenfunctions and eigenvalues e^{-in\phi}.