The discrete Hartley transform and the sine and cosine family

The discrete Hartley transform and the sine and cosine family

• Define the integral Hartley transform with the cas kernel and connect it with the Fourier transform by S=e-io and \psi=\mathrm{Re}\,S-\mathrm{Im}\,S.
• Compute the DHT of short sequences, use its reciprocity (the same formula for forward and inverse) and the shift, reversal, difference and quadratic-content theorems.
• Perform convolution with the DHT (two real multiplications instead of four) and obtain the power spectrum directly from the DHT.
• Implement the fast Hartley algorithm: bit-reversal permutation, the even/odd decomposition formula, and compare multiplication counts with the FFT.
• Distinguish DFT, DHT, DST, DCT-1 and DCT-2, and explain why DCT-2 is used to compress images.

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

A strictly reciprocal real transform

Given a real waveform V(t) we define \psi(\omega)=(2\pi)^{-1/2}\int V(t)(\cos\omega t+\sin\omega t)\,dt. Clearly \psi is a sum of double-sided cosine and sine transforms, so the inverse has the same form: V(t)=(2\pi)^{-1/2}\int\psi(\omega)(\cos\omega t+\sin\omega t)\,d\omega. These relations were presented by Hartley in 1942 and never disappeared from the technical literature but have not been well known. Forward and inverse are identical, and if V(t) is real so is \psi(\omega) (Bracewell, p. 293). In this chapter Bracewell drops the factor (2\pi)^{-1/2} and moves to the familiar notation with 2\pi\times frequency.

From Hartley to Fourier: reflect and add

Write \psi(\omega)=e(\omega)+o(\omega) with e,o the even and odd parts: e(\omega)=(2\pi)^{-1/2}\int V\cos\omega t\,dt and o(\omega)=(2\pi)^{-1/2}\int V\sin\omega t\,dt. Then the Fourier transform S(\omega)=(2\pi)^{-1/2}\int V(\cos\omega t-i\sin\omega t)dt=e(\omega)-io(\omega) (3): from \psi we extract S by simple reflections and additions. Conversely \psi(\omega)=\mathrm{Re}\,[S(\omega)]-\mathrm{Im}\,[S(\omega)] (4): from S we get \psi as the sum of the real part and the sign-reversed imaginary part (pp. 293 to 294).

Example: V(t)=e^{-t} for t>0

Dropping the factor (2\pi)^{-1/2}: the Fourier transform is S(\omega)=\dfrac{1-i\omega}{1+\omega^2} and the Hartley transform \psi(\omega)=\dfrac{1+\omega}{1+\omega^2} (Bracewell, p. 294, Fig. 12.1); \psi is real and clearly unsymmetrical. Numerical integration of \int_0^\infty e^{-t}\mathrm{cas}\,\omega t\,dt gives \psi(0) = 1, \psi(1) = 1, \psi(2) = 0.6 and \psi(-1) = 0. At \omega=1: \mathrm{Re}\,S = 0.5 and \mathrm{Im}\,S = -0.5, so \psi=\mathrm{Re}-\mathrm{Im} = 1. From the even and odd parts of \psi we recover S, as in the book's construction.

Reciprocity checked numerically

With V(t)=e^{-\pi(t-0.5)^2} (an off-centre Gaussian, neither even nor odd) and kernel \mathrm{cas}\,2\pi st: the Fourier transform is e^{-\pi s^2}e^{-i\pi s}, so \psi(s)=e^{-\pi s^2}[\cos\pi s+\sin\pi s]. Applying the same formula a second time to \psi returns V: at t=0.3 we have V = 0.8819 and the reconstructed value deviates by 1e-16. No factor i or \pm pairs as in the Fourier case; this is what Hartley called more symmetrical.

The view in the complex plane

Plot F(s)=R+iI as a locus in the complex plane (Fig. 2.10). Rotate the axes by -45^\circ: R'(s)+iI'(s)=[R(s)+iI(s)]e^{i\pi/4}. The component in the R' direction is the Hartley transform (times 1/\sqrt2), while the S' component adds no information if f is real. The locus can be reconstructed from R and I together; the absence of redundancy due to the symmetries of R and I permits the real f to be recovered from the R' direction alone, or from any direction except the R and S axes (p. 307). Measured at \omega=1 of the example: R' = 0.7071, and \psi/R' = 1.4142. For complex data the S' direction handles the remaining information and the Hartley transform becomes complex (Millane, 1994).

Physical meaning: phase from amplitude

The radiation field of an antenna is the Fourier transform of the aperture distribution, and a lens converts an aperture distribution into another similarly expressible, so it might seem that Hartley is a mere mathematical construct. But an intimate relationship exists: whatever property one transform possesses, a corresponding property is possessed by the other. The first apparatus used an optical beamsplitter and Michelson interferometer technique (Bracewell et al., 1985); later microwave elements extended the wavelength range (Villasenor and Bracewell, 1988). In the Hartley plane phase is encoded by amplitude alone; of interest since at ultraviolet and x-ray wavelengths phase is important but difficult to measure (p. 307).

Ralph Hartley and the history

Hartley was a Rhodes scholar, a Fellow of the Institute of Radio Engineers, in charge of telephone line research at the Bell Telephone Laboratories (1918 to 1929); to generations of electrical engineers his name was known through the Hartley oscillator, a textbook example of a sine-wave source with a single triode and feedback from a tapped inductance. The special issue of Proc. IEEE (March 1994, Olejniczak and Heydt) has hundreds of papers. When the fast Hartley algorithm appeared in 1984 there was some incredulity about the claimed speed, reality and reversibility; but hundreds of papers appeared, chips were manufactured and the physical significance was demonstrated in the laboratory with microwaves and visible light. It soon became clear that the relationship with Fourier amounts to full equivalency (Bracewell, pp. 296 to 297, 308 to 309).

Self-check, part 1

• What is the Hartley transform of e^{-t}H(t) at \omega=-1?
• How do you get the Fourier transform from the Hartley transform?
• Why is the Hartley transform of a real function real and invertible?

Hints: (1+\omega)/(1+\omega^2)=0; even part minus i times odd part; because the cas kernel is real and the inverse has the same form.

The DHT definition

For a discrete variable \tau=0,\ldots,N-1 and real f(\tau), H(\nu)=N^{-1}\sum_{\tau=0}^{N-1}f(\tau)\,\mathrm{cas}(2\pi\nu\tau/N), where \mathrm{cas}\,\theta=\cos\theta+\sin\theta. Compare the DFT F(\nu)=N^{-1}\sum f(\tau)e^{-i2\pi\nu\tau/N}. The inverse: f(\tau)=\sum_{\nu=0}^{N-1}H(\nu)\,\mathrm{cas}(2\pi\nu\tau/N) (Bracewell, p. 295). The proof uses the orthogonality \sum_{\nu}\mathrm{cas}(2\pi\nu\tau/N)\,\mathrm{cas}(2\pi\nu\tau'/N)=N if \tau=\tau' and 0 otherwise. Measured with N=16: the deviation of the orthogonality sum is 5e-14, and transforming and inverting random data returns the data with deviation 6e-15.

Quasi-symmetry and the strictly symmetric version

The DHT is quasi-symmetrical: forward and inverse differ only by the factor N^{-1} familiar from the DFT, and it is real. To make it strictly symmetrical, N^{-1} can be changed to N^{-1/2} in both directions, as Hartley did with his integrals. But in computing it is often convenient to stop short of multiplying by the constant, letting it be absorbed later into a normalisation or graphics scale factor. With N^{-1/2} the matrix \mathbf C_{\nu\tau}=N^{-1/2}\mathrm{cas}(2\pi\nu\tau/N) is symmetric and \mathbf C^2=\mathbf I: measured N=8, \|\mathbf C^2-\mathbf I\| = 2e-15. So \mathbf C is an involution and all its eigenvalues are \pm1 (p. 296).

From the DHT to the DFT and back

Split H(\nu)=E(\nu)+O(\nu) with even part E(\nu)=\tfrac12[H(\nu)+H(N-\nu)] and odd part O(\nu)=\tfrac12[H(\nu)-H(N-\nu)]. Then F(\nu)=E(\nu)-iO(\nu); conversely H(\nu)=\mathrm{Re}\,F(\nu)-\mathrm{Im}\,F(\nu) (p. 296). The power spectrum and Hermitian relations can accordingly be read directly from H. Measured with random data N=16: F=E-iO deviates from the FFT by 8e-16, and H=\mathrm{Re}\,F-\mathrm{Im}\,F deviates by 1e-15. For a DHT from existing software: run the complex FFT and subtract the imaginary parts from the real parts, using only half of the values (p. 324).

Short pairs: N=2, N=4, N=8

N=2: \{a\ b\} has DHT \tfrac12\{a+b\ \ a-b\}; for \{1,2\} we get 1.5 and -0.5 (p. 309). N=4: \{1,2,3,4\} has DHT \{2.5,-1,-0.5,0\} (p. 310). N=8: \{1,2,3,4,5,6,7,8\} has DHT \{4.5,\,{-1.7071},\,{-1},\,{-0.7071},\,{-0.5},\,{-0.2929},\,0,\,0.7071\} (p. 314; the book prints rounded values): the measurements are H(1) = -1.7071, H(3) = -0.7071, H(5) = -0.2929, H(7) = 0.7071, and H(6) = 0. Check: the sum of the data equals N\times H(0) = 36, and the sum of the DHT values equals f(0) = 1.

A 16-point example: the truncated exponential

Take f(0)=0.5 (the mean of the two values at the discontinuity of V(t)) and f(\tau)=e^{-\tau/2} for 1\le\tau\le15: a representation of the continuous function e^{-t/2}H(t) by N=16 samples (Fig. 12.2). The DHT closely resembles samples of the continuous transform taken at intervals 1/16; the small discrepancies are due to truncation and aliasing, exactly as with the DFT. Measured: H(0) = 0.1275, H(1) = 0.1386 against the continuous value \psi(1/16)/16 = 0.138, with \psi(s)=(0.5+2\pi s)/(0.25+4\pi^2s^2). This is the Hartley version of the DFT-versus-FT problem of module 14.

The binomial example: a smooth hump

The binomial sequence 1,6,15,20,15,6,1 represents samples of a smooth pulse; to obtain the simplest result the peak is placed at \tau=0: f=\{20,15,6,1,0,\ldots,0,1,6,15\} (N=16, even about 0). The DHT therefore equals the real DFT: H(\nu)=\{4,\ 3.56,\ 2.49,\ 1.32,\ 0.5,\ 0.12,\ 0.01,\ 0,\ \ldots\} and symmetric H(N-\nu)=H(\nu): a smooth pulse peaking at \nu=0 (Fig. 12.3). The measurements: H(1) = 3.56, H(2) = 2.49, H(4) = 0.5. (The copy of the book used here prints a 0 at \nu=4, whereas the symmetry H(12)=0.5 requires 0.5; this may be a scanning or printing error.)

Why N reals can replace 2N reals

At first sight it seems strange: the DFT has N complex numbers, 2N real numbers in all. But the Hermitian property means redundancy by a factor of 2: N/2 real numbers suffice for the cosine transform and N/2 for the sine transform, together forming N DHT coefficients with no degeneracy due to symmetry. The cas function is a sine wave shifted 45^\circ, so it responds equally to cosine and sine components. With a kernel \sqrt2\sin(\theta+\alpha) with arbitrary \alpha the responses are unequal but no information is lost unless \alpha=0,\pi/2,\ldots; the inverse uses a combination of \sin\theta and \cos\theta with \alpha-dependent weights (Bracewell, p. 298; the weights 1/(\sqrt2\sin\alpha) and 1/(\sqrt2\cos\alpha) below are our result). With \alpha=0.6: the data are recovered with deviation 5e-15.

Self-check, part 2

• What is the DHT of \{1,2,3,4\}?
• How does the inverse DHT formula differ from the forward one?
• Why does the DFT of a real sequence carry half redundant information?

Hints: \{2.5,-1,-0.5,0\}; only in the factor N^{-1}; because F(N-\nu)=\overline{F(\nu)}.

Theorems that match exactly

There is a Hartley theorem for every Fourier theorem; some match exactly: \sum H(\nu)=f(0) and \sum f(\tau)=NH(0). Reversal: f(-\tau) has DHT H(-\nu), with arguments outside 0\ldots N-1 interpreted modulo N. Quadratic content (Parseval): \sum|f(\tau)|^2=N\sum|H(\nu)|^2, like the DFT but with no complex conjugates, only real numbers (Bracewell, pp. 300 to 301). Measured for \{1,\ldots,8\}: \sum f^2 = 204 and N\sum H^2 = 204; the reversal theorem deviates by 1e-15 for random data. Cross-check with the DFT: \sum|F|^2=\tfrac12\sum[H(\nu)^2+H(N-\nu)^2] coincides with \sum H^2.

The shift and first-difference theorems

Unlike Fourier, the shift reads: the DHT of f(\tau+a) equals H(\nu)\cos(2\pi\nu a/N)-H(-\nu)\sin(2\pi\nu a/N) (p. 300), needed for the fast Hartley algorithm. The first-difference theorem: f(\tau+1)-f(\tau) has DHT [\cos(2\pi\nu/N)-1]H(\nu)-\sin(2\pi\nu/N)H(N-\nu) (p. 301). Measured for random N=16, a=3: the shift theorem deviates by 2e-15, and the difference theorem by 1e-15. Mind the sign: swapping the minus for a plus gives an error of about 0.7, so the theorem must be written with the right sign.

The general convolution theorem

If f=f_1*f_2 is the cyclic convolution then H(\nu)=H_1(\nu)H_{2e}(\nu)+H_1(-\nu)H_{2o}(\nu), where H_2=H_{2e}+H_{2o} are the even and odd parts (p. 299; as in the DFT convolution theorem f_1*f_2\supset NF_1F_2 there is also a factor N when H contains N^{-1}). Measured with two random N=16 sequences: the relation holds with the factor N, deviation 1e-14. If f_2 is even then H_{2o}=0 and H=NH_1H_2: one real product per element. Convolution by the DFT needs one complex multiplication per element, i.e. 4 real multiplications; by the DHT only 2 in the general case.

The convolution procedure with the DHT

In image processing and digital filtering one function is often even (a symmetric filter). Then f_1*f_2=\mathrm{DHT}[(\mathrm{DHT}\,f_1)\times(\mathrm{DHT}\,f_2)]: take two DHTs, multiply the two real sequences term by term, and take one more DHT (with a factor N depending on the normalisation convention). Beforehand extend both sequences to double length with zeros to make room for the result; the zeros may lead, trail or bracket the data. If one sequence is much shorter, as in ordinary digital filtering, direct convolution is recommended. With \{1,2,3,4\} and the even filter \{0.25,0.5,0.25\} padded to 16: the result matches direct convolution with deviation 4e-15; the result \{0.25,\ 1,\ 2,\ 3,\ 3.25,\ 1\} has largest element 3.

The two-dimensional DHT

An image f(\tau_1,\tau_2) as an M\times N matrix has the DHT as a real matrix H(\nu_1,\nu_2)=M^{-1}N^{-1}\sum\sum f\,\mathrm{cas}[2\pi(\nu_1\tau_1/M+\nu_2\tau_2/N)] and the inverse has the same formula without the factor (p. 299). It can be broken into N one-dimensional transforms. Benefit: the full spectral analysis of an image can be presented in a single diagram (Fourier needs two: real and imaginary parts, or amplitude and phase); one real array is stored. The two-dimensional power spectrum is the same for both transforms. Measured with a 6\times5 image: H=\mathrm{Re}\,F-\mathrm{Im}\,F with F the 2D DFT, deviation 6e-16; and transforming twice returns the image with deviation 4e-15. Three dimensions: Hao and Bracewell (1987); four: Buneman (1987).

The cas-cas transform

In two dimensions one may use the kernel \mathrm{cas}\,2\pi\nu_1\tau_1\,\mathrm{cas}\,2\pi\nu_2\tau_2, which is separable (Perkins, 1987). The cas-cas transform is also reciprocal (Bracewell, 1983) and is related to the standard DHT the way choosing between Fourier components \cos2\pi(\nu_1\tau_1+\nu_2\tau_2), corrugations with parallel null loci, and \cos\cos with null loci at right angles. Separability has not proved a basis for preference, so the kernel \mathrm{cas}[2\pi(\nu_1\tau_1+\nu_2\tau_2)] became the standard (p. 300). Measured on an 8\times8 image: the cas-cas transform applied twice returns the image (deviation 6e-15) but differs from the standard DHT: largest difference 1.52.

The power spectrum directly from the DHT

If the power spectrum |F(\nu)|^2=R^2+X^2 is wanted one need not go through complex numbers: |F(\nu)|^2=\tfrac12[H(\nu)^2+H(N-\nu)^2] (p. 314). For \{8,7,6,5,4,3,2,1\} (Problem 4): H=\{4.5,\ 1.7071,\ 1,\ 0.7071,\ 0.5,\ 0.2929,\ 0,\ -0.7071\}, the even part E(1) = 0.5, the odd part O(1) = 1.2071, so F(1)=E-iO = 0.5 - 1.2071i and |F(1)|^2 = 1.7071, matching the direct formula (2e-15). Deviation over all \nu for random data: 3e-16.

Exercises of the chapter

Problem 2: the sequence \{1,2,\ldots,2^P\} always has exactly one DHT element equal to 0. From P=2 (index 3) and P=3 (index 6) the rule is \nu=3N/4; for N=16 the index 12 and for N=64 the index 48. Proof: H(3N/4)\propto\sum(\tau+1)[\cos(3\pi\tau/2)+\sin(3\pi\tau/2)]; grouping four consecutive terms gives a zero sum. Problem 8: \{1234\} repeated four times has a nonzero DHT only at \nu=0,4,8 with values 2.5, -1, -0.5. Problem 10: \{-2,2,-2,2,\ldots\} is a Nyquist wave, its DHT has only H(4) = -2, the power spectrum concentrated at \nu=N/2.

Self-check, part 3

• How does the DHT shift theorem differ from the DFT one?
• When does convolution by the DHT need only one real product per element?
• How do you obtain the power spectrum from the DHT?

Hints: an extra term H(-\nu)\sin; when one function is even; \tfrac12[H(\nu)^2+H(N-\nu)^2].

Motivation and merits

A sequence of N real numbers has a DHT that is a sequence of the same length and is also real; applying the same formula a second time returns the original sequence. No need to manage real and imaginary parts in two arrays, or interleaved in one array of double length, and no inverse that differs from the forward. If the real and imaginary parts of the DFT are wanted they come from the even and odd parts of the DHT; the power spectrum directly. One-way FFT programs for real data skip the conjugate coefficients but the output is complex, so a separate inverse program must be stored, and neither can accept its own output. Bracewell states the FHT is faster than the complex FFT by a factor depending on circumstances (p. 308).

The even/odd decomposition formula

For N=2: \{a\ b\} has DHT \tfrac12\{a+b\ \ a-b\}. Split f into the even-indexed sequence (DHT H_1) and the odd-indexed one (DHT H_2), each of length N/2. Using the similarity theorem (stretching: the DHT repeats) for the even part and the shift theorem with a=-1 for the odd part: H(\nu)=\tfrac12\left[H_1(\nu)+H_2(\nu)\cos\tfrac{2\pi\nu}N+H_2(N-\nu)\sin\tfrac{2\pi\nu}N\right] with H_1,H_2 understood periodic with period N/2; in the text the factors \tfrac12 are postponed to the end (pp. 309 to 310). Measured: for random N=16 the formula matches the direct DHT, deviation 4e-15. In the first stages after permutation \cos,\sin take only the values 0, 1, -1 so only additions and subtractions remain.

Bit-reversal permutation

Permutation rearranges the sequence to bisect progressively down to pairs: for 8 elements the order is f(0),f(4),f(2),f(6),f(1),f(5),f(3),f(7). Write the index i in binary (000 to 111), reverse the digits and reconvert to decimal to give the new index j. Half of the indices (0, 2, 5, 7) remain unchanged because their binary representations are symmetrical, so skipping them saves time; the permutation reduces to two swaps 1\leftrightarrow4 and 3\leftrightarrow6. Scatter diagrams of j versus i break down into replicated cells (the kite for N=4, the frog for N=8) (pp. 321 to 323). Measured: the number of fixed indices is 4 for N=8 and 32 for N=1024 (2^{\lceil P/2\rceil}); for N=8 there are 2 swapped pairs.

A hand example: \{1,2,\ldots,8\}

Table 12.3: step 1 separates \{1357\} and \{2468\}, step 2 splits into \{15\}, \{37\}, \{26\}, \{48\}; the permuted sequence f_0=\{1,5,3,7,2,6,4,8\}. Each pair \{a\ b\} gives \{a+b\ \ a-b\} (dropping the factor \tfrac12): f_1=\{6,-4,10,-4,8,-4,12,-4\}, with 6=1+5, -4=1-5, 10=3+7, -4=3-7... After two stages combined with the cos and sin of eighths of a turn, 8H=\{36,\ -13.657,\ -8,\ -5.657,\ -4,\ -2.343,\ 0,\ 5.657\}; dividing by 8 gives H=\{4.5,-1.7,-1,-0.7,-0.5,-0.3,0,0.7\} (pp. 311 to 314). Measured: f_1(0) = 6, f_1(1) = -4, 8H(1) = -13.6569. Check: the sum of the data equals 8 times the first element of the DHT.

Counting multiplications: FHT against FFT

Implement the decomposition formula recursively and count the multiplications by factors other than 0,\pm1. For N=1024 the FHT uses 14344 real multiplications; the halving FFT has 3586 genuine complex multiplications, i.e. 14344 real multiplications (four real per complex), equal in this raw form: the advantage of the FHT is that it computes N real numbers (half the storage and output) instead of N complex numbers, and real input does not need the redundant Hermitian half of the spectrum. Bracewell reports by experiment on a personal computer: the figure of merit Y ranges from 3.1 to 5.8 for a translated FFT program and about 1.3 when the FHT is used to obtain the Fourier transform (p. 315): running-time numbers on one machine, not reproducible here.

The stripe diagram: time by component

Bracewell analyses running time against N in the stripe diagram (Fig. 12.5). Three kinds of component: fixed overhead that dies out when normalised by NP; precalculation of sines and cosines, proportional to N, so its stripe thins like P^{-1}; and permutation, N reassignments repeated P times, proportional to NP, so its stripe does not die out: any effort to streamline permutation is worthwhile (like replacing J/2 by 0.5\times J). The nondecaying zone has thickness 0.038 s on the author's computer, so the limiting time is 0.038NP; near P\approx10 it is about one third higher. Stages 1 and 2 use the short equations of Table 12.2 (no trigonometry) for speed (pp. 315 to 317). The radix-4 algorithm saves about 25 percent of the running time, for N a power of 4 (16, 64, 256, ...) (p. 323).

The matrix representation and the DFT factorisation

One can write \mathbf H=\mathbf L_P\cdots\mathbf L_2\mathbf L_1\,\mathbf P\,\mathbf f (with \mathbf P the permutation matrix and \mathbf L_i sparse matrix operators of level i), and proceed to the DFT \mathbf F=\mathbf W\mathbf f by multiplying with a (complex) conversion matrix; combined this gives a new factorisation of the DFT matrix \mathbf W with W=e^{-i2\pi/N} (pp. 317 to 320). \mathbf L_1 and \mathbf L_2 have only 2N nonzero elements compared with 3N for the higher levels; \mathbf L_1 has sine elements on lines of slope 45^\circ. Measured: the 16\times16 cas matrix has 224 nonzero elements out of 256 (the zeros are exactly where \mathrm{cas}=0).

Self-check, part 4

• What is the bit-reversal order of 8 elements?
• How many indices stay fixed in the permutation of 8 elements?
• Why is permutation the component that does not die out in the stripe diagram?

Hints: 0,4,2,6,1,5,3,7; 4 (indices 0, 2, 5, 7); because it is proportional to NP.

The DFT-type cosine transform is not invertible

The discrete cosine transform \sum f(\tau)\cos(2\pi\nu\tau/N) (and likewise the sine) cannot be inverted to recover f(\tau) because components of the form \sin(2\pi\nu\tau/N) are integrated out and leave no trace. They are invertible when attention is restricted to even f (f(\tau)=f(N-\tau) for \tau=1\ldots N-1), and the integral cosine and sine transforms are invertible if f(x) is zero for negative x (Bracewell, pp. 301 to 302). Measured: the sequence \{1,2,3,4,5,6,7,8\} and its reversal give the same cosine transform (difference 2e-14) though they are different sequences (data difference 6).

The invertible DST and DCT-1 with half-integer frequencies

Change the frequencies to 0,\tfrac1{2N},\tfrac1N,\tfrac3{2N},\ldots: the fundamental frequency is 1/2N and the fundamental period 2N. The DST: F_s(\nu)=\tfrac2N\sum_{\tau=1}^{N-1}f(\tau)\sin(\pi\nu\tau/N), inverse f(\tau)=\sum_{\nu=1}^{N-1}F_s(\nu)\sin(\pi\nu\tau/N); applying the sum \sum f\sin twice returns \tfrac12Nf(\tau). The value f(0) leaves no trace since \sin is zero there. DCT-1 (Rao and Yip, 1990) is strictly symmetric with N+1 values: F(\nu)=\sqrt{2/N}\,k_\nu\sum_{\tau=0}^Nk_\tau f(\tau)\cos(\pi\nu\tau/N), k=1/\sqrt2 at both ends and 1 elsewhere (pp. 301 to 302). Measured N=16: the DST inverse deviates by 4e-15; applying it twice gives the factor 8; DCT-1 applied twice deviates by 2e-15.

Boundary-value problems: the string

A violin string of length L drawn aside at a noncentral point: the boundary conditions y(0)=y(L)=0 render cosine terms unnecessary so only \sin(\pi x/L),\sin(2\pi x/L),\ldots are needed to synthesise the string shape. If instead the condition is [dy/dx]_{x=0}=0 cosine components alone may suffice. It can be done with the full Fourier transform but is not customary since it yields both kinds (Bracewell, pp. 302 to 303). Our example (a string plucked at a=0.3, L=1): y(x)=x/a for x<a and (1-x)/(1-a) for x>a, sine coefficients b_n=\dfrac{2\sin n\pi a}{n^2\pi^2a(1-a)}: b_1 = 0.7807, b_2 = 0.2294; a 300-term sum gives y(0.5) = 0.7143 against the exact 0.7143, deviation 8e-08. Waves in a canal with closed ends analyse into cosines alone.

Differential equations and dimensions

Fields of physics involving the wave equation, Laplace's equation, the diffusion equation and other fundamental differential equations in two, three and four dimensions provide a wide range of application for numerical analysis using the discrete sine and cosine transforms (p. 303). The reason: natural boundary conditions choose the basis functions, and if the sine functions are solutions of the eigenproblem the transform turns differential operators into multiplication (modules 13 and 19: the eigensolutions of \phi''+\lambda\phi=0 are \sin k\pi t). A numerical illustration: the 50\times50 second-difference matrix has eigenvalues 2-2\cos(k\pi/51) and eigenvectors \sin(k\pi j/51): the deviation between the numerical eigenvalues and the formula is 2e-15.

Data compression and transform encoding

An image of a million samples each with 256 levels: a fourfold reduction in bits results if 2\times2 pixel blocks are replaced by one mean pixel, but the image blurs, describable as convolution followed by coarser sampling. Empirically many images have little content at high spatial frequencies, suggesting transmitting the Fourier transform of the image with the higher frequencies suppressed or deemphasised (for example quantised to fewer levels), and inverting at the receiver: called transform encoding. Further gain comes from one-dimensional encoding (run-length, TIFF, delta-row) and segmenting the image into parts of uniform statistics; text, ASCII and digitised speech benefit little (Bracewell, pp. 303 to 304). Measured: 2\times2 blocks reduce the bits by a factor 4; and a further reduction of 32 since newsprint offers only about 8 grey levels: 256/8 = 32.

DCT-2: the transform of image compression

Ahmed, Natarajan and Rao (1974) introduced {}^{\mathrm{II}}F(\nu)=\sum f(\tau)\cos[\pi\nu(2\tau+1)/2N] (normalised by \sqrt{2/N}\,a(\nu) with a(0)=1/\sqrt2), and it entered an international standard (JPEG). DCT-2 analyses into frequencies 0,\tfrac1{2N},\tfrac1N,\ldots like DCT-1 but the samples are shifted by half a unit, so the basis functions are no longer purely cosinusoidal; the basis vectors remain orthogonal. Bracewell stresses that this transform is not strictly invertible (applying it twice does not return the object): measured \|\mathbf C^2-\mathbf I\| = 1.99 but \mathbf C\mathbf C^T=\mathbf I deviates by 2e-15. The shift and convolution theorems become opaque. The claimed advantage: empirical performance on images whose statistics approximate spatially uniform first-order Markov statistics (p. 304).

Why DCT-2 suits first-order Markov data: energy compaction

A first-order Markov process has autocorrelation r(k)=\rho^{|k|} (module 14, AR(1)). With \rho=0.95 and block N=16, sort the coefficient variances (the diagonal of \mathbf TR\mathbf T^T) in decreasing order and keep the four largest: the fraction of energy retained is 0.9557 for DCT-2, 0.9364 for the normalised DHT, 0.9559 for the Karhunen-Loève transform (the optimum bound, module 13), and only 0.25 with no transform (4 samples). So DCT-2 nearly reaches the optimum bound; this is why it is empirically preferred for compression (a numerical test of the remark on p. 304; the exact values are our result).

A feel for indices: cosine, sine and dc in the DHT

Bracewell recalls that a frequency f is represented at two locations \nu=\pm Nf mod N; cosine components have equal values, sine components equal values with opposite signs (pp. 305 to 306). With N=8: \{1,0,-1,0,1,0,-1,0\} (\cos2\pi\tau/4) has DHT \tfrac12 at \nu=2 and 6: 0.5 and 0.5; \{0,1,0,-1,0,1,0,-1\} (\sin) has DHT 0.5 at \nu=2 and -0.5 at \nu=6; adding dc 2: \{3,2,1,2,3,2,1,2\} has DHT 2 at \nu=0. With the other transforms, DCT-1 and DCT-2 of these very sequences give different numbers because they are not Fourier: the frequency meaning of a DCT spectrum must be read with care (p. 306).

Self-check, part 5

• Why is the DFT-type cosine transform not invertible?
• How does the DST become invertible?
• Why does DCT-2 suit first-order Markov data?

Hints: the sine components integrate out; use frequencies 0,\tfrac1{2N},\tfrac1N,\ldots and drop \tau=0; it packs nearly all the energy into a few coefficients like the KLT.

Takeaways

• Hartley: \psi=\mathrm{Re}\,S-\mathrm{Im}\,S with kernel \mathrm{cas}=\cos+\sin; forward and inverse have the same form; it is the component of the Fourier locus along a 45^\circ axis.
• DHT: N real numbers in, N real numbers out; F=E-iO, H=\mathrm{Re}F-\mathrm{Im}F; \{1,2,3,4\}\to\{2.5,-1,-0.5,0\}.
• Theorems: shift H\cos-H(-\nu)\sin; Parseval \sum f^2=N\sum H^2; convolution H_1H_{2e}+H_1(-\nu)H_{2o}; power spectrum \tfrac12[H(\nu)^2+H(N-\nu)^2].
• FHT: bit-reversal permutation then the even/odd decomposition H=\tfrac12[H_1+H_2\cos+H_2(N-\nu)\sin]; the raw multiplication count equals the FFT but with real numbers only.
• DST and DCT-1 are invertible (half-integer frequencies); DCT-2 (JPEG) is orthogonal but not an involution and packs the energy of a first-order Markov process almost like Karhunen-Loève.