WO2000044104A1 - Methods of digital filtering and multi-dimensional data compression using the farey quadrature and arithmetic, fan, and modular wavelets - Google Patents

Abstract

Methods are presented for calculating the wavelet filters and their inverses which rely on a new method for sampling (UA, TA, A) either digital or analog data. These methods combine and extend to give novel procedures for non-reversible multi-dimensional data compression. For selected applications, this procedure improves achievable compression factors by an estimated one to three orders of magnitude and is well-suited to picture build-up or other iterative refinement. Combining these wavelet filters and their inverses with previous theoretical work furthermore provides novel methods for calculating Fourier and other transforms. In a preferred embodiment used to calculate the Fourier transform (1) and its inverse (34) applied to digital input data, the method replaces the fast Fourier transform and its inverse and provides improvements in achievable accuracy. The new sampling method is inherently multiscale, and the invention thereby obviates the usual Nyquist constraint on the meaningful bandwidth in terms of the number of samples. Finally, the invention provides novel and efficient analog-to-digital and digital-to-analog interface.

Description

METHODS OF DIGITAL FILTERING AND MULΗ-DIMENSIONAL DATA COMPRESSION USING THE FAREY QUADRATURE AND ARITHMETIC, FAN, AND MODULAR WAVELETS

Field of the Invention

It is well-known that digital filtering of digitally represented samples of analog data can be simplified and improved by performing the required processing in some transformed domain. Using both a transform and its inverse to simplify the digital or other manipulation of data, one first transforms the data from its initial form to an intermediate form, then applies a specified manipulation, and finally applies the inverse transform to recover the data in its initial form from its intermediate form. By way of example and not by way of limitation: the initial form of the data may be digital samples of analog data, and the intermediate form of the data may be the coefficients of the Fourier transform; the initial form of the data may be the coefficients of the Fourier transform, and the intermediate form of the data may be the coefficients in some wavelet or other expansion; the initial form of the data may be an analog signal, and the intermediate form of the data may be some digital representation of it; the initial form of the data may be digital samples of an analog signal, and the intermediate form of the data may be a quantization of the digital samples. For instance, transform coding methods for data compression in signal or data processing are realized in practice as transform of data/ quantization, storage or transmission, de-quantization/reconstruction. Another class of examples consists of finite-impulse response digital filters where the input data is filtered using the indirect calculation of convolution given by Fourier transform/multiplication/inverse Fourier transform. For analog data, sampling and filtering together determine fidelity and speed of manipulation.

The invention disclosed here provides a new method for sampling analog or digital input data which may be used to calculate new wavelet transforms as well as their inverse reconstruction algorithms. This immediately provides novel and efficient methods for data compression based on this new method of non-linear transform coding. In combination with these wavelet filters, the invention also provides new methods for calculating various classical transforms including the Fourier transform and its inverse. This immediately provides novel and efficient methods for digital filtering. By way of example and not of limitation, practical applications of the methods include non- speech audio compression, speech compression, speech recognition, speech synthesis, voice printing, audio filtering for hearing aids, still two-dimensional image compression, moving video compression, video compression for purposes of telephony, precision Fourier analysis, precision trigonometry, denoising, interpolation, medical imaging, geological imaging for recovery of oil or other resources, and other emission/detection apparatus. In specific applications, it may be necessary to store streams of input data in a buffer for separate processing or manipulation while concurrently collecting further input data whether analog or digital. It may also be useful in practice in particular applications to calculate and archive or otherwise store specific coefficients or constants and recover them from memory at run time rather than computing them at run time.

Description of Prior Art

In the landmark paper "Orthonormal bases of compactly supported wavelets, Communications in Pure and Applied Mathematics 41 , pp. 909-996 (1988), Daubechies presented a class of v. jlets which have been effectively applie* 1 numerous contexts; see, for instance, US Patents 5,453,945 by Tucker et al., 5,606,575 by Williams, as well as Daubechies' monograph "Ten Lectures on Wavelets" SIAM (1992) and its references. In another landmark paper, "An algorithm for the machine calculation of complex Fourier series" , Mathematics of Computation 19, pp. 297-301 (1965), J. W. Cooley and J. W. Tukey described a simplified algorithm for the calculation of Fourier series on the computer. The utility of this fast Fourier transform FFT is well known with thousands of important applications and implementations; see, for instance, US Patents 3,871,577 by Avellar et al., 3,881,100 by Works et al., 4,051,357 by Bonnerot, and "The DFT, an owner's manual for the discrete Fourier transform" , SIAM (1995) by W. Briggs and Van Emden Henson and its references. A requirement for both the FFT method and for the multiscale filtering of the Daubechies' wavelets is that the sampling of input data must be equally spaced though further methods of adaptive refinement of grids have also been developed in each context.

Other examples of transforms whose efficient digital implementations are important and to which the invention is relevant include the transforms of Hubert, Haar, Laplace, Bessel, Laguerre, Hermite, Chebyshev, Hotelling, Mersenne, and Fermat; see, for instance, US Patents 3,891,443 by Lynch et al. and 4,093,994 by Nussbaumer.

Invention Summary

Three new families of wavelets, respectively called arithmetic wavelets, fan wavelets, and modular wavelets, are defined on the standard circle of radius one in the complex plane. Any complex-valued function defined on this circle may be written essentially uniquely as a complex linear combination of each family, and the approximation of specified input data as a finite linear combination for each family leads to a corresponding wavelet filter. Each filter depends upon a new quadrature on the circle, called the Farey quadrature, which relies on a novel method of non-equally-spaced sampling. The output of each wavelet filter plays the role in the invention of an intermediate representation between input data and any one of a number of useful output transformations of data which may be computed from the intermediate quantity.

With either analog or digital input data, it is convenient to think of the given input data as spatial data. The inverse wavelet transform or reconstruction algorithm requires calculating spatial data from wavelet data and admits an especially convenient implementation: the values of the spatial function at certain non-equally-spaced grids of arbitrarily fine spacing may be computed exactly using a specified finite set of wavelet coefficients. This reconstruction algorithm combines with the wavelet filter into a binary cascade giving an efficient procedure for data compression which is especially well-suited to discontinuous input data and to iterative refinement of output data. Furthermore, the methods extend directly to procedures for compression of multi-dimensional data. Computer coding for the reconstruction algorithm itself is sufficiently abbreviated and low-level that it might be transmitted together with compressed data, for instance, for down-loading from satellite to home computer.

For the calculation of classical transforms of spatial input data, it is convenient to think of the desired output data, for instance, the Fourier coefficients, as frequency data. The transform from spatial to frequency data is accomplished in two main steps, the first of which approximates the input data by wavelets using the wavelet filter. The second step is the calculation ?ourier coefficients or other frequency dat lTe rrST- wavelets and requires formulas derived in previous theoretical work. Similarly, the calculation of the inverse transform depends upon the reconstruction algorithm and known formulas which express frequency data in terms of the wavelets. The calculation of a transform or its inverse admits an elegant implementation as a binary cascade.

Brief Description of Figures

FIG la illustrates the standard circle C of radius one in the complex plane, where C bounds the standard disk D of radius one. The complex numbers +1 and —1 are also indicated, and the chord of C with these endpoints is labeled by the matrix

J o

0 1

FIG lb illustrates a chord in the top semi-circle of C labeled by a matrix A = I , ) , together with two further "descendant" chords of C labeled by the matrix products

UA and TA as indicated, where U = ( ) and T = I _ ) . The three chords determine a triangle inscribed in C, and the vertices of this triangle are indicated as complex numbers.

FIG lc illustrates a chord in the bottom semi-circle of C labeled by a matrix A = , ) , together with two further "descendant" chords of C labeled by the matrix products U^~1A and T^~lA as indicated, where U^~l = ( __ι i ) ^and T^~l — _{n 1} 1. The three chords determine a triangle inscribed in C, and the vertices of this triangle are indicated as complex numbers.

FIG 2 depicts a classical figure in mathematics called "the Farey tesselation in Klein's model of hyperbolic geometry" . It is constructed beginning from the chord in Figure la by recursively taking descendants as in Figure lb and Figure lc, respectively, on the top and bottom semi-circle of C. The figure itself is classical, but the sampling of input data on C at the endpoints of the chords of the Farey tesselation in increasing generation is a novel aspect of the invention disclosed here.

FIG 3a illustrates the circle C, the chord labeled / and the two triangles in the Farey tesselation on either side of this chord. The vertices of these triangles are indicated inside the circle as complex numbers. Outside the circle are given four explicit combinations of cosine and sine functions, one such function next to each circular arc determined by the vertices, where θ is the usual angular coordinate on the circle. A function ϋι(θ) is uniquely determined by the figure, where on each such circular arc, ϋι{θ) takes the values of the nearby combination of cosine and sine functions. The other parts Figures 3b-3e likewise determine analogous functions ϋ_(θ) for other matrices A = I j . These functions taken together form the system of wavelets on which the invention is based, and they are unique to this invention. FIG 3b illustrates th< lord in the top semi-circle of C labeled th the matrix A, where a > c. As before, the figure determines a function on the circle denoted "0A (&) .

FIG 3c illustrates the chord in the top semi-circle of C labeled with the matrix A, where c >= a. As before, the figure determines a function on the circle denoted ~0 (Θ) -

FIG 3d illustrates the chord in the bottom semi-circle of C labeled with the matrix A, where α > — c. As before, the figure determines a function on the circle denoted X^~ _A {Θ) .

FIG 3e illustrates the chord in the bottom semi-circle of C labeled with the matrix A =, where — c >— a. As before, the figure determines a function on the circle denoted A (Θ).

FIG 4 illustrates the notation near the chord labeled by the matrix A = ( , J in the top semi-circle of C; this notation will be useful in several subsequent discussions.

FIG 5 provides a flow chart for the main driving routine in calculating the Fourier transform.

FIG 6 provides a flow chart for the main recursive routine in calculating the Fourier transform.

FIG 7 provides a flow chart for the procedure of updating coefficients in calculating the Fourier transform.

FIG 8 provides a flow chart for the procedure of processing terminal arrows of the recursion in calculating the Fourier transform.

FIG 9 provides a flow chart for the main driving routine in calculating the inverse Fourier transform.

FIG 10 provides a flow chart for the main recursive routine in calculating the inverse Fourier transform.

FIG 11 illustrates the region U consisting of all complex numbers with positive imaginary part. Map the disk D to U via the function z f→ i(z + l)/(z — 1); if two points in C are the endpoints of a chord in Klein's model of the Farey tesselation as in Figure 2, then draw the semi-circle in U passing through the corresponding two points which is perpendicular to the real axis. This produces another classical figure in mathematics, called the "Farey tesselation in the upper half-space model of hyperbolic geometry" , which is indicated in Figure 11.

Description of Preferred Embodiments

In order to disclose the methods, it is necessary to develop some notation and terminology. To this end, consider the standard disk D of radius one about the origin in the complex plane with its usual complex coordinate z — x + iy, where i = /^~~ϊ and x, y are real numbers, together with its boundary circle C as illustrated in Figure la. An arrow is by definition an oriented chord of the circle C labeled by a two-by-two matrix I , ) , where the initial point of the chord is given by the complex number £ ^ € C, and the final >int of the chord is given by ^ e C. A Dcial arrow called the doe (short for "distinguished oriented edge" ) is given by the diameter of C with endpoints — 1, +1 labeled by the matrix I = I „ . j as in Figure la. Also required subsequently are the matrices 5 = 1 I and Uⁿ — ( j , Tⁿ = ( „ j , for any integer n, which satisfy the identities S² — —I, SUS = — ^_1, STS — —U^~ of matrix products.

Inductively define two separate binary cascades, one of top arithmetic arrows and another of bottom arithmetic arrows, as follows:

Basis: The doe, which is illustrated in Figure la, is both a top arithmetic arrow and a bottom arithmetic arrow.

Top Induction: As in Figure lb, given the top arithmetic arrow labeled by A = I , 1 , there are also top arithmetic arrows labeled by

_ττ . ( a b \ _ -, _A f a + c b + d

UA = [ , , , , and TA = a + c b + d I \ c a

Bottom Induction: As in Figure lc, given the bottom arithmetic arrow llaabb_'eled by A = ( , J , there are also bottom arithmetic arrows labeled by

U-^'A = [ c — ^a a d , —^h b A J and T-^'A = [ ^{a ~} c ^{C b ~} a ^d

It is a classical theorem in mathematics that the chords underlying any two arithmetic arrows meet only at their endpoints, if at all, and that the family of all such chords decomposes the disk D into an infinite collection of triangles with vertices in the circle C as is illustrated in Figure 2; this figure is called "the Farey tesselation in Klein's model of hyperbolic geometry" . There is no suitable reference for this procedure in the literature, and it is for completeness that the explicit recursive definition of arithmetic arrows has been given here. An arithmetic arrow is said to be of generation g if it arises from g applications of the inductive clause starting from the doe in case of either top or bottom arrows; the doe has generation zero by convention.

Consider the infinite set Q C C consisting of the endpoints of chords underlying all arithmetic arrows, and say a point of Q is of generation g if it is the endpoint of the chord underlying an arrow of generation g and no less; —1 and +1 are of generation zero by convention. There is a canonical enumeration of the generation g points of Q clockwise in C starting from — 1 and the associated enumeration of Q itself by increasing generation. For example, the ordering begins with:

— 1 — 2z -2 - . 2 - i 1 - 2. -1 - 3- -1, +1, +ι, -ⁱ, __{1 + 2i} , 2 ' 2 +7' Ϊ+ 2Ϊ' -1 + 3* ' ^{' " •}

Given input data on the circle C, the sampling of it at the points Q C. C in this order of increasing generation will be referred to as the Farey quadrature. In practice, one might interpolate given ta as required at the sampling points of tl ^τarey quadrature or restrict the Farey quadrature to a circular arc in C. (In the subsequent discussion of mathematical basis, the sample points of the Farey quadrature will be seen to correspond to the rational points in the real line enumerated in increasing order, where the generation is the length of a corresponding continued fraction expansion.)

Notice that the doe separates the triangle with vertices — 1, +1, — i from the triangle with vertices —1, +1, +i. In the same way, an arbitrary arithmetic arrow separates two triangles complementary to all the arithmetic arrows in D. The endpoints of chords of these triangles in various cases are illustrated in Figures 3a-3e as functions of the matrix A = I , j labeling the arithmetic arrow. Figure 3a depicts the doe itself with A — I, Figure 3b and 3c depict top arrows with a > c and c > a, respectively, and Figure 3d and 3e depict bottom arrows with a > —c and — c > α, respectively. The notation inside the circles indicates the endpoints of the chords as complex numbers, and the notation outside the circles in each case of Figures 3a-3e is yet to be explained.

To this end, define a trigonometric function to be a real-valued function f(θ) = α cos (9 + β sin θ + 7, where θ is the usual angular coordinate on the unit circle C and , β, 7 are real numbers. Next, define a piecewise trigonometric function on the circle to be a function f{θ) so that there is a finite collection of circular arcs I_j C C, for j = 1 , . . . , J, any two of which meet at common endpoints in C, if at all, together with a collection of trigonometric functions t_j(θ) for j = 1, . . . , J, so that f(θ) = tj(θ) if θ € I₃. In other words, C is decomposed into a finite number of non-overlapping circular arcs, and on each such arc, / takes the values of a trigonometric function.

There is a particular family of piecewise trigonometric functions, one such arithmetic wavelet t5 (#) for each matrix A labeling an arithmetic arrow. These functions are defined in Figures 3a-3e, where the endpoints of the circular arcs are as described before, and each trigonometric function is written outside the circle next to its corresponding circular arc. There are thus five cases in the definition of arithmetic wavelet corresponding to the five cases in Figures 3a-3e already discussed, and in each case, the arithmetic wavelet itself is a piecewise trigonometric function defined using exactly four circular arcs. Arithmetic wavelets are a fundamental new ingredient of this invention and no doubt seem entirely ad hoc and unmotivated at this juncture of the discussion. The mathematical basis for them will be given later, but it is worth pointing out now that each arithmetic wavelet A {Θ) is a once-continuously differentiable function taking value zero at the points —1, +1, — i of the circle C, that is, for -9 = 0, — π/2, — π, except for ϋu-ι {θ) and ϋτ-ι (θ) which take the value zero only at the points ±1.

In order to explicate the definition of arithmetic wavelets, notice that any two- by-two matrix A = I , 1 determines a self-mapping MA '■ C — C of the circle C as follows:

_n -_{r /} ' ^{i _} (ta + b) - i(tc + d)

MA ' — t + ij {ta + b) + i(tc + d) ' and in fact, the image of the endpoints, —1 and +1, of the doe under M_A are the endpoints of the arithmetic arrow with label A. Construct from the function ϋj(θ) defined on C in Figure -. -he corresponding vector field on C given ϋj(θ) β_g , where ^■ g is the usual unit tangent vector field to C.

There is an explicit formula to be given presently for the transformation of suitable vector fields under the change of coordinates on C given by MA - Namely, the vector field t(θ)βg = [a cos θ + β sin #+7] transforms under MA to the vector field t^A(θ)βg, where i [2(cd - ab)β + (ci² - b²)(a + 7) + (a² - c²)( - 7)] cos θ t^A(θ) = j + [(ad + bc)β + bd(a + 7) - ac( - 7)] sin θ \ .

+ \ [2(cd + ab)β + (d² + b²)(a + 7) - (α² + c²)( - 7)]

More generally, given a piecewise trigonometric function f(θ) taking the values of the trigonometric function t_j(θ) on the circular arc I_j as before, the corresponding vector field f(θ)β_ρ transforms under MA to the vector field which on the circular arc MA (IJ) is given by t (θ)β_g .

Changing coordinates on C by M_A using this formula transforms the vector field 'ϋj(θ) β_g to another vector field s9_A{θ)βg, thereby defining the function $A {&) on C. Finally except for A = U^~x and A = T^~1, the corresponding arithmetic wavelet is given by

^■&A(Θ) = O_A(Θ) - [a cos θ + β sin θ + 7], where a, β, 7 are chosen so that ϋ(π) = ϋ(0) = ϋ(—π/2) = 0; the two special cases are ϋ_u-ι (θ) = ϋ_u-ι (θ) + 2 s θ and _ _τ-ι (θ) = ϋ_τ-χ (θ) - 2 sin θ

(which are defined in this manner to guarantee a symmetry of the upper and lower semi-circles of C). This is the abstract definition of arithmetic wavelets which is made explicit in Figure 3.

It is furthermore convenient to define

( ϋ_A(θ); ύ ^' A ≠ U- T- A{θ) = l ϋ_υ-ι (θ); ύ ^' A = U- ( ϋ_τ-ι (θ); if J = T-\ so that for every label A on an arithmetic arrow, ι9_>ι(π) = A (0) — ^A (—^/2) = 0.

The arithmetic wavelets enjoy an important renormalization property, namely,

^~>BA (M_A(Θ)) = M_A' (θ) ϋ_B(θ), where the prime denotes the usual derivative and where it is required that either both matrices B and A are matrix products of U⁺¹, T⁺¹ or both matrices B and A are matrix products of £/^-1, -T^~1.

Just as a function / = f(θ) may be expressed in its representation as a Fourier series / ~ ^~__n c_ne^{ιnθ '}, so too it may be expressed as a series / « ^~~ CA^A (Θ) of arithmetic wavelets, which together form a linearly independent set. The sum is over all arithmetic arrows, and the coefficients CA are called the arithmetic wavelet coefficients. The calculation of arithmetic wavelet c -icients from the input data / is called th rithmetic wavelet filter. The arithmetic wavelet coefficients are not quite uniquely determined by / as will be explained.

The reconstruction algorithm or inverse wavelet transform provides a method for calculating function values at points of the circle from arithmetic wavelet coefficients. To determine the value of the sum ^~_ CA A at a specified point q G Q C C, draw a line segment from q to the origin in -D; this segment is contained in a finite collection of triangles complementary to all arithmetic arrows, whose boundary edges taken together form a finite set of chords; the corresponding finite set of arithmetic arrows A enumerates the only coefficients β_A which affect the value of ___A _s ^e ^_A at q. This is an important finiteness property of arithmetic wavelets which is crucial to the efficiency of the reconstruction algorithm.

This completes the general discussion of arithmetic wavelets, the first family of new wavelets disclosed here, and the two further families of new wavelets are next defined.

For each label A on an arithmetic arrow, there are two fan wavelets

as well as two modular wavelets

ΨSA (Θ) = ]^~~ n ϋ_τ-_nA(θ). n>0

It follows immediately from the definitions that there are the following relations among the wavelets for any label A on an arithmetic arrow.

ΦA(Θ) = ΨU- IA(Θ) - ΨA(Θ) , ΦSA(Θ) = ΨU-ISA(Θ) - ΨSA(Θ), ϋ_A(θ) = ΦA (Θ) - ΦUA (Θ) = ΦSA (Θ) - ΦUSA (Θ)

= ΨU- I A (Θ) - 2ΨA (Θ) + ΨUA (Θ) = 1>U-ISA (Θ) - ^~ΨSA (Θ) + ΨUSA (Θ) .

On the other hand as is not at all obvious from the definitions, φ and ψi are given by the following explicit formulas as piecewise trigonometric functions

2 sin θ; if 0 < θ < π, φi = { -2 cos θ - 2 sin θ - 2; if π < θ < 3π/2, -2 cos θ + 2 sin θ + 2; if 3ττ/2 < θ < 2π, and

2 - 2 cos if O ≤ ø ≤ π, φι(θ)- 0; if π < θ < 2τr.

In fact, modular wavelets arise from ψ in exactly the same manner that arithmetic wavelets arise from ϋj, namely, Φι{θ)β_g transforms under MA to the vector field A {Θ)

— ,{θ) — [a cos θ + β sin θ + 7], where a, γ are chosen to guarantee that φ_A(0) = V_A (^) = _Λ (- TΓ/2) = 0: transforming ψι{θ) β_g using the change of coordinates given by s^ likewise gives rise to Φ_SA (Θ) ■ The fan wavelets also arise from φj in the analogous manner.

In contrast to arithmetic wavelets, neither the fan nor modular wavelets are linearly independent, and in each case, there is one linear relation for each arithmetic arrow, namely,

ΨU- ^IA ^{~ ~}ΦA + ΦUA = ΦU- ^SA - 2φsA + ΦUSA- On the other hand, the following collection of modular wavelets forms a basis

{ΦA - A labels a top arithmetic arrow and c > a} U {Φ_SA - A labels a top arithemtic arrow and a > c) U {ΦA '• A labels a bottom arithemtic arrow and — c > a) U {Φ_SA '■ A labels a bottom arithemtic arrow and a > — c}; there is an analogous basis for the fan wavelets.

Both fan and modular wavelets enjoy precisely the same finiteness property as arithmetic wavelets. Furthermore, both fan and modular wavelets enjoy renormalization properties, namely,

ΦBA (M_A (Θ)) = M_A' (θ) φ_B(θ),

ΦBA(M_A (Θ)) = M_A' (θ) φ_B(θ), where the requirements on A, B are as before.

This completes the necessary abstract definitions. The methods will first be presented using arithmetic wavelets for real-valued input functions with extensions to modular and fan wavelets and complex-valued or multi-dimensional input/output data to be described later. As a point of notation throughout these discussions, there will be the tacit assumption that the entries in the matrix denoted A are given by

The Fourier Transform

Included in the section entitled "Source Codes" is an implementation in the computer language C of a preferred embodiment of the method described in this section employing arithmetic wavelets.

Suppose that f(θ) is some specified real-valued input function defined on the circle. The calculation of the Fourier transform of / proceeds by two main steps:

Stepl : Choose a finite set S of arithmetic arrows and approximate / « _2A __S ^e _{A A} calculating the coefficients e^ in a manner to be described in terms of the values of / at specified points using the Farey quadrature.

Step 2: Substitute into the expression in Step 1 the known Fourier expansions ^~A {Θ) ~ __n ^c _n ^{e%nθ m} order to derive the approximation ^cn — / _, C-A C_n . The coefficients require n Step 2 are known theoretically and

for n 0, ±1 by

π i (n³ — n) c^Λ = -[(c - a)² + (b - d)²

(b + d) i(a + c)^~

-{(c + a)² + (b + d)²

[{b + d) +i(a + c)_

Furthermore, the coefficients c^A,c _r.A^A, „cA_₁ are given by

πc^ = Θ_h(c² -d²± 2icd) + θ_r[(b - d)² - (a - c)² + 2i(b - d)(a - c)]/2 + θ_t(a²-b²±2iab) + θ_t[(b + d)² - (α + c)² + 2i(b + d)(a + c)]/2,

πc^A = 2θ_h(c² + d²) - θ_r[(b - d)² + (α - c)²] + 20_f(α² + έ>²) - ^[(6 + d)² + (α + c)²] where the angles are

2cό _ _^__ 2(c. + -i)(α + c) - tan Vd Al _- Λι ) ' ^{~tan {}{b + d)²-{a + c)²

θ_t = tan^_1( ^2ab \ a --.- ^'(b-d)(a-c) ₂), ^r = tan (^_ .₂__<, -₂). b² α^ (b-d)²-(a-c)²'

It may be useful in practice in particular applications to calculate and archive or otherwise store the constants required in the previous expressions for c^A and recover them from memory at run time rather than computing them at run time.

The Fourier coefficients co.c±i of / thus require special treatment which does not present additional challenges. For simplicity, suppose that the input function f(θ) has these three coefficients fixed by the condition that f(θ) = 0 for θ — 0, — π/2, —-. More precisely, if f(θ) is any real- valued function on the circle, then let f(θ) denote its normalization defined by f(θ) = f(θ) + a cos θ + β sin θ + 7, where a, β, 7 are chosen so that / takes value zero at 0, —π/2, — π. The specified simplifying assumption on / thus amounts to the condition that f = f.

The procedure for calculating wavelet coefficients is slightly different depending upon whether the corresponding arithmetic arrow is a top or bottom arrow. However, the antipodal map θ 1— . - + θ transforms the bottom of the circle C to the top in such a way that it suffices to discuss only the case of top arrows. More precisely, the algorithm for bottom arrows arises from that for top arrows by:

• replacing the input data value f(θ) with the input data value f(θ + π),

• replacing the matrix A in the described manipulations for top arrows by the conjugate matrix SAS for bottom arrows, and

• replacing the trigonometric function cos θ + 3sin θ + 7 for top arrows by the trigonometric function —αcos θ — 3sin θ + 7 for bottom arrows. The algorithm and cor itcr coding for bottom arrows is theref. derived without difficulty from that for top arrows, and the subsequent discussion is specialized without loss of generality to top arrows.

In addition to the function f(θ) defined on the circle, suppose that there is also specified a bandwidth of interest N (i.e. the method should calculate the Fourier coefficients c_n only in the range \n \ < N) and a tolerance EPS (where contributions to the specified Fourier coefficients are regarded as negligible if they are of magnitude less than EPS). There are numerous variants of the method which depend upon other parameters which may be tailored to a given application, and some these variants will be discussed further.

Steps 1 and 2 are merged into a single binary cascade which keeps track of an ongoing approximation to the Fourier coefficients as follows. A cascade element is called an arrow- structure and is defined to consist of the specification of

• an arithmetic arrow A of generation g,

• the wavelet coefficient e^ on ϋ_A in the approximation of /,

• the coefficients , β, 7 of a trigonometric function defined as follows. The end- points of the chord corresponding to A decompose C into two circular arcs, exactly one of which contains none of the points —1, +1, —i. The sum ∑B≠A ^e _{B B} over all arithmetic arrows B of generation at most g except for A itself is given on this interval by the particular trigonometric function cos θ + β sin θ + 7.

The coefficients , β, η in an arrow-structure keep a lagged/updated running tally of the overall effect of what has come before it in the cascade; as a result of this technique, the method requires essentially no memory other than the storage of the running approximation to Fourier coefficients and the stack required for the recursive computation in the cascade of arrow-structures.

The next technical point involves a regularization in the calculation of arithmetic wavelet coefficients; this regularization is required because the approximation to / as a linear combination ∑..₄es &A&_A is not uniquely determined. Given an arrow-structure in the specified notation and referring to Figure 4, when sampling the input data at the point [_{fe d}|-^ ° _g] in the Farey enumeration with corresponding angular coordinate ΘQ, the one real input datum /(#o) must determine the two coefficients eu_A and ex A - According to the formulas presented in Figures 3a-3e, these new coefficients are constrained by

(6 + d)² - (a + c) 21 2(a + c)(b + d) f(θo) = (b + d)² + (a + c)² + β (c. + i) + (α + c) _ + 7

+ 4(euA - C A) + 8 β_A sgn (α — c)

(b + d)² + (α + c)² (6 + d)² + (a + c)² ' where sgn (α — c) is a sign defined to be +1 if α > c and —1 otherwise. There remains, however, one real degree of freedom in the specification of e_ _A, CTA , and it is eliminated by demanding the best least-squares fit to the values of / at the next generation of points _(2b+-.₎-._(2α+_c) ' _(d+_2d)- _(α+_2c) ^with re pective angular coordinates θι , θ₂, as illustrated in Figure 4. Updating the lagged trigonometric coefficients in accordance with the formulas presented in Figures 3a-3e produces trigonometric functions Vι(θ) and V_(θ) which give the values of the ongoing approximations at θ_\ and #₂, respectively, and whose coefficients are e icit inhomogeneous linear functions of e\ ex A , namely,

Vi θi) = vie + vm euA + ^v*t ^eTA. for Ϊ = 1, 2.

The coefficients are explicitly described most concisely in pseudo-code, where e = e_A, alp = α, bet = β, gam = 7, bpd = b + d and ape = a + c: if (a > c){ axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; ayl=2*(b*b-a*a); bxl=bet+(c*d-a*d-b*c-a*b)*4*e; byl-=4*a*b; cxl=gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; cyl=-2*(a*a+b*b); ax2=alp+4*(d*d-c*c)*e; ay2=2*(c*c-d*d); bx2=bet+8*e*c*d; by2=-4*c*d; cx2=gam-4*e*(c*c+d*d); cy2=2*(c*c+d*d);

} else{ axl=alp+4*e*(a*a-b*b); ayl=2*(b*b-a*a); bxl=bet-8*e*a*b; byl=4*a*b; cxl=gam+4*e*(a*a+b*b); cyl=-2*(a*a+b*b); ax2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; ay2=2*(c*c-d*d); bx2=bet+(a*d-a*b+b*c+c*d)*4*e; by2=-4*c*d; cx2=gam+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e; cy2=2*(c*c+d*d);

} vlc=cxl+ (((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl

+2*(b+bpd)*(a+apc)*bxl)/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); vlt=cyl+(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*ayl

+2*byl*(b+bpd)*(a+apc))/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); vlu=8./((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); v2c=cx2+(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ax2

+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); v2u=cy2+(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ay2

+2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); v2t=-8./((d+bpd)*(d+bpd)+(c+apc)*(c+apc));

It is easy to minimize [/(#ι) - Vi (#ι)] + [/(#.>) - ^2(^2)] as a function of euA, &TA subject to the given constraint, and this procedure provides a suitable regularization. Weighted least-squares .s or other statistical treatments using i it data values at further generations allow alternative regularizations.

Two other simple regularizations which have also been useful are: to require that euA = 0 if c > α and e _A = 0 if α > c, or, to require that CUA/CTA = —[a(a + c) + b(b + d)}/[c(a + c) + d(b + d)\; the latter regularization is related to imposing differentiability of the approximation at θo ■

The final technical point involves the stopping criteria and terminal processing for the cascade of arrow-structures. The basic stopping parameter is NVAN, where a branch of the cascade terminates whenever there have been NVAN consecutive generations of offspring whose contributions to all coefficients c_n in the bandwidth N have been of magnitude at most EPS. There are further parameters governing stopping criteria which may be introduced depending upon the scale and resolution of features in the input data, including a minimum or maximum number of generations, the size of the angle subtended by the chord of a terminal arrow-structure, and so on.

Another useful for several generation is the o priori estimate

de when NVAN consecutive generations of coefficients β_A have satisfied the condition that \e_A\ < ^ EPS |αc + k-|⁵/⁴.

When NVAN generations of offspring contribute negligibly to the specified bandwidth, the cascade is terminated with those offspring of greatest generation. Suppose that A labels the arrow of such a terminal offspring, adopt the notation in Figure 4, and let J . C denote the circular arc with endpoints fzjf , 3₃^ which contains none of the points —1, —i, +1 € C. The final update r. = alpi cos θ + beti sin θ + gam-, for i = 1, 2, of the running trigonometric function on J in accordance with Figures 3a-3e is prescribed in pseudo-code with notation as before as follows: if (a > c){ alpl=alp+e*2*(a*(a+c)+d*(d-b)-b*(b+d)-c*(c-a)); betl=bet+e*4*(c*(d-b)-a*(d+b)); gaml=gam+e*2*(a*(a+c)+c*(a-c)+b*(b+d)-d*(d-b)); al_P2=alp+e*4*(d-c)*(d+c); bet2-=bet+e*8*c*d; gam2=gam-e*4* (c*c+d*d) ;

} else{ alpl=alp+e*4*(a-b)*(a+b); betl=bet-e*8*a*b; gaml=gam+e*4*(a*a+b*b); alp2=alp+e*2*(a*(a-c)+b*(d-b)-c*(c+a)+d*(b+d)); bet2=bet+e*4*(a*(d-b)+c*(b+d)); gam2=gam+e*2* (a* (a-c)-c* (a+c)-b* (d-b)-d*(b+d)) ;

}

Now serially letting r denote the trigonometric functions τ_\ and r₂, again adopt the notation of Figure 4, where A labels the arrow corresponding to r. The updated r is compared with another trigonometric function σ defined by extrapolation which is uniquely determined fr ,-manding that it agree with / at the thr< .ngles comprising the next two generations, θ = ΘQ, Θ_\ , Θ₂ in the notation of Figure 4. The Fourier coefficients of the piecewise trigonometric function which takes values σ(θ) — τ(θ) on the circular arc J and value zero otherwise are easily computed using standard formulas and added to those of the ongoing approximation.

Other algebraic or statistical samplings of higher generation offspring allow alternative final processings.

For instance, one might save computational expense and simply terminate the cascade with no final processing at all.

Furthermore, the terminal arrow-structures could be stored for restart or iterative refinement capabilities.

One might alternatively record in an arrow-structure different transforms of the trigonometric function defined by α, β, 7. For instance, it is useful to take advantage of the renormalization property of wavelets and change these coefficients by transforming the vector field [ cos θ + β sin θ + η]β_g by the change of coordinates M_A ^l. This technique avoids the necessity of calculating certain large integers attendant to the calculation of wavelet coefficients as will be illustrated later.

One favorable aspect of the method is its advantageous mix of floating-point and integer operation types. Moreover, except for the coefficients c₀, c±i, the implementation of the algorithm is purely algebraic, that is, requires only addition and multiplication. Furthermore, by its very nature as a cascade, the algorithm is amenable to parallelization and efficient hardware implementation.

This completes the description of the implementation of the method of calculating Fourier coefficients using arithmetic wavelets.

Several points will next be addressed which are special to the implementation of the methods using modular wavelets. Again there is a basic Step 1, the calculation of the modular wavelet transform where the input function is approximated as f(θ) ∞ ___A 9AΦA (θ) _> and a basic Step 2, substitution of known expressions for the Fourier coefficients of the modular wavelets; these two steps are again conveniently merged into a binary cascade of arrow-structures in practice.

To elucidate Step 1 for modular wavelets, refer again to Figure 4 for the notation near the arithmetic arrow labeled by the matrix A. The single value of the input function f{θo) together with the ongoing calculation of the updated trigonometric function cos θ + β sin θ + 7 in the arrow-structure this time uniquely determines the modular wavelet coefficient

since ΦB (ΘQ) = 0 for every labeling matrix B of generation greater than A; there is thus no need for regularization with modular wavelets.

As to Step 2, recall that a basis of modular wavelets was given before by a collection ΦB (Θ) of wavelets where A is the label on an arithmetic arrow and either B = A or B = SA. In each case, 3 can calculate that φs (θ) has Fourier e: nsion

„.n0

ΦB (Θ) d% + d ^■Bi, + d + 9. Σ + n yZ + n ,-sι n° n

N >1 where

and _7r

Thus, the method of calculating Fourier transforms using modular wavelets is easily derived from the method using arithmetic wavelets. In fact, it is a simple matter to incorporate the two minor modifications just described and alter the included source code employing arithmetic wavelets to produce source code employing modular wavelets.

In the same way for fan wavelets, there is again a Step 1 and Step 2, which are merged into a binary cascade of arrow-structures. As for arithmetic wavelets, Step 1 for fan wavelets again requires a regularization scheme. As to Step 2, the formulas for Fourier coefficients of fan wavelets can be derived from those given before for modular wavelets using the identity φ_A(θ) = ΦA (Θ) — ΦUA (Θ) mentioned before. Again, the included source code is easily modified to employ fan wavelets.

Flow Charts for the Fourier Transform

Flow charts for various procedures used in the calculation of Fourier transforms are presented in Figures 5-8. For definiteness, the flow charts and source code refer to arithmetic wavelets without renormalization, but the renormalizing source code will also be included for completeness. For purposes of brevity in these flow charts, the arrow-structures defined before are referred to in these figures simply as "arrows" .

In Figure 5 is presented a flow chart for the main driving routine. The input data is normalized as indicated in program segment 1. There are separate recursions established in program segments 2 and 3 for the top and bottom of the circle respectively. Each recursion is estab- .ed with a call to the subroutine gener, w^' h then calls itself in turn. The Fourier coefficients are stored internally with an overall suppressed factor of π, and output data normalization of multiplying by - is accomplished in program segment 4. Output Fourier coefficients are finally displayed before exiting in program segment 5. Program segments 2 and 3 are entirely independent and could be performed in parallel; more generally, each of program segments 2 and 3 could be decomposed further into multiple parallel procedures.

A flow chart for the main recursive routine gener is given in Figure 6. The procedure starts with a test in program segment 6 to determine if:

• the argument envy >= NVAN, in which case there have been NVAN consecutive generations of negligible contributions, and

• the argument generation >= MING, in which case there have been a required minimum number of iterations in the recursion.

If both of these inequalities hold, then the procedure passes to program segment 15, where the cascade is terminated and the output data corrected with a call to the subroutine prune.

In the contrary case that the recursion should continue, the procedure passes to program segment 7, where the descendant arrows in the cascade are determined using the least- squares fit to the next generation of data as described before to compute the regularized wavelet coefficients of the descendants. These are combined with integer calculations to update the lagged trigonometric functions. The procedure then passes to program segment 8, where the ongoing approximations to Fourier coefficients are updated to include contributions from the two descendant arrows with a call of the subroutine charles for each descendant. The procedure continues with a test in program segment 9. If the contribution calculated in the subroutine charles for either descendant to any Fourier coefficient in the bandwidth was non-negligible, then the recursive argument envy is set to zero in program segment 10. In the contrary case that both contributions to all Fourier coefficients in the bandwidth were negligible, the control parameter envy is increased by one in program segment 13.

In any case, the procedure passes to program segment 11, where there is a test on the number of generations. In case too many iterations of the recursion have occurred, it is terminated in program segment 12, and suitable warning of non-convergence of the method is written to the output. In the contrary case that the maximum number of generations has not been reached, program segment 14 establishes the recursion by calling gener once for each of the descendant arrows with the updated control parameter envy and an incremented generation.

In Figure 7 is presented a flow chart for the subroutine charles. Calling the subroutine charles has the effect of updating the ongoing approximations to the Fourier coefficients for a single argument arrow and returning a flag which keeps track of whether any such contribution has been non-negligible. The flag is cleared in program segment 16, and several preliminary calculations are accomplished in program segment 17. The procedure passes to program segment 18, where it is determined whether to calculate the 0,+l,-l Fourier coefficients. If these are to be calculated, then the procedure passes to program segment 19, which calls a subroutine to update these three Fourier coefficients. Program segments 20 and 21 set the flag as required depending upon whether these three contributions are non-negligible. In any case, the proce re passes to program segment 22. If the yp over the bandwidth is complete, then all Fourier coefficients have been updated, so the subroutine charles returns the flag in program segment 23. In the contrary case, the current Fourier coefficient is updated in program segment 24. The flag is set as required depending on whether the current contribution was non-negligible in program segments 25 and 26. The iteration over bandwidth is established by the update in program segment 27 and the return to program segment 22.

In Figure 8 is presented a flow chart for the subroutine prune, which modifies the array of Fourier coefficients when terminating the cascade. The final update of lagged trigonometric functions to produce τι , T₂ is performed in program segment 28. There is one such function for each possible descendant of the argument arrow. The procedure passes to program segment 29, where further input data is sampled in order to compute two trigonometric extrapolations σι, σ , one such extrapolation for each possible descendant of the argument arrow. Each function σ, — Tj, for i — 1, 2, is truncated as described before, and the Fourier coefficients of the truncated functions are added to the ongoing approximations of Fourier coefficients in program segment 30. Program segments 31 and 32 implement the calculation of 0, +1,-1 Fourier coefficients if desired, and the subroutine prune is terminated with the return in program segment 33.

The Inverse Fourier Transform

Included in the section entitled "Source Codes" is an implementation in the computer language C of a preferred embodiment of the method described in this section employing arithmetic wavelets.

Since the method for calculating the inverse Fourier transform is similar in spirit and execution to that for the Fourier transform described in detail before, it need only be briefly discussed. The input data includes the specification of a collection of Fourier coefficients c_n in a given bandwidth N.

There are again two main steps in the calculation:

Stepl : Calculate arithmetic wavelet coefficients e from the given Fourier coefficients c„.

Step 2: Use the wavelet coefficients from Step 1 and the reconstruction algorithm to output values of the function ^~__A CA^A {Θ) at the Farey quadrature points on the circle.

For Step 1, there is again an explicit formula from previous theoretical work. Deep mathematical results prove that for n ψ 0, ±1,

e^inθ = cos nθ + i sin nθ

₌ _ [_{c + c e}« _{+ c}n _ιe-i_{fl] +} * ∑L_{ξn _{+ η}n_{) +} ^ _{ξn _ _{ηn) A {θ}s_{jt A} s. / S )

where the sum is over all arithmetic arrows, the underlying chord of which has endpoints ξ, η £ C, and

Substituting this expression for e^τnθ into the given Fourier expansion f(θ) ~ ∑_n c_ne^ιnθ gives

/(0) « c₀ + c[e^iθ + c'__{ι e}-^iθ + ∑ e_A ϋ_A(θ),

where e_A \ n ξⁿ + η-) + ^ ( - Vⁿ) \, n≠ Σ η + ξ

Cn 0,_l v - ξ and c'_e = c_t — ^~~~ c_n c", for I = 0, ±1. ≠O.±l

This expression f(θ) « c₀ + cie^iθ +

of /(6¹) constitutes Step 1 in the calculation of inverse Fourier transform.

In fact, the calculation of wavelet coefficients e_A in Step 1 can be implemented using mostly integer operations to improve run-time dramatically. To see this, notice that the formulas ξ = £z^, η = ^^ for the endpoints of the chord underlying the arrow labeled by the matrix A = I , J may be substituted into the expression for e_A to yield

_ ⁱ ^ " ^eA — /__ ^{Cn β}A > n≠0,±l where

_ [(bd - ac) + i(ad + bc)]ⁿ (a² + b²)ⁿ(c² + d²)^{n X} n{(bd + ac + i)ⁿ + (bd + ac - i)ⁿ) + i(bd + ac) {(bd + ac + i)ⁿ - (bd + ac - i)ⁿ)

In the obvious implementation of this formula to calculate e_A using integer operations, the intermediary integer expressions grow too large to be practical for n large if the arrow corresponding to A is of large generation.

This difficulty is overcome by wavelet renormalization, where these attendant large integer calculations are obviated entirely as is illustrated in the included source code.

Step 2 in the method for calculating inverse Fourier transforms depends upon the reconstruction algorithm to output function values on the circle taken by the linear combination f(θ) m c'₀ + c e^%θ + d__le^~'lθ + __A e_A ϋ_A^ ). The finiteness property of the reconstruction algorithm allows the exact calculation of these function values at the sample points Q . C in their ordering determined by the Farey quadrature. The two steps are again conveniently merged into a single binary cascade as follows. There is only one con 1 parameter, which is called SCALE, v. ;e the method is required to determine at least one output value in each circular arc subtending an angle SCALE. Thus, as SCALE is decreased the grid of output values on the circle is refined. A binary cascade is established using arrow-structures which is similar to that for the Fourier transform before. However, the wavelet coefficients are given in this case in closed form by the formula for e^ (so no least-squares fit is required). Because of the finiteness property of the reconstruction algorithm, exact output values are determined at each stage of the recursion. The cascade is terminated with an arrow whenever the chords underlying both descendant arrows subtend angles less than SCALE. Output data can be written to an array either when calculated to store data in the ordering of the Farey quadrature or when terminating the cascade to store data in the standard ordering on the circle. Furthermore, the terminal arrow-structures could be stored for restart or iterative refinement capabilities.

This completes the description of the implementation of the method of calculating the inverse Fourier transform using arithmetic wavelets.

Several points will next be addressed which allow the extension of the methods to modular and fan wavelets. Again there is a basic Step 1, the calculation of wavelet coefficients from Fourier coefficients, and a Step 2, the reconstruction algorithm; in each case of fan or modular wavelets, these two steps are again conveniently merged into a binary cascade of arrow-structures in practice. The finiteness condition on fan or modular wavelets mentioned before renders Step 2 entirely analogous to that for arithmetic wavelets. As to Step 1, simply substitute the identities ϋ_A = ΦA — ΦUA — ΦUA — 2 A +

described before into the expression above for e^ιnθ in terms of arithmetic wavelets to immediately derive corresponding expressions in terms of fan and modular wavelets.

Thus, the method of calculating inverse Fourier transforms using fan or modular wavelets is easily derived from the method using arithmetic wavelets. In fact, it is a simple matter to incorporate the minor modifications just described and alter the included source code employing arithmetic wavelets to produce source code employing fan or modular wavelets.

Flow Charts for the Inverse Fourier Transform

Flow charts for various procedures used in the calculation of inverse Fourier transforms are presented in Figures 9-10. For definiteness, the flow charts and source code refer to arithmetic wavelets without renormalization, but the re-normalizing source code will also be included for completeness. For purposes of brevity in these flow charts, the arrow-structures defined before are referred to in these figures simply as "arrows" .

In Figure 9 is given a flow chart for the main driving routine. Modified Fourier coefficients c₀' , c'_l , c'_₁ are computed in program segment 34. The formula given in Step 1 for the inverse Fourier transform is applied in program segment 35 to calculate the arithmetic wavelet coefficients in terms of Fourier coefficients for a particular family of nine arrows. These nine arithmetic wavelet coefficients are required to initialize the trigonometric functions for recursions established in program segment 36 with calls to the recursive subroutine gener. Recursions with initial conditions corresponding to A = U, T, T^~2, U^{~ l} , T-^l , T-^lU-^χ, U^~2 are undertaken in program segment 36. These recursions, as well as fi ler possible subrecursions derived from m, are independent and may be undertaken in parallel. This elaborate driving routine is convenient because e^ιnθ is expressed in Step 1 of the inverse Fourier transform method in terms of the normalized arithmetic wavelets _A which differ from the arithmetic wavelets X9A themselves only for A = U^~1, T^{~ 1}.

The procedure passes to program segment 37 to perform an overall correction of the function values using the modified Fourier coefficients c'_Q, c'₁, c_₁. Output data is displayed and the procedure terminated in program segment 38.

In Figure 10 is given a flow chart for the main recursive routine gener. In the cascade, the two descendant arrows of the argument arrow are generated in program segment 39 employing Step 1 of the inverse Fourier transform method to calculate the wavelet coefficients of the descendants. These expressions are then used to update the lagged trigonometric functions in program segment 39. A test is performed in program segment 40 to determine if the chord underlying each descendant arrow subtends a sufficiently small angle, as estimated by the comparison of an integer quantity with SCAT; the relation with the control parameter SCALE discussed before is given by SCAT = [exp sinh^_1(l/SCALE)]². In case one of the chords subtends too large an angle, the procedure passes to program segment 41. The recursion is established in program segment 41 with two calls to gener, where the arguments are given by the two current descendant arrows. In the contrary case, the procedure performs the final update of the lagged trigonometric functions in program segment 42, then stuffs the output array with the new function values in program segment 43, and finally returns in program segment 44.

Data Compression

As with any method for data compression based on transform coding, five basic manipulation are required, namely:

• wavelet filtering,

• quantization,

• storage and retrieval or transmission

• de-quantization, and

• reconstruction.

The wavelet filter has already been fully disclosed as Step 1 for the calculation of the Fourier transform, and the wavelet inverse filter or reconstruction algorithm has likewise already been described as Step 2 for the calculation of the inverse Fourier transform. Application-specific quantization is done according to psychovisual or psychoacoustic thresholds. Quantization and storage are furthermore merged with wavelet filtering into a single binary cascade as before, and retrieval or transmission are likewise merged with reconstruction into another single binary cascade. The reconstruction algorithm is exact. Furthermore, source code for it may be transmitted along with compressed data since the coding is sufficiently abbreviated and low-level.

As is well-known, useful compression by transform coding requires an efficient specifi- cation of the basis eler. .ts. This is especially advantageous for .' new wavelets: a top arrow of generation g labeled by the two-by-two integral matrix A is specified by g bits, namely, by writing A uniquely as the matrix product of factors U or T. For example, the matrix A = ( _ _ } factors as

so the coefficient for the wavelet ϋ_A is coded by the binary sequence TUUUTTU.

The method is also especially well-suited to progressive picture build-up or other iterative refinement as follows. The Farey quadrature determines a linear ordering on the set Q of quadrature points in the circle C as before. The ordered subset Q\{— i, ±1} C Q C C is put into bijective correspondence with the collection of all arithmetic arrows (where \ denotes the set-theoretic relative complement); as illustrated in Figures la-lc, the arithmetic arrow labeled by A = I , 1 corresponds to the complex number

[ttT-Zt_c) 6 C, for a top arrow; jϋlg-ϋ_αl_c) ^{6 g}> for a bottom arrow; i _ C, for the doe.

Thus, the ordering from the Farey quadrature on Q\{— i, ±1} determines a corresponding ordering on the set of all arithmetic wavelets. Storage, retrieval, transmission, and reconstruction of wavelet coefficients is accomplished at a specified input or output data resolution in this canonical ordering since energy compacts to the lesser wavelet coefficients. In the case of real-time compression and transmission or other manipulation not requiring retrieval, the method can be further improved by the manipulation of previously computed wavelet coefficients in parallel with the calculation of subsequent ones.

Complex and Multi-Dimensional Data

There is nothing remarkable about the extension of the invention to complex-valued input data f(θ) = u(θ) + i v(θ), where u and v are real- valued functions on the circle C. Concentrating for definiteness on arithmetic wavelets, one simple approach is to separately and independently transform u and v,

u(θ) « ^~~~ e_A ϋ_A(θ) and v(θ) ∞ ∑ f_A $_A(Θ),

using the techniques already described in order to calculate

f(θ) ^ ^~~ (e_A + if_A) ϋ_A(θ).

Another direct approach is to allow complex values for the input function / and solve for complex wavelet coefficients using the same algebraic formulas as before. -Turning now to multi-. lensional input/output data, suppose th there are -R > 1 input functions each of which has M > 1 independent variables, and define the torus C^M to be the iW-fold Cartesian product

C^M = C x • • • x C .

M times

Coordinates on C^M are given by -tuples θ = (#]. , #2_> • • • . ΘM) of angles, and the input data is specified by an -ft-tuple of multivariable real- or complex- valued functions

F(θ)

f²(θ). . . . , f^R(θ)).

Using the bijection described in the previous section on data compression, each quadrature point Z G Q\{—i, ±l} uniquely determines a corresponding arithmetic wavelet ϋz(θ); by convention, one also defines ϋ-ι(θ) = ϋ±ι(θ) — 1. Concentrating on arithmetic wavelets for definiteness, define the corresponding multiwavelet _ (θ) = ϋ_Zl ψ ϋz ) ^{■ ■ ■} Z_M (ΘM) .

As before, the input function F must be normalized. To this end, define a trigonometric monomial to be a product of the form X^l

^{■ • •} X^M (ΘM) , where each X^l(θ) is given by one of X*(θ) = cos θ, X^{(θ) = sin θ, or Xⁱ(θ) = 1. The values of F at the points in {— i, -£1}^M uniquely determine coefficients in a linear combination (θ) of trigonometric monomials so that (θ) = F(θ) for each point of {— i, ±l}^M. Just as in the one-dimensional wavelet filter already discussed, now one must consider the normalization F(θ) = F(θ) — v(θ), so F is zero at each point of {— i, ±1}^M .

Adopting the Farey quadrature in each factor circle C of the torus C^M, there is an induced lexicographic ordering on the set Q^M\{—i, ±l}^M C C^M of potential sample points Z € C^M . One separately approximates each coordinate function f (θ) « ∑ e_z ϋ^~z(θ), for i = 1, 2, . . . , R,

of F(θ) as a finite linear combination of multiwavelets with application-specific zonal sampling for optimum energy compaction; this calculation may employ least-squares or other regularization in each factor as well as renormalization as for the one-dimensional wavelet filter. The methods are again elegantly implemented in practice as M nested binary cascades. The final processing for the cascades of arrow-structures in the calculation of multi-dimensional Fourier transforms may involve adding suitable linear combinations of trigonometric monomials based on further sampling, or one may simply truncate with no final processing at all, as in the one-dimensional case.

These same remarks apply verbatim to fan and modular wavelets.

Mathematical Basis

"The decorated Teichmϋller space of punctured surfaces" , Communications in Mathematical Physics 13, pp. 299-339 (1987), written by the author of this patent application, began the study of certain abstract geometric coordinates called lambda lengths in the context of hype ->lic geometry on Riemann surfaces. Th_* leas developed in this and other related papers have found applications in high energy physics. A more recent publication in this series of about 20 papers is entitled "Universal constructions in Teichmiiller Theory" , Advances in Mathematics 98, pp. 143-215 (1993). Here it was shown how to generalize lambda lengths to provide coordinates on a certain space of homeomorphisms (i.e. continuous bijections whose inverses are also continuous) of the circle, where there is one generalized lambda length for each arithmetic arrow. In further recent work described in "The Lie algebra of homeomorphisms of the circle" , Advances in Mathematics 140, pp. 282-322 (1998), written jointly with Feodor Malikov, there is described a representation of the coordinate deformations of these generalized lambda lengths as explicit vector fields on the circle. These vector fields are called "normalized elementary vector fields" , and their study led to other vector fields on the circle called "fans" and "hyperfans" . These three families of vector fields correspond, respectively, to the arithmetic, fan, and modular wavelets described here via the identification of the vector field f(θ) β_g with the function f(θ); there was, however, no discussion of wavelets, sampling of data, approximation of functions, or any algorithms in the published literature.

Each arithmetic wavelet is once-continuously differentiable on the circle, compactly supported, and localized in space. Arithmetic wavelets are not compactly supported in frequency, but instead, the frequency profile of an arithmetic wavelet is given algebraically by the formula for c^A used in Step 2 in the calculation of Fourier transforms; a non-compactly supported localization in frequency follows directly from this. The asserted formula for c^A can be derived without much difficulty but in several cases directly from Figures 3a-3e integrating twice by parts the usual expression for Fourier coefficients using the fact ϋ_A is once-continuously differentiable. More difficult to derive is the formula given for the exponential functions e^τnθ in terms of arithmetic wavelets which was used in Step 1 of the calculation of inverse Fourier transforms; this subtle computation with lambda lengths is given in "The Lie algebra of homeomorphisms of the circle".

Similar remarks apply to fan and modular wavelets, but fan wavelets are only continuous (they are not differentiable), and modular wavelets are not even continuous.

In order to explain the failure of orthonormality of arithmetic wavelets and illuminate the regularization used in the arithmetic wavelet filter, fix some matrix A labeling an arithmetic arrow whose underlying chord has initial point q in the circle C. The sequence UⁿA of matrices label the arithmetic arrows whose underlying chords also have q as initial point, where n denotes an arbitrary integer. There is one infinite linear dependence

for each q e Q, and this explains the additive degree of freedom in the calculation of eu_A and er_A in the wavelet filter which is handled by regularization as was discussed before. This additive degree of freedom is also manifest in the initial specification e = 0 of wavelet coefficient on the doe in the wavelet filter.

In order to better understand the Farey quadrature, it is convenient to transform the disk D to the upper h plane U = {u + iv : v > 0} by means of ' transform

K : D → U

. 2 + 1

2 H J - . z — 1

This function K is a homeomorphism on the interior of D which maps the unit circle C bijectively to the real-axis {u + iv : v = 0} together with an additional point oo at infinity, where K(—l) = 0, K(—i) = 1, and K(+l) = oo. Furthermore, as q € Q varies over all quadrature points, the real number K(q) varies over all rational numbers, where the rational number corresponds to the quadrature point ?^ ² _ Q. In order to illustrate the Farey tesselation itself in U, whenever q\, q₂ £ Q are the endpoints of the chord underlying an arithmetic arrow in D, draw the semi-circle in U with real endpoints 9₁, - ₂; in the degenerate case that one of the quadrature points is infinite, say the point q₂ = 00, then draw the vertical ray in U with endpoint q_\ . Applying this procedure to each arithmetic arrow produces the classical model of the Farey tesselation in U that is indicated in Figure 11. This figure illustrates the relationship between the Farey tesselation and the indicated circle packing, where the circle in the figure that is tangent to the real-axis at the rational point p/q has diameter q^~2. For further information on number-theoretic aspects of this classical figure, see, for instance, Ford's book "Automorphic Functions" , Chelsea (1972).

Source Codes

In order to adequately disclose the methods, it is the purpose of this paragraph to include source codes for their implementations in the computer language C. Source code is presented for each of the following two implementations.

• awft is the implementation of a preferred embodiment of the method for computing the Fourier transform of real-valued input data defined on the circle. The algorithm depends upon a bandwidth N, tolerance EPS, generation cut-off NVAN, and minimum generation MING as described before. A further flag SMODE=l determines that the output data will include the 0,+l,-l Fourier coefficients, and SMODE=0 otherwise. (The only configuration in either code which requires any library routines is SMODE=l, which requires math.h.)

• awift is the implementation of a preferred embodiment of the method for computing the inverse Fourier transform of specified complex-valued Fourier coefficients defined in a specified bandwidth N. The output data is controlled by a single parameter SCAT, which is related to the parameter SCALE discussed before by SCAT=[exp sinh ^" (-1)(1/SCALE)] ^" 2. (SCAT varies inversely with SCALE and may be compared with conveniently derived approximate quantities.)

The comparison of the subroutines tgener in awft and in awift illustrates the easy transition between source codes for real and for complex data. The corresponding complex-awft and real-awift are therefore easily derived from the included source codes. Furthermore, the met J for compression of one-dimensional inj data amounts to the first step of awft, then quantization/de-quantization, then the second step of awift. The included source codes together thus also implement the method for compression of one-dimensional data since source code for it may be extracted from the included source codes. Finally, it is straight-forward to write the further driving routine required for multi-dimensional data compression, so the source codes included with this patent application are sufficient to readily implement the method for multi-dimensional data compression as well.

The included implementations have not been optimized. (For instance, straight-forward optimization has improved run-time by as much as 70% over the included awft code; on the other hand, optimization destroys intelligibility of the coding in relation to the underlying algorithm, and this explains the inclusion of unoptimized source code.) Another substantial improvement in run-time over the included source codes can be achieved by factoring so as to replace various manipulations of complex numbers with manipulations of integers as discussed before.

Moreover, for clarity, the included source codes do not take advantage of renormalization, but for completeness, to each of awft and awift is appended source code to replace the subroutines tgener and bgener by corresponding subroutines which do employ renormalization.

/ This is an implementation in C of a preferred embodiment //of the invention for the calculation of the Fourier transform

#include <stdio.h> #include <math.h> //math.h is required only for the calculation of 0,+ 1 ,- 1 Fourier coeffs #define GEN 1000 //the maximum allowed generation #define N 50 //output bandwidth #defιne COR 1 //irrelevant parameter #define EPS l.e-12 //tolerance, where contributions to Fourier coeffs of norm less than //EPS are regarded as negligible

#define PI 3.141592653589793

#define NVAN 1 //terminate the recursion when NVAN consecutive

//generations have contributed negligibly

#define MING 15 //calculate at least MING generations

#define SMODE 1 //SMODE=l includes the calculation of 0,+l,-l Fourier coeffs //SMODE=0 is faster and does not require math.h

struct complex{ /a complex number double x; double y;

struct trip{ //convenient to have a real triple double alp; double bet; double gam;

}; struct edge { //an arrow-structure int a,b,c,d; double e,alp,bet,gam;

int count=0; //a running count of the total number of samples of input data struct complex fourier[2*N+l]; //the complex Fourier coeffs stored in the order

//n=-N,...,-l,0,l....,N, so fourier(n) is (N+n)th coeff double abar,bbar,cbar; //coeffs in trig function abar^cos theta + bbar*sin theta + cbar //which agrees with input function for theta = -PI, -PI/2, 0 void main(void){ void genertop(struct edge *,int,int,double); //recursive routine for top of circle void generbot(struct edge *,int,int,double); //recursive routine for bottom of circle void normalout(void); //output normalization routine void normalin(void); //input normalization routine double fbar(int.int); //fbar(p.q) returns f(theta) - (abar*cos theta + bbar*sin theta + cbar) //where theta is the argument of the complex number (p-iq)/(p+iq) struct edge doe[l]; //initial edge for recursive calls int i; nor alinQ; //this calculates the normalization parameters abar, bbar, cbar doe->a=doe->d=l; //initialize the matrix of the doe to the identity doe->b=doe->c=0; doe->alp=doe->bet=doe->gam=doe->e=0; //initialize wavelet and lagged trig coeffs genertop(doe, 1 ,0,fbar(- 1,1)); //begin recursion on top of the circle generbot(doe, 1 ,0,0.) ; //begin recursion on bottom of the circle normaloutQ; //this divides each output Fourier coeff by an overall factor PI //which was suppressed during the recursion and determines //negative from positive Fourier coeffs using reality of input printfC'with EPS=%e, NVAN=%d, COR=%d, and MING=%d, the count is: %d\n" EPS,NVAN,COR,MING,count); for(i=0;i<=N;i++) printf("c%3d = %e+(i)%e\n",i,fourier[N+i].x,fourier[N+i].y);

void genertop(struct edge *p,int g.int envy .double fc){ //recursion for top of circle

//input pointer to a current edge p of generation g, where there //have so far been envy consecutive generations of negligible //contributions, and fc is the current normalized input sample struct edge twofer[2], // the descendant edges of p

struct complex temp; void genertop(struct edge *,int,int,double); //recursive function struct complex tgener(struct edge *,struct edge *,double);

//tgener(p,twofer,f) stuffs twofer with the descendants of p //using the updated normalized input sample f returning the //next two normalized samples as its real and imaginary parts int charles(struct edge *); //charles(p) updates the Fourier coeffs with the

//contribution from edge p returning 1 if they //are negligible, and returning 0 otherwise

void tprune(struct edge *); //this implements the processing of a terminal edge if(envy>=NVAN && g>= MING){ //if appropriate, then terminate cascade tprune(p); return;

temp=tgener(p,twofer,fc); //otherwise, generate descendants of p if(charles(twofer)*charles(twofer+l)) //if both descendants contribute negligibly, envy++; //then increment envy else envy=0; //otherwise reset envy to zero if(g>=GEN) { //problem: non-convergence after GEN generations printfC'fall through with gen=%d for a,b,c,d=%d, d,%d,%d and e=%e\n", g,p->a,p->b,p->c,p->d,p->e); return;

else{ //otherwise, continue recursion with descendants genertop(twofer,g+ 1 ,envy,temp.x); genertop(twofer+ 1 ,g+ 1 ,envy,temp.y);

void generbot(struct edge *p,int g.int envy.double fc){ //recursion for bottom of circle struct edge twofer[2]; //this is entirely analogous to genertop struct complex temp; void generbot(struct edge *,int,int,double); struct complex bgener(struct edge *,struct edge *,double); int charles(struct edge *); void bprune(struct edge *); if(envy>=NVAN && g >= MING){ bprune(p); return;

temp=bgener(p,twofer,fc); if(charles(twofer)*charles(twofer+l)) envy++; else envy=-0; if(g>=GEN){ printfC'fall through with gen=%d for a,b,c,d=%d,%d,%d,%d and e= e\n", g,p->a,p->b,p->c,p->d,p->e); return;

else{ generbot(twofer,g+l,envy,temp.x); generbot(twofer+ 1 ,g+ 1 ,envy,temp.y);

}

int charles(struct edge * p) ( //charles(p) updates the Fourier coeffs using p

//and returns a control parameter flag=l if int n.flag; //contribution is negligible and flag=0 otherwise struct complex makecpx(int.int); //input p,q returns the complex number (p-iq)/(p+iq) struct complex mult(struct complex.struct complex); //complex multiplication int a,b,c,d; //local variables for matrix values double e; //local variable for wavelet coeff struct complex pzl,pzr,pzf,pzb; //l,r,f,b labels left,right,front,back vertices

//of Farey tesselation near the edge p, and struct complex plm,prm,pfm,pbm; //in a loop below on index n, pzx=(pxm)^Λn //for x=l,r,f,b struct complex tmp,mul,cn; int sixupc(int,int,int,int,double); //this updates the Fourier coeffs 0,+l,-l //provided SMODE=l

a=p->a; //stuff local values b=p->b; c=p->c; d=p->d; e=p->e; pzr.x=b-d; //set-up for loop over Fourier coeffs pzr.y=c-a; prm=makecpx(b-d,a-c) ; pzr=mult(pzr,pzr);

pzl.x=b+d; pzl.y=-c-a; plm=makecpx(b+d,a+c); pzl=mult(pzl,pzl);

pzf.x=d; pzf.y=-c; pfm=makecpx(d,c); pzf=mult(pzf,pzf);

pzb.x=b; pzb.y=-a; pbm=makecpx(b,a); pzb=mult(pzb,pzb); flag=l; //initialize flag tmp.x=0; if(SMODE) flag=sixupc(a,b,c,d,e); //update 0,+l,-l modes and re-set flag to 0 if any

for(n=2;n<=N;n++){ //loop over positive Fourier coeffs n>l tmp.y=e/(n-n*n*n); pzr=mult(pzr,prm); //recursively calculate powers pzl=mult(pzl,plm); pzf=mult(μ. ,pfm); ρzb=muit(pzb,pbm) ;

cn.x=-pzr.x+2*pzf.x+2*pzb.x-pzl.x; cn.y=-pzr.y+2*pzf.y+2*pzb.y-pzl.y; mul=mult(cn,tmp); fourier[n+N].x+=mul.x; //update nth Fourier coeff fourier[n+N].y+=mul.y; if(flag==l && mul.x*mul.x+mul.y*mul.y>=EPS) //if contribution to nth Fourier flag=0; //coeff is non-negligible,

//then set flag to zero } return(flag);

int sixupc(int a.int b.int c,int d.double e){ //this updates the 0,+ 1,-1 Fourier coeffs //using the input matrix entries a,b,c,d //and wavelet coeff e double tp,tm,tr,tl; //p,m,r,I labels the terminal,initial,right,left vertices of Farey

//tesselation near edge, and tx denotes corresponding angle double dtemp; double tb; //tb=l on the top and tb=-l on the bottom of the circle int fl; //this is the local version of the flag in charles and is the //value returned to charles by sixupc struct complex makecpx(int,int), temp; fl=l; //set the flag to 1 tb=l.; //default to the top of the circle if(b+c<0) //if on the bottom of the circle tb=-l.; //then re-set tb=-l.

//calculate the various angles temp=makecpx(d,-c); tp=acos(tb*temp.x); temp=makecpx(b,-a);

temp=makecpx(b+d,-(a+c)); tl=acos(tb*temp.x); temp=makecpx(b-d,c-a); tr=acos(tb*temp.x); dtemp=e*( //contribution to 0th Fourier coeff

2.*tp*(c*c+d*d)+2.*tm*(a*a+b*b)

-tl*((a+c)*(a+c)+(b+d)*(b+d))

-tr*((a-c)*(a-c)+(b-d)*(b-d))

); if(dtemp*dtemp>=EPS) //if this contribution is non-negligible, fl=0; //then re-set the flag fourier[N].x+=dtemp; //up-date the 0th Fourier coeff temp.x=e*( //contribution to real part of 1st tp*(c*c-d*d)+tm*(a*a-b*b) //Fourier coeff

+tr*((b-d)*(b-d)-(a-c)*(a-c))/2. +tl*((b+d)*(b+d)-(a+c)*(a+c))/2.

); fourier[N+ 1 ].x+=temp.x; //update real part of 1st Fourier coeff

temp.y=e*( //contribution to imag part of 1st

2.*(tp*c*d4-tm*a*b) //Fourier coeff

-tr*(b-d)*(a-c)

-tl*(b+d)*(a+c)

); fourier[N+l).y+=temp.y; //update imag part of 1st Fourier coeff if(temp.x*temp.x+temp.y*temp.y>=EPS) //if this contribution is non-negligible, fl=0; //then re-set the flag to zero return(fl);

struct complex mult(struct complex u.struct complex v)( //complex multiplication struct complex temp; temp.x=u.x*v.x-u.j .y; temp.y=u.x*v.y+v.x*u.y; return(temp);

struct complex makecpx(int p.int q){ //input ρ,q returns the complex number (p+iq)/(p+iq) struct complex temp; double den; den=p*p+q*q; temp.x=(p*p-q*q)/den; temp.y=-2*p*q/den; return(temρ);

void normalout(void){ //output normalization double fz=0.,fι=0.,fo=0.; double alp,bet,gam; in n; for(n=0;n<=N;n++){ fourier[N+n].x=fourier[N+n].x/PI; //include factor PI suppressed before fourier[N+n] .y=fourier[N+n] .y/PI; fourier[N-n].x=fourier[N+n].x; //calculate negative from positve fourier[N-n].y=-fourier[N+n].y; //Fourier coeffs

)

void normalin(void){ //input normalization double fz=0.,fi=0.,fo=0.; //values of input function for angles -PI, -PI/2, 0 int i; double f(int,int); //f(p,q) returns the un-normalized input sample //at the angle corresponding to the point on the //circle given by the complex number (p-iq)/(p+iq) for(i=0;i<=2*N;i++) fourier[i].y=fourier[i].x=0; fz=f(0, l); fi=f(l,0); fo=f(l,l); abar=(fι-fz)/2; //abar*cos theta + bbar*sin theta + cbar takes the same values cbar=(fi+fz)/2; //as f for the angles theta=-PI, -PI/2, 0 bbar=cbar-fo; fourier[N].x=cbar*PI; //stuff the Oth Fourier coeff scaled by PI to fourier[N].y=0.; //recover it after re-scaling by 1/PI in normalout fourier[N+ 1 ] .x=abar*PI/2; //stuff the 1st Fourier coeff scaled by PI to fourier[N+ 1 ].y=-bbar*PI/2; //recover it after re-scaling by 1/PI in normalout

double fbar(int p, int q){ //input p,q to fbar produces the nornalized value struct complex temp; //of the input data at the point (p-iq)/(p+iq) struct complex makecpx(int.int); double f (int.int); temp=makecpx(p.q); return(f(p,q) - (abar*temp.x + bbar*temp.y+cbar));

} double f(int p,int q){ //f is the input data, where f(p,q) is the double fsinm(int,int,int); //value taken at (p-iq)/(p+iq) return(fsinm(p,q,6));

double fsinm(int p.int q, int m){ //for example, the input function is int n; //here taken to be sin(6*theta) struct complex zeta.meta; struct complex makecpx(int,int),mult(struct complex.struct complex); zeta=makecpx(p,q) ; meta.x=l.; meta.y=0; for(n=l ;n<=m;n++) meta=mult(zeta,meta); return(meta.y);

struct complex tgener(struct edge *pc,struct edge *pn,double fc){

//tgener(p, twofer.fc) stuffs twofer with the descendants //of p for the top of the circle, where fc is the current //normalized input sample, and tgener returns a complex nun-oer. //the real and imag parts being the normalized input values //required for the least-squares regularization double fbar(int,int); //fbar(p,q) returns the normalized input values int a,b,c,d,bpd,apc; //local variables for matrix entries, bpd=b+d, apc=a+c double e,alp,bet,gam; //local variables for wavelet and lagged trig coeffs double axl ,ayl,ax2,ay2; //index 1,2 refers to fιrst,second (counter-clockwise) double bxl,byl,bx2,by2; //descendants, prefix a,b,c refers to alp,bet,gam double ex 1 ,cy 1 ,cx2,cy2; //update procedures below define indices x,y double sigl,taul,sig2,tau2,chi; //quantities used in the least-squares regularization double eu.et; //wavelet coeffs of first and second descendant edges double fn.fnp; //normalized input samples

double vlc,vlu,vlt,v2c,v2u,v2t;

//coeffs in the inhomogeneous linear expression vic+viu*eu+vit*et, i=l,2, //of ongoing approximation to normalized input at the respective sample points //((2b+d)+i(2a+c))/((2b+d)-i(2a+c)), ((b+2d)+i(a+2c))/((b+2d)-i(a+2c))

struct complex temp; struct edge *pnp;

pnp=pn; //initialize pointers pnp++; a=pc->a; //initialize matrix entries b=pc->b; c=pc->c; d=pc->d; apc=a+c; bpd=b+d; pn->a=a; //stuff matrix of first descendant pn->b=b; pn->c=apc; pn->d=bpd; pnp->a=apc; //stuff matrix of second descendant pnp->b=bpd; pnp->c=c; pnp->d=d; alp=pc->alp; bet=pc->bet; gam=pc->gam; e=pc->e; //initialize wavelet coeff

temp.x=fn=fbar(-(b+bpd),a+apc); //store new normalized input samples as the temp.y=fnp=fbar(-(d+bpd),c+apc); //real and imag parts of temp to return count+=2; //update count of total sample points if(a>c){ //prepare for stuffing alp,bet,gam in case a>c axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; ayl=2*(b*b-a*a); bxl=bet+(c*d-a*d-b*c-a*b)*4*e; byl=4*a*b; cxl=gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; cyl=-2*(a*a+b*b); ax2=alp+4*(d*d-c*c)*e; ay2=2*(c*c-d*d); bx2=bet+8*e*c*d; by2=-4*c*d; cx2=gam-4*e*(c*c+d*d); cy2=2*(c*c+d*d); chi=2*e+( //chi=eu-et

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam)

-alp*((b+d)*(b+d)-(a+c)*(a+c))

-bet*2*(a+c)*(b+d)

)/4.;

} else{ //prepare for stuffing alp,bet,gam in case o=a

axl=alp+4*e*(a*a-b*b); ayl=2*(b*b-a*a); bxl=bet-8*e*a*b; byl=4*a*b; cxl=gam+4*e*(a*a+b*b); cyl=-2*(a- a+b*b); ax2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; ay2=2*(c*c-d*d); bx2=bet+(a*d-a*b+b*c4-c*d)*4*e; by2=-4*c*d; cx2=gam+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e; cy2=2*(c*c+d*d);

chi=-2*e+( //chi=eu-et

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam)

-alp*((b+d)*(b+d)-(a+c)*(a+c))

-bet*2*(a+c)*(b+d)

)/4.;

//in any case, calculate the coeffs vic,viu,vit, for i=l,2 vlc=-fn+cxl+

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl

+2*(b+bpd)*(a+apc)*bxl)/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); vlt=cyl+

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*ayl

+2*byl*(b4-bpd)*(a+apc))/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); vlu=8./((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); v2c=-fnp+cx2+

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ax2

+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); v2u=cy2+

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ay2

+2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); v2t=-8./((d+bρd)*(d+bpd)+(c+apc)*(c+apc)); sigl=vlc+chi*vlu; taul=vlt+vlu; sig2=v2c+chi*v2u; tau2=v2t+v2u; et=-(sigl*taul+sig2*tau2)/(tau l*taul+tau2*tau2); //calculate et and eu eu=chi+et; pn->e=eu; //stuff new wavelet coeffs pnp->e=et; pn->alp=ax 1 +ay 1 *et; //stuff first lagged trig coeffs pn->bet=b l+byl*et; pn->gam=cx 1 +cy 1 *et; pnp->alp=ax2+ay2*eu; //stuff second lagged trig coeffs pnp->bet=bx2+by2 *eu ; pnp->gam=cx2+cy2*eu; return (temp);

struct complex bgener(struct edge *pc,struct edge *pn,double fc){

//bgener(p,twofer,fc) stuffs twofer with the descendants //of p for the bottom of the circle, where fc is the current //normalized input sample and returns a complex number, //the real and imag parts being the normalized input values //required for the least-squares regularization

//this is entirely analogous to tgener, and only the differences //will be noted here double fbar(int.int); int a,b,c,d,bpd,apc; double e,alp,bet,gam; double axl,ayl,ax2,ay2; double bxl, by I,bx2,by2; double cxl,cyl,cx2,cy2; double sigl,taul,sig2,tau2,chi; double vlc,vlu,vlt,v2c,v2u,v2t; double eu,et,fn,fnp; struct complex temp; struct edge *pnp; pnp=pn; pnp++; a=pc->d; //differs from tgener in that d,-c,-b,a replaces a,b,c, b=-(pc->c); c=-(pc->b); d=pc->a; apc=a+c; bpd=b+d; pn->a=bpd; //differs from tgener in that d,-c,-b,a replaces a,b,c,d ρn->b=-apc; pn->c=-b; pn->d=a; pnp->a=d; //differs from tgener in that d,-c,-b,a replaces a,b,c,d ρnp->b=-c; pnp->c=-bpd; pnp->d=apc; alp=pc->alp; bet=pc->bet; gam=pc->gam;

e=pc->e;

temp.x=fn=fbar(a+apc,b+bpd); //differs from tgener in that temp.y=fnp-=fbar(c+apc,d+bpd); //(-q.p) replaces (p,q) count+=2; if(a>c){ axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; ayl=2*(b*b-a*a); bx l=bet+(c*d-a*d-b*c-a*b)*4*e; byl=4*a*b; cxl=gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; cyl=-2*(a*a+b*b); ax2=alp+4*(d*d-c*c)*e; ay2=2*(c*c-d*d); bx2=bet+8*e*c*d; by2=-4*c*d; cx2=gam-4*e*(c*c+d*d); cy2=2*(c*c+d*d); chi=2*e+(

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam)

-alp*((b+d)*(b+d)-(a+c)*(a+c))

~bet*2*(a+c)*(b+d)

)/4.;

else{ axl=alp+4*e*(a*a-b*b); ayl=2*(b*b-a*a); bxl=bet-8*e*a*b; byl=4*a*b; ex 1 =gam+4*e*(a*a+b*b); cyl=-2*(a*a+b*b); ax2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; ay2=2*(c*c-d*d); bx2=bet+(a*d-a*b+b*c+c*d)*4*e; by2=-4*c*d; cx2=gam+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e; cy2=2*(c*c+d*d);

chi=-2*e+(

((a+c)*(a+c)+(b+d)*(b+d))*(fc-gam)

-alp*((b+d)*(b+d)-(a+c)*(a+c))

-bet*2*(a+c)*(b+d)

)/4.;

vlc=-fn+cxl+

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl

+2*(b+bpd)*(a+apc)*bxl)/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); vlt=cyl+

(((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*ayl

+2*by 1 *(b+bpd)*(a+apc))/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); vlu=8./((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); v2c=-fnp+cx2+

(((d+bpd)*(d+bpd)-(c+apc)*(c+apc))*ax2

+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); v2u=cy2+

(((d+bpd)*(d4-bpd)-(c+apc)*(c+apc))*ay2

+2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); v2t=-8./((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); sigl=vlc+chi*vlu; taul=vlt+vlu; sig2=v2c+chi*v2u; tau2=v2t+v2u; et=-(sig 1 *tau 1 +sig2*tau2)/(tau 1 *tau 14-tau2*tau2); eu=chi+et; pn->e=eu; pnp->e=et; pn->alp=ax 1 +ay 1 *et; ; pn->bet=bx 1+by 1 *et; pn->gam=cx 1 +cy 1 *et; pnp->alp=ax2+ay2*eu; pnp->bet=bx2+by2*eu; pnp->gam=cx2+cy2*eu; return (temp);

void tprune(struct edge *p){ //this routine implements the processing of terminal //edge p in the cascade for the top of the circle double alp, bet, gam.e; //local copies of lagged trig coeffs and wavelet coeff double alpl,alp2,betl,bet2,gaml,gam2; //trig coeffs of the trigonometric function sigma at

//the first and second (counter-clockwise) samples struct trip mob(struct complex.struct complex, struct complex, double,double,double); //mob(zl,z2,z3,fl,f2,f3) returns the triple with entries alp,bet,gam //where alp*cos theta + bet*sin theta + gam takes the respective values // fl,f2,f3 at the points zl,z2,z3 on the circle given as complex numbers double fbar(int,int); //returns normalized input values void fixup(struct trip,int,int,int,int,int); //adjusts the output Fourier coeffs using sigma-tau struct complex makecpx(int,int); //input p,q returns (p-iq)/(p+iq) void sixupp(struct trip,int,int,int,int,int); //adjusts the 0,+l,-l Fourier coeffs provided //SMODE=l int a,b,c,d; //local variables for the matrix entries struct trip tripl,trip2; //first and second trig triples a=p->a; //stuff local variables b=p->b; c=p->c; d=p->d; e=p->e; alp=p->alp; bet=p->bet; gam=p->gam; if (a>c){ //final update of trig fields in case a>c alpl=alp+e*2*(a*(a+c)+d*(d-b)-b*(b+d)-c*(c-a)); betl=bet+e*4*(c*(d-b)-a*(d+b)); gaml=gam+e*2*(a*(a+c)+c*(a-c)+b*(b+d)-d*(d-b)); alp2=alp+e*4*(d-c)*(d+c); bet2=bet+e*8*c*d; gam2=gam-e*4*(c*c+d*d);

elsef //final update of trig fields in case o=a alp 1 =alp+e*4*(a-b)*(a+b) ; betl=bet-e*8*a*b; gaml=gam+e*4*(a*a+b*b);

alp2=alp+e*2*(a*(a-c)+b*(d-b)-c*(c+a)+d*(b+d)); bet2=bet+e*4*(a*(d-b)+c*(b+d)); gam2=gam+e*2*(a*(a-c)-c*(a+c)-b*(d-b)-d*(b+d));

tripl=mob( //calculate sigma for first descendant makecpx(-(3*b+d),3*a+c), makecpx(-(2*b+d),2*a+c), makecpx(-(3*b+2*d),3*a+2*c), fbar(-(3*b+d),3*a+c), fbar(-(2*b+d),2*a+c), fbar(-(3*b+--*d),3*a+2*c)

); trip2=mob( //calculate sigma for second descendant makecpx(-(2*b+3*d),2*a+3*c), makecpx(-(b+2*d),a+2*c), makecpx(-(b+3*d),a+3*c), fbar(-(2*b+3*d),2*a+3*c), fbar(-(b+2*d),a+2*c), fbar(-(b+3*d),a+3*c)

); count+=6; //update count of total sample points tripl.alp=(tripl.alp-alpl)/COR; //difference of first trig functions trip 1.bet=(trip 1.bet-bet 1 )/COR; trip 1.gam=(trip 1.gam-gam 1 )/COR; trip2.alp=(trip2.alp-alp2)/COR; //difference of second trig functions trip2.bet=(trip2.bet-bet2)/COR; trip2.gam=(trip2.gam-gam2)/COR; fixup(tripl,b+d,a+c,b,a,l); //adjust Fourier coeffs with first trig function difference fixup(trip2,d,c,b+d,a+c,l); //adjust Fourier coeffs with second trig function difference if(SMODE){ //if SMODE=l, then adjust 0,-1-1,-1 Fourier coeffs as well sixupp(trip 1 ,b+d,a+c,b,a, 1 ); sixupp(trip2,d,c,b+d,a+c,l); }

void bprune(struct edge *p){ //this routine implements the processing of terminal

//edge p in the cascade for the bottom of the circle

//this is entirely analogous to tgener, and only the //differences will be noted here double alp, bet, gam,e; double alp 1 ,alp2,bet 1 ,bet2,gam 1 ,gam2; struct trip mob(struct complex.struct complex, struct complex, double.double.double); double fbar(int.int); void fixup(struct trip, int,int,int,int,int); struct complex makecpx(int.int); void sixupp(struct trip,int,int,int,int,int); int a,b,c,d; struct trip tripl,trip2; a=p->d; //differs from tprune in that d,-c,-b,a replaces a,b,c,d b=-(p->c); c=-(p->b); d=p->a; e=p->e; alp=p->alp; bet=p->bet; gam=p->gam; if (a>c){ alpl=alp+e*2*(a*(a+c)+d*(d-b)-b*(b+d)-c*(c-a)); betl=bet+e*4*(c*(d-b)-a*(d+b)); gaml=gam+e*2*(a*(a+c)+c*(a-c)+b*(b+d)-d*(d-b)); alp2=alp+e*4*(d-c)*(d+c); bet2=bet+e*8*c*d; gam2=gam-e*4*(c*c+d*d);

else{ alpl=alp+e*4*(a-b)*(a+b); betl=bet-e*8*a*b; gaml=gam+e*4*(a*a+b*b); alp2=alp+e*2*(a*(a-c)+b*(d-b)-c*(c+a)+d*(b+d)); bet2=bet+e*4*(a*(d-b)+c*(b+d)); gam2=gam+e*2*(a*(a-c)-c*(a+c)-b*(d-b)-d*(b+d));

tripl=mob( //arguments of mob differ from those in tprune makecpx(3*a+c,3*b+d), makecpx(2*a+c,2*b+d), makecpx(3*a4-2*c,3*b+2*d), fbar(3*a+c,3*b+d), fbar(2*a+c,2*b+d), fbar(3*a+2*c,3*b+2*d)

); trip2=mob( //arguments of mob differ from those in tprune makecpx(2*a+3*c,2*b+3*d), makecpx(a+2*c,b+2*d), makecpx(a+3*c,b+3*d), fbar(2*a+3*c,2*b+3*d), f ar(a+2*c,b+2*d), fbar(a+3*c,b+3*d)

); count+=6; tr-pl.alp=(tripl.alp+alpl)/COR; //differs from tprune in that alp.bet tripl.bet=(tripl.bet+betl)/COR; //is replaced by -alp,-bet trip 1.gam=(tri p 1.gam-gam 1 )/COR; trip2.alp=(trip2.alp+alp2)/COR; //differs from tprune in that alp,bet triρ2.bet=(triρ2.bet+bet2)/COR; //is replaced by -alp,-bet trip2.gam=(trip2.gam-gam2)/COR; fixup(tripl,b+d,a+c,b,a,-l); //last argument of fixup differs from tprune fιxuρ(trip2,d,c,b+d,a+c,-l);Q //last argument of fixup differs from tprune if(SMODE){ //if SMODE=l, adjust the 0,+l,-l Fourier coeffs sixupp(tripl,b+d,a+c,b,a,-l); //last argument of sixupp differs from tprune sixupp(trip2,d,c,b+d,a+c,-l); //last argument of sixupp differs from tprune

struct trip mob(struct complex zl, struct complex z2,struct complex z3, double fl,double f2,double D){ //mob(zl,z2,z3,fl,f2,f3) returns the triple with entries alp,bet,gam //where a!p*cos theta + bet*sin theta + gam takes the respective values // fl,f2,f3 at the points zl,z2,z3 on the circle given as complex numbers double det,aa,bb,cc,dd; struct trip temp; aa=(zl.x-z3.x); bb=(zl.y-z3.y); cc=(z2.x-z3.x); dd=(z2.y-z3.y); det=aa*dd-bb*cc; //general principles guarantee that det is non-zero aa=aa det; bb=bb/det; cc=cc/det; dd=dd/det; temp.alp=dd*(fl-f3)-bb*(f2-f3); temp.bet=aa*(f2-f3)-cc*(f 1 -f3); temp.gam=f3-z3.x*temp.alp-z3.y*temp.bet; return(temp); void fixup(struct trip delt,int dm.int cm.int dp.int cp.int tb){

//fixup delt,dm,cm,dp,cp,tb) adjusts the Fourier coeffs by adding the Fourier //coeffs of delt=T^ΛC-T with support truncated to be the interval with endpoints //(dm+i*cm)/(dm-i*cm), (dp+i*dp)/(dp-i*cp)

//let tm, tp be the corresponding angles

// tb=+l on the top, and tb=-l on the bottom of the circle int n; double cpp,cmp,cpm,cmrn; //in loop below on n, cpx=[cos(l+n)tx]/(n+l), cmx=[cos(l-n)tx]/(l-n), for x=p,m double spm,smm, spp.smp; //in loop below on n, spx=[sin(l+n)tx]/(l+n), smx=[sin(l-n)tx]/(l-n), for x=p,m double crm,crp,srm,srp; //in loop below on n, crm=[cos n*tx]/n, sπτι= [sin n*tx]/n, for x=p,m double alp,bet,gam; //local variables for entries of delt double ctp, ctm, cnp.cnm; //in loop below on n, ctx=cos tx, cnx=cos n*tx, for x=p,m double stm,snm,stp,snp; //in loop below on n, stx=sin tx, snx=sin n*tx, for x=p,m double cbufp,sbufp,cbufm,sbufm; struct complex makecpx(int,int),temp;

alp=delt.alp; //stuff local variables bet=delt.bet; gam=delt.gam; temp=makecpx(dm,-cm); ctm=tb*temρ.x; stm=tb*temρ.y; temp=makecpx(dp,-cp) ; ctp=tb*temp.x; stp=tb*temp.y; cnp=ctp*ctp-stp*stp; snp=2*stp*ctp; cnm=ctm*ctm-stm*stm; snm=2*stm*ctm; for(n=2;n<=N;n++){ //loop over Fourier coeffs cbufp=ctp*cnp-stp*snp; //update using formulas for cos sbufp=stp*cnp+ctp*snp; //or sin of a sum of angles cbufm=ctm*cnm-stm*snm; sbufm=stm*cnm+ctm*snm; cpp=cbufp/( 1 +n) ; //update for p cmp=(cbufp+2*stp*snp)/(l-n); spp=sbufp/(l+n); smp=(sbufp-2*ctp*snp)/(l-n); cpm=cbufm/(l+n); //update for m cmm=(cbufm+2*stm*snm)/(l-n); spm=sbuf m/( 1 +n) ; smm=(sbufm-2*ctm*snm)/( 1 -n) ; crp=2.*cnp/n; crm=2.*cnm/n; srp=2.*snp/n; srm=2.*snm/n; cnp=cbufp; //update for next iteration snp=sbufp; cnm=cbufm; snm=sbufm; fourier[N+n].x+=( //adjust real part of nth Fourier coeff

+alp*(smp+spp-smm-spm)

+bet*(-cpp-cmp+cpm+cmm)

+gam*(srp-srm)

)/4.; fourier[N+n].y+=( //adjust imag part of nth Fourier coeff

+alp*(cpp-cmp-cpm+cmm)

+bet*(spp-smp-spm+smm)

+gam*(c -crm)

)/4.;

} return;

void sixupp(struct trip triple.int dm.int cm.int dp,int cp,int tb)(

//sixup(delt,dm,cm,dp,cp,tb) adjusts the 0.+1.-1 Fourier coeffs by adding the 0,-4-1,-1 Fourier //coeffs of delt=T^ΛC-T with support truncated to be the interval with endpoints //(dm+i*cm)/(dm-i*cm). (dp+.-"dp)/(dp-i*cρ)

//let tm, tp be the corresponding angles

// tb=+l on the top, and tb=-l on the bottom of the circle

double alp,bet,gam; //local variables for entries of delt double com,sim,cop,sip,tm,tp; //cox=cos tx, six=sin tx, for x=p,m struct complex temp; struct complex makecpx(int.int);

temp=makecpx(dm,-cm); //calculate com,sim,tm com=tb*temp.x; sim=tb*temp.y; tm=acos(com); temp=makecpx(dp,-cρ) ; //calculate cop.sip.tp cop=tb*temp.x; sip=tb*temp.y; tp=acos(cop);

alp=tb*triple.alp; //adjust alp.bet for the bottom of the circle bet=tb*triple.bet; gam=triple.gam; fourier[N].x+=( //update 0th Fourier coeff gam*(tp-tm)+alp*(sip-sim)-bet*(cop-com) )/2.; fourier[N+l].x+=( //update real part of 1st Fourier coeff alp*(tp-tm)/2. +gam*(sip-sim)

-bet*(cop*cop-sip*sip-com*com+sim*sim)/4. +alp*(cop*sip-com*sim)/2. )/2.;

fourier[N+l].y+=( //update imag part of 1st Fourier coeff

-bet*(tp-tm)/2. +gam*(cop-com)

+alp*(cop*cop-sip*sip-com*com+sim*sim)/4. +bet*(cop*sip-com*sim)/2. )/2.; //Here are the subroutines replacing tgener and bgener above in an //implementation which takes advantage of renormalization

struct complex tgener(struct edge *pc,struct edge *pn,doub!e fc)

double fbar(int.int); int a,b,c,d,bpd,apc; double e,alp,bet,gam; double axl,ayl,ax2,ay2; double bxl,byl,bx2,by2; double cxl,cyl,cx2,cy2; double sigl,taul,sig2,tau2,chi; double vlc,vlu,vlt,v2c,v2u,v2t; double deriv(int,int,int,int,int,int); double dc,dn,dnp; double eu,et,fn,fnp; struct complex temp; struct edge *pnp; extern int count; pnp=pn; pnp++;

a=pc->a; b=pc->b; c=pc->c; d=pc->d; apc=a+c; bpd=b+d; pn->a=a; pn->b=b; pn->c=apc; pn->d=bpd; pnp->a=apc; pnp->b=bpd; pnp->c=c; pnp->d=d; alp=ρc->alp; bet=ρc->bet; gam=pc->gam; e=pc->e;

temp.x=fn=fbar(-(b+bpd),a+apc); temp.y=fnp=fbar(-(d+bpd),c+apc); temp.n=0; count+=2; dc=deri v(a,b,c,d,- 1,1); dn=deriv(a,b,c,d,-l,2); dnp=deriv(a,b,c,d,-2, 1); fc=fc/dc; fn=fn dn; fnp=fnp/dnp;

if(a>c){ a l=alp+4.*e; bxl=bet-4.*e; cxl=gam; ax2=alp+4.*e; bx2=bet; cx2=gam-4.*e; chi=2.*((fc-gam-bet)/4.+e);

else{

axl=alp+4.*e; bxl=bet; cxl=gam+4.*e; ax2=alp+4.*e; bx2=bet+4.*e; cx2=gam; chi=2. *((fc-gam-bet)/4.-e) ;

vlc=-fn+cxl+(4.*bxl-3.*axl)/5.; vlt=-4./5.; vlu=8./5.; v2c=-fnp+cx24-(4.*UAi+3.*ax2)/5.; v2u=4./5.; v2t=-8./5.; sigl=vlc+chi*vlu; taul=vlt+vlu; sig2=v2c+chi*v2u; tau2=v2t+v2u; et=-(sigl*taul*(l)+sig2*tau2*(l))/(taul *taul*(l) +tau2*tau2*(l)); eu=chi+et; pn->e=eu; pnp->e=et; alp=axl-2.*et; bet=bxl; gam=cxl-2.*et; pn->alp=(2.*bet+alp+gam)/2.; pn->bet=(bet-alp+gam); pn->gam=(2.*bet-alp+3.*gam)/2.; alp=ax2-2.*eu; bet=bx2; gam=cx2+2.*eu; pnp->alp=(-2. *bet+alp-gam)/2. ; pnp->bet=alp+bet+gam; pnp->gam=(2.*bet+alp+3.*gam)/2.;

return (temp);

} double deriv(int a.int b.int c.int d,int p.int q){ double x,y,nep; x=p*p-q*q; y=-2*p*q; x=x/(p*p+q*q); y=y/(p*p+q*q); nep=-(a*a+b*b+c*c-κl*d)+x*(a*a+b*b-c*c-d*d)-y*(2*a*c+2*b*d); return(-2./nep); struct complex bgener(struct edge *pc,struct edge *pn,double fc) double fbar(int,int); int a,b,c,d,bpd,apc; double e,alp,bet,gam; double axl,ayl,ax2,ay2; double bxl,byl,bx2,by2; double ex 1 ,cy 1 ,cx2,cy2; double sigl,taul,sig2,tau2,chi; double vlc,vlu,vlt,v2c,v2u,v2t; double eu,et,fn,fnp; extern int count; double deriv(int,int,int,int,int,int); double dc,dn,dnp; struct complex temp; struct edge *pnp; pnp=pn; pnp++; a=pc->d; b=-(pc->c); c=-(pc->b); d=pc->a; apc=a+c; bpd=b+d; pn->a=bpd; pn->b=-aρc; pn->c=-b; pn->d=a; pnp->a=d; pnp->b=-c; pnp->c=-bpd; pnp->d=apc; alp=pc->alp; bet=pc->bet; gam=pc->gam; e=ρc->e; temp.x=fπ=fbar(a+apc,b+bpd); temp.y=fnρ=fbar(c+apc,d+bpd); temp.n=0; count+=2; dc=deri v(d,-c,-b,a, 1,1); dn=deriv(d,-c,-b,a,2, 1); dnp=deri v(d,-c,-b,a, 1 ,2); fc=fc/dc; fn=fn/dn; fnp=fnp/dnp; if(a>c){ axl=alp+4.*e; bxl=bet-4.*e; cxl=gam; ax2=alp+4.*e; bx2=bet; cx2=gam-4.*e; chi=2.*((fc-gam-bet)/4.+e);

else{ axl=alp+4*e; bxl=bet; cxl=gam+4*e; ax2=alp+4*e; bx2=bet+4*e; cx2=gam; chi=2*((fc-gam-bet)/4-e);

vlc=-fn+cxl+(4.*bxl-3.*axl)/5.; vlt=-4./5.; vlu=8./5.; v2c=-fnp+cx2+(4. *bx2+3. *ax2)/5. : v2u=4./5.; v2t=-8./5.; sigl=vlc+chi*vlu; taul=vlt+vlu; sig2=v2c+chi*v2u; tau2=v2t+v2u; et=-(sigl *taul*(l) . _-g2*-au2*(l))/(taul *-aul *(l) +tau2*tau2*(l)); eu=chi+et; pn->e=eu; pnp->e=et; alp=axl-2.*et; bet=bxl; gam=cxl-2.*et; pn->alp=(2.*bet+alp+gam)/2.; pn->bet=(bet-alp+gam) ; pn->gam=(2. *bet-alp+3. *gam)/2. ; alp=ax2-2.*eu; bet=bx2; gam=cx2+2.*eu; pnp->alρ=(-2. *bet+alp-gam)/2. ; pnp->bet=alp+bet+gam; pnp->gam=(2.*bet+alp+3.*gam)/2.; return (temp);

/ This is an implen._.-tation in C of a preferred embodiment

//of the invention for the calculation of the inverse Fourier transform

#include <stdio.h> #define SCAT 50 //SCAT=[exp sinh^Λ(-l)(l/SCALE)]^Λ2 #define N 50 //input bandwidth #define MINMODE 0 //minimum frequency of input #defιne PI 3.141592653589793 #defιne CHUNK 5000 //chunk of memory allocated as required

struct complex! //a complex number double x; double y;

struct edge ( //a complex arrow-structure int a,b,c,d; struct complex e,alp,bet,gam;

}; struct pqxy{ //the basic data type of the output, where int p,q; //p,q corresponds to the complex number (p-iq)/(p+iq) double x,y; //at which the output function takes the complex value x+iy

struct complex pfourier[2*N+l]; struct complex *fourier=pfourier+N; //the input Fourier coeffs //indexed -N,...,-1,0,1,...,N struct complex czer.cone.cmon; //the modified 0,+l,-l Fourier coeffs struct pqxy graf[CHUNK]; //output spatial values in

//clockwise-order around the circle //starting from the complex number -1 int outcount=0; //running tally of the total //number of output values

void main(void){ void genertop(struct edge *,int); //recursive routine for top of circle void generbot(struct edge *,int); //recursive routine for bottom of circle void normalout(void); //output normalization routine void normalin(void); //input normalization routine

struct complex cogi(struct edge *); //returns complex wavelet coeff of argument struct complex ei,emu,emt,eu-et: //wavelet coeffs on I,U^Λ{-1 },T^Λ{ -1 ) ,U,T struct complex emtt,emtu,emut,emuu; //wavelet coeffs on T^Λ{-2),

//T^Λ{-l }U^Λ{-l },U^Λ{-l }T^Λ{-l ),U^Λ{-2} struct edge startf 1]; //initial edge for each of //the various recursive calls int i;

//begin input data for(i=0;i<=2*N;i++){ pfourier[i].x=0.; pfourier[i].y=0.;

fourier[-5].x=.5;

//end input data

normalinO; //this stuffs czer,cone,cmon with the //modified 0,-1-1,-1 Fourier coefficients graf[outcount].p=0; //initial normalized output value graffoutcount] .q= 1 ; graf[outcount].x=0.; graf[outcount].y=0.; outcount++;

//calculate the wavelet coeffs for the various required edges of small generation

start->a=start->d= 1 ; start->c=0; start->b=-l; emt=cogi(start); start->a=start->d=l ; start->c=0; start->b=l; et=cogi(start); start->a=start->d=l ; start->c=-l; start->b=0; emu=cogi(start); start->a=start->d= 1 ; start->c=l; start->b=0; eu=cogi(start); start->a=start->d= 1 ; ei=cogi(start); stan->a-=start->d=l ; start->b=0; start->c=-2; emuu=cogi(start); start->a=start->d= 1 ; start->b=-2; start->c=0; emtt=cogi(start); start->a=l ; start->b=start->c=- 1 ; start->d=2; emut=cogi(start); start->a=2; start->d=l; start->b=start->c=- 1 ; emtu=cogi(start);

//serially perform the recursions for the initial edges in counter-clockwise order

start->a=start->d=l; //initialize start->c-=l ; start->b=0; start->e=eu; start->alp.x=2.*ei.x-2.*et.x; start->alp.y=2.*ei.y-2.*et.y; start->bet.x=2.*(emt.x-emu.x)-2.*ei.x; start->bet.y=2.*(emt.y-emu.y)-2.*ei.y; start->gam.x=2.*ei.x-2.*et.x; start->gam.y=2.*ei.y-2.*et.y; genertop(start,l); //recursion for U start->a=start->d=l; //initialize start->c=0; start->b=l; start->e=et; start->alp.x=2.*ei.x-2.*eu.x; start->alp.y=2.*ei.y-2.*eu.y; start->bet.x=2.*(emt.x-emu.x)+2.*ei.x; start->bet.y=2.*(emt.y-ernu.y)4-2.*ei.y; start->gam.x=-2.*ei.x+2.*eu.x; start->gam.y=-2.*ei.y+2.*eu.y; genertop(start,l); //recursion for T start->a=start- d=l; //initialize start->c=0; start->b=-2; start->e=emtt; start->alp.x=-(-2.*ei.x+2.*emu.x-4.*emt.x+2.*emut.x); start->alp.y=-(-2.*ei.y+2.*emu.y-4.*emt.y+2.*emut.y); start->bet.x=-(2.*e....-2.*emu.x+2.*emt.x); start->bet.y=-(2.*ei.y-2.*emu.y+2.*emt.y); start->gam.x=2.*ei.x-2.*emu.x+4.*emt.x-2.*emut.x; start->gam.y=2.*ei.y-2.*emu.y+4.*emt.y-2.*emut.y; generbot(start,l); //recursion for T^Λ { -2 } start->a=l; //initialize start->d=2; start->b=start->c=- 1 ; start->e=emut; start->alp.x=-(-2.*ei.x+2.*emu.x+2.*emt.x); start->alp.y=-(-2.*ei.y+2.*emu.y+2.*emt.y); start->bet.x=-(2.*ei.x-2.*emu.x-6.*emt.x+4.*emtt.x); start->bet.y=-(2.*ei.y-2.*emu.y-6.*emt.y+4.*emtt.y); start->gam.x=2.*ei.x-2.*emu.x-6.*emt.x+4.*emtt.x; start->gam.y=2.*ei.y-2.*emu.y-6.*emt.y+4.*emtt.y; generbot(start,l); //recursion for U^Λ{-1 }T^Λ{-1 }

start->a=2; //initialize start->d=l; start->b=start->c=- 1 ; start->e=emtu; start->alp.x=-(-2.*ei.x+2.*emu.x+2.*emt.x); start->alp.y=-(-2.*ei.y+2.*emu.y+2.*emt.y); start->bet.x=-(-2.*ei.x+6.*emu.x+2.*emt.x-4.*emuu.x); start->bet.y=-(-2.*ei.y+6.*emu.y+2.*emt.y-4.*emuu.y); start->gam.x=-2.*ei.x+6.*emu.x+2.*emt.x-4.*emuu.x; start->gam.y=-2.*ei.y+6.*emu.y+2.*emt.y-4.*emuu.y; generbot(start,l); //recursion forT^Λ{-l }U^Λ{-1 } start->a=start->d=l; //initialize start->b=0; start->c=-2; start->e=emuu; start->alp.x=-(-2.*ei.x-4.*emu.x+2.*emt.x+2.*emtu.x); start->alp.y=-(-2.*ei.y-4.*emu.y+2.*emt.y+2.*emtu.y); start->bet.x=-(-2.*ei.x-2.*emu.x4-2.*emt.x); start->bet.y=-(-2.*ei.y-2.*emu.y+2.*emt.y); start->gam.x=-2.*ei.x-4.*emu.x+2.*emt.x+2.*emtu.x; start->gam.y=-2.*ei.y-4.*emu.y+2.*emt.y+2.*emtu.y; generbot(start,l); //recursion for U^Λ{-2}

//recursions complete, and graf is stuffed with normalized output data in counter-clockwise order normaloutO; //un-normalize by altering

//each output value using the //modified 0,4-1,-1 Fourier coeffs

cmon.x,cmon.y,czer.x,czer.y,cone.x,cone.y); for (i=0;i<outcount;i++) printf("(%-. -.d) %e+(i)%e\n", graf[i].p,graf[i].q,graf[i].x,graf[i].y);

void genertop(struct edge *p,int g)( //recursion for top of the circle struct edge twofer[2],*ped; //recursion generates twofer[2] from p, where

//the entries of twofer are the descendants of p struct complex temp; struct complex makecpx(int.int); //input p,q returns (p-iq)/(p+iq) double ct,st; //cos.sin for various angles void genertop(struct edge *,int); //recursive function void tgener(struct edge *,struct edge *); //tgener(p,twofer) stuffs twofer[2] //with the descendants of p int subtend(struct edge *),pp,qq; //returns 1 if both argument edges //subtend angles less than SCALE, //otherwise returns 0 tgener(p,twofer); //tgener(p,twofer) stuffs twofer //with the two descendants of p if(subtend(twofer)) { //if both descendants subtend //angles less than SCALE, then //write two output values and //terminate the recursion ped=twofer; pp=-(ped->d); //stuff output values for the first qq=ped->c; //(counter-clockwise) sample point graffoutcount] .p=pp; graf[outcount].q=qq; temp=makecpx(pp,qq); ct=temp.x; st=temp.y; graffoutcount]. x=ped->alp.x *ct+ped->bet.x*st

+ped->gam.x+(ped- >e.x)*4./(pp*pp-K}q*qq); graffoutcount] .y=ped->alp.y ^; *ct+ped->bet.y*st

+ped->gam.y+(ped- >e.y)*4./(pp*pp+qq*qq); outcount++; //update number of output values ped++; pp=-(ped->d); //stuff output values for the second qq=ped->c; //(counter-clockwise) sample point graffoutcount] .p=pp; graf[outcount].q=qq; temp=makecpx(pp,qq); ct=temp.x; st=temp.y; graf[outcount].x=ped->alp.x*ct+ped->bet.x*st

+ped->gam.x; graf[outcount].y=ped->alp.y*ct+ped->bet.y*st

+ped->gam.y; outcount++; //update number of output values return; //terminate the recursion

genertop(twofer,g+ 1 ); //in the contrary case that one of genertop(twofer+l,g - ; //the descendant edges subtends -. //large angle, then continue the //recursion in counter-clockwise sense void generbot(struct edge *p,int g){ //recursion for the bottom of

//the circle, this is entirely //analogous to genertop, and struct edge twofer[2],*ped; //only the differences will //be noted here struct complex temp, makecpx(int,int); double ct,st; void generbot(struct edge *,int); void bgener(struct edge *,struct edge *); int subtend(struct edge *),pp,qq; bgener(p,twofer); if(subtend(twofer)) ( ped=twofer; pp=-ped->b; //differs from genertop in that qq=ped->a; //b,a replaces d,c grafloutcount] .p=pp; graffoutcount] .q=qq; temp=makecpx(pp,qq); ct=temρ.x; st=temp.y; graffoutcount].x=-ped->alp.x*ct-ped->bet.x*st //differs from

+ped->gam.x+(ped->e.x)*4./(pp*pp+qq*qq); // genertop in that graffoutcount].y=-ped->aIp.y*ct-ped->bet.y*st //alp and bet are

+ped->gam.y+(ped->e.y)*4./(pp*pp+qq*qq); // multiplied by -1 outcount++; ped++; pp=-ped->b; //differs from genertop in qq=ped->a; //that b,a replaces d,c graf[outcount].p=pp; graffoutcount] .q=qq; temp=makecpx(pp,qq); ct=temp.x; st=temp.y; graf[outcount].x=-ped->alp.x*ct-ped->bet.x*st //differs from genertop

+ped->gam.x; //in that alp and bet are graf[outcount].y=-ped->alp.y*ct-ped->bet.y*st //multiplied by -1

+ped->gam.y; outcount++; return; } generbot(twofer,g+- 1 ); generbot(t wof er+ 1 ,g+ 1 ) ;

struct complex muit(struct complex u.struct complex v){ //multiplication struct complex temp; //of complex numbers

return(temp);

struct complex makecpx(int p,int q){ //input p,q returns struct complex temp; //the complex number double den; //(p-iq)/(p+iq) den=p*p+q*q; temp.x=(p*p-q*q)/den; temp.y=-2*p*q/den; return(temp);

int subtend(struct edge *p)( //subtend returns 1 if each of int nn; //of the two edges in the //array p subtends an angle nn=(p->a)*(p->c)-ι-(p->b)*(p->d); //approximately less than if(nn*nn<SCAT) //SCALE and returns 0 otherwise return(O);

P++; nn=(p->a)*(p->c)-t-(p->b)*(p->d); if(nn*nn<SCAT) retum(O); return(l);

void normalin(void){

//this routine calculates the modified struct complex cz[4],co[4],cm[4]; //0,+l,-l Fourier coeffs czero,cone,cmon //from the input Fourier coeffs struct complex temp; int n; czf0].x=-l //intialize array for czero calculation cz[0].y=0. cz[l].x=0. cz[l].y=0. cz[2].x=-l cz[2].y=0. cz[3].x=0. cz[3].y=0. co[0].x=0.; //initialize array for cone calculation co[0].y=0.; co[l].x=-l. co[l].y=0.; co[2].x=0.; co[2].y=l .; co[3].x=0.; co[3].y=0.; cm[0].x=0.; //initialize array for cmon calculation cm[0].y=0.; cm[l].x=0.; cm[l].y=0.; cm[2].x=0.; cm[2].y=-l.; cm[3].x=-l.; cm[3].y=0.;

cone=fourier[l]; cmon=fourier[-l];

temp=mult(fourier[n],cz[n%4]); czer.x+=-temp.x; //contribution to czer from the czer.y+=-temp.y; , //positive Fourier coeffs temp=mult(fourier[-n],cz[(4*N-n)%4]);

czer.x+=-temp.x; //contribution to czer from the czer.y+=-temp.y; //negative Fourier coeffs temp=mult(fourier[n],co[n%4]); cone.x+=-temp.x; //contribution to cone from the cone.y4-=-temp.y; //positive Fourier coeffs temp=mult(fourier[-n],co[(4*N-n)%4]); cone.x+=-temp.x; //contribution to cone from the cone.y+=-temp.y; //negative Fourier coeffs temp=mult(fourier[n],cmfn%4]); cmoπ.x+=-temp.x; //contribution to cmon from the cmon.y+=-temp.y; //positive Fourier coeffs temp=mult(fourier[-n],cm[(4*N-n)%4]); cmon.x+=-temp.x; //contribution to cmon from the cmon.y+=-temp.y; //negative Fourier coeffs

void normalout(void)f //this routine alters the int n,m; //output values in the struct complex temp,templ,makecpx(int,int); //array graf by adding double ct,st; //czero+cone*z+cmon*z^Λ(-l), //where z=(p-iq)/(p-> for (n=0;n<outcount;n++){ temp=makecpx(graf[n].p,graf[n].q); templ=mult(temp,cone); temp.y=-temp.y; temρ=mult(temp,cmon); graf[n].x+=czer.x+templ.x+temp.x; graf[n].y+=czer.y+templ.y+temp.y;

struct complex cogi(struct edge *p){ //cogi(p) calculates the wavelet coeff for the

//arrow underlying edge p from the Fourier coeffs int a,b,c,d,n; struct complex makecpx(int.int), mult(struct complex, struct complex); struct complex temp2,templ,temp,crun,cbl,cml,cb2,cm2; a=p->a; b=p->b; c=p->c; d=p->d; crun.x=0.; crun.y=0.; cm 1 =makecpx(b,-a); cm2=makecpx(d,-c); cbl=cml; cb2=cm2; for(n=2;n<=N;n++){ cbl=mult(cbl,cml); //updated nth power of cml cb2=mult(cb2,cm2); //updated nth power of cm2 if(n<MINMODE)//skip contribution to wavelet coeffs continue; //from low-frequency Fourier coeffs temp.x=n; temρ.y=b*d+a*c; templ=mult(cb 1 ,temp); temp.y=-temp.y; temp2=mult(cb2,temp); temp.x=temp 1.x+temp2.x; temp.y=templ.y+temp2.y; templ=mult(fourier[n],temp); crun.x+=templ.x; //contribution to wavelet coeff from crun.y+=templ.y; //nth Fourier coeff for n positive temp.x=-temp.x; temp2=mult(fourier[-n],temρ); crun.x+=temp2.x; //contribution to wavelet coeff from crun.y+=temp2.y; //nth Fourier coeff for n negative

} temp=crun; crun.x=-temp.y/4.; //overall multiple of i/4 crun.y=temp.x/4.; return(crun);

void tgener(struct edge *pc,struct edge *pn){ //tgener(p,twofer) stuffs twofer

//with descendants of p for top //of circle int a,b,c,d,bpd,apc; //matrix entries a,b,c,d, //bpd=b+d and apc=a+c struct complex e,alp,bet,gam; //wavelet coefficient e, trig coefficients alp.bet.gam struct complex axl,ax2; //index 1,2 refers to first.second struct complex bxl,bx2; //(counter-clockwise) descendants struct complex cxl,cx2; //prefix a,b,c refers to alp,bet,gam int ayl,ay2,byl,by2,cyl,cy2; //update procedure below defines x,y struct complex eu,et; //wavelet coeffs of first, second //descendant edges struct complex temp; struct complex cogi(struct edge *),mult(struct complex.struct complex); struct complex makecpx(int.int); double ct,st; struct edge *pnp; pnp=pn; //initialize pointers pnp++; a=pc->a; //initialize matrix entries b=pc->b; c=pc->c; d=pc->d; apc=a+c; bpd=b+d; pn->a=a; //stuff matrix of first descendant pn->b=b; pn->c=apc; pn->d=bpd; pnp->a=apc; //stuff matrix of descendant pnp->b=bpd; pnp->c=c; pnp->d=d; alp=pc->alp; //initialize lagged trig coeffs bet=pc->bet; gam=pc->gam; e=pc->e; //initialize wavelet coeff if(a>c){ //prepare for stuffing alp,bet,gam in case a>c ax 1.x=alp.x+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.x; ax 1.y=alp.y+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.y ;

ayl=2*(b*b-a*a); bxl.x=bet.x+(c*d-a*d-b*c-a*b)*4*e.x; bx 1.y=bet.y+(c*d-a*d-b*c-a*b)*4*e.y ; byl=4*a*b; ex 1.x=gam.x+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.x; ex 1.y=gam.y+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.y; cyl=-2*(a*a+b*b); ax2.x=alp.x+4*(d*d-c*c)*e.x; ax2.y=alp.y+4*(d*d-c*c)*e.y; ay2=2*(c*c-d*d); bx2.x=bet.x+8*e.x*c*d; bx2.y=bet.y+8*e.y*c*d; by2=-4*c*d; cx2.x=gam.x-4*e.x*(c*c+d*d); cx2.y=gam.y-4*e.y*(c*c+d*d); cy2-=2*(c*c-κl*d);

elsef //prepare for stuffing alp,bet,gam in case o=a ax 1.x-=alρ.x+4*e.x*(a*a-b*b); ax 1.y=alρ.y+4*e.y*(a*a-b*b); ayl=2*(b*b-a*a); bx 1.x=bet.x-8*e.x*a*b; bx I .y=bet.y-8*e.y *a*b; byl=4*a*b; ex 1.x=gam.x+4*e.x*(a*a+b*b); ex 1.y=gam.y+4*e.y *(a*a+b*b); cyl=-2*(a*a+b*b); ax2.x=alp.x+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.x; ax2.y=alι- _r(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.y; ay2=2*(c*c-d*d); bx2.x=bet.x+(a*d-a*b+b*c+c*d)*4*e.x; bx2.y=bet.y+(a*d-a*b+b*c+c*d)*4*e.y; by2=-4*c*d;

cx2.x=gam.x+(a*a+b*b •2*(a*c+b*d)-c*c-d*d)*2*e.x; cx2.y=gam.y+(a*a-ι-b*b ■2*(a*c+b*d)-c*c-d*d)*2*e.y; cy2=2*(c*c+d*d); } pn->e=eu=cogi(pn); //calculate wavelet coeffs pnp->e=et=cogi(pnp); //of the descendants pn->alp.x=ax 1.x+ay 1 *et.x; //update first alp pn->alp.y=axl.y+ayl *et.y; pn->bet.x=bx 1.x+by 1 *et.x; //update first bet pn->bet.y=bx 1.y+by 1 *et.y ; pn->gam.x=cx 1.x4-cy 1 *et.x; //update first gam ρn->gam.y=cx 1.y+cy 1 *et.y ; pnp->alp.x=ax2.x+ay2*eu.x; //update second alp pnp->alp.y=ax2.y+ay2*eu.y; pnp->bet.x=bx2.x+by2*eu.x; //update second bet pnp->bet.y=bx2.y+by2*eu.y; pnp->gam.x=cx2.x+cy2*eu.x; //update second gam pnp->gam.y=cx2.y+cy2*eu.y;

void bgener(struct edge *pc,struct edge *pn)( //bgener(p,twofer) stuffs twofer

//with descendants of p for //bottom of circle, int a,b,c,d,bpd,apc; //it is entirely analogous to struct complex e,alp,bet,gam; //tgener, and only the //differences will be noted here struct complex axl,ax2; struct complex bxl,bx2; struct complex cxl,cx2; int ayl,ay2,byl,by2,cyl,cy2; struct complex eu,et; struct complex ten._r , struct complex cogi(struct edge *),mult(struct complex.struct complex); struct complex makecpx(int.int); struct edge *pnp; double ct.st; pnp=pn; pnp++; a=pc->d; //differs from tgener in that d,-c,-b,a b=-(pc->c); //replaces a.b.c.d c=-(pc->b); d-=pc->a; apc=a+c; bpd=bκl; pn->a=bρd; //differs from tgener in that d,-c,-b,a pn->b=-apc; //replaces a,b,c,d pn->c=-b; pn->d=a; pnp->a=d; //differs from tgener in that d,-c,-b,a pnp->b=-c; //replaces a,b,c,d pnp->c=-bpd; pnp->d=apc; alp=pc->alp; bet=pc->bet; gam=pc->gam; e=pc->e; if(a>c)[ axl.x=alp.x+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.x; ax 1.y=alp.y+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e.y; ayl=2*(b*b-a*a); bxl.x=bet.x+(c*d-a*d-b*c-a*b)*4*e.x; bxl.y=bet.y+(c*d-a*d-b*c-a*b)*4*e.y; byl=4*a*b; cxl.x=gam.x+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.x; cxl.y=gam.y+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e.y; cyl=-2*(a*a+b*b); ax2.x=alp.x+4*(d*d-c*c)*e.x; ax2.y=alp.y+4*(d*d-c*c)*e.y; ay2=2*(c*c-d*d); bx2.x=bet.x4-8*e.x*c*d; bx2.y=be- ,-8*e.y*c*d; by2=-4*c*d; cx2.x=gam.x-4*e.x*(c*c+d*d); cx2.y=gam.y-4*e.y*(c*c+d*d); cy2=2*(c*c+d*d); } elsef ax 1.x=alp.x+4*e.x*(a*a-b*b); ax 1.y=alp.y+4*e.y*(a*a-b*b); ayl=2*(b*b-a*a); bxl.x=bet.x-8*e.x*a*b; bx 1.y=bet.y-8*e.y*a*b; byl=4*a*b; ex 1.x=gam.x+4*e.x*(a*a+b*b); cxl .y=gam.y+4*e.y*(a*a+b*b); cyl=-2*(a*a+-b*b); ax2.x=alp.x+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.x; ax2.y=alp.y+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e.y; ay2=2*(c*c-d*d); bx2.x=bet.x+(a*d-a*b+b*c+c*d)*4*e.x; bx2.y=bet.y+(a*d-a*b+b*c+c*d)*4*e.y; by2=-4*c*d; cx2.x=gam.x+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.x; cx2.y=gam.y+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e.y; cy2=2*(c*c+d*d);

pn->e=eu=cogi(pn); pnp->e=et=cogi(pnp); pn->alp.x=ax 1.x+ay 1 *et.x; pn->alp.y=ax 1.y+ay 1 *et.y ; pn->bet.x=bx 1.x+by 1 *et.x; pn->bet.y=bx 1.y+by 1 *et.y ; pn->gam.x=cx 1.x+cy 1 *et.x; pn->gam.y=cx 1.y+cy 1 *et.y ; pnp->alp.x=ax2.x+ay2*eu.x; pnp->alp.y=a.\2.y+ay2*eu.y; pnp->bet.x=bx2.x+by2*eu.x; pnp->bet.y=bx2.y+- _*eu.y; pnp->gam.x=cx2.x+cy2*eu.x; pnp->gam.y=cx2.y+cy2*eu.y;

//Here are the subroutines replacing tgener and bgener above in an //implementation which takes advantage of renormalization

void tgener(struct edge *pc,struct edge *pn)

int a,b,c,d,bpd,apc; struct complex e,alp,bet,gam; extern int outcount; extern struct pqxy graffCHUNK]; struct complex axl,ax2; struct complex bxl,bx2; struct complex cxl,cx2; struct complex eu.et; struct complex cogi(struct edge *); double ct.st; struct edge *pnp; pnp=pn; pnp++; a=pc->a; b=pc->b; c=pc->c; d=pc->d; apc=a+c; bpd=b+d; pn->a=a; pn->b=b; pn->c=apc; pn->d=bpd; pnp->a=apc; pnp->b=bpd; pnp->c=c; pnp->d=d; alp=pc->alp; bet=pc->bet; gam=pc->gam; e=pc->e; if(a>c){

bxl .y=bet.y-4.*e.y; cxl.x=gam.x; cxl.y=gam.y; ax2.x=axl.x; ax2.y=axl.y; bx2.x=bet.x; bx2.y=bet.y; cx2.x=gam.x-4.*e.x; cx2.y=gam.y-4.*e.y;

else{ axl.x=alp.x+4.*e.x; axl.y=alρ.y+4.*e.y; bxl.x=bet.x; bxl.y=bet.y; cxl.x=gam.x+4.*e.x; ex 1.y=gam.y4-4.*e.y ; ax2.x=axl.x; ax2.y=axl.y; bx2.x=bet.x+4.*e.x; bx2.y=bet.y+4.*e.y; cx2.x=gam.x; cx2.y=gam.y;

pn->e=eu=cogi(pn); pnp->e=et=cogi(pnp); alp.x=axl.x-2.*et.x; bet.x=bxl.x; gam.x=cx 1.x-2.*et.x; alp.y=axl.y-2.*et.y; bet.y=bxl.y; gam.y=cxl.y-2.*et.y; pn->alp.x=(2.*bet.x+alp.x+gam.x)/2.; pn->bet.x=(bet.x-alp.x+gam.x); pn->gam.x=(2.*bet.x-alp.x+3.*gam.x)/2.; pn->alp.y=(2.*bet.y+alp.y+gam.y)/2.; pn->bet.y=(bet.y-alp.y+gam.y); pn->gam.y=(2.*bet.y-alp.y+3.*gam.y)/2.; alp.x=ax2.x-2.*eu.x; bet.x=bx2.x; gam.x=cx2.x+2.*eu.x; alp.y=ax2.y-2.*eu.y; bet.y=bx2.y; gam.y=cx2.y+2.*eu.y; pnp->alp.x=(-2.*bet.x+alp.x-gam.x)/2.; pnp->bet.x=alp.x+bet.x+gam.x; pnp->gam.x=(2.*bet.x+alρ.x+3.*gam.x)/2.: pnp->alp.y=(-2.*bet.y+alp.y-gam.y)/2.; pnp->bet.y=alp.y+bet.y+gam.y; pnp->gam.y=(2.*bet.y+alρ.y+3.*gam.y)/2.;

double deriv(int a.int b.int c.int d.int p.int q){ double x,y,nep; x=p*p-q*q; y=-2*p*q; x=x/(p*p+q*q); y=y/(p*p+q*q); nep=-(a*a+b*b+c*c+d*d)+x*(a*a+b*b-c*c-d*d)-y*2.*(a*c+b*d); return(-2./nep);

void bgener(struct edge *pc,struct edge *pn)

int a,b,c,d,bpd,apc; struct complex e,alp,bet,gam; extern int outcount; extern struct pqxy graffCHUNK]; struct complex axl,ax2; struct complex bxl,bx2; struct complex cxl,cx2; struct complex eu.et; struct complex cogi(struct edge *); double ct.st; struct edge *pnp; pnp=pn; pnp++;

a=pc->d; b=-(pc->c); c=-(pc->b); d=pc->a; apc=a+c; bpd=b+d; pn->a=bpd; pn->b=-apc; pn->c=-b; pn-xl=a; pnp->a=d; pnp->b=-c; pnp->c=-bpd; pnp->d=apc; alp=pc->alp; bet=pc->bet; gam=pc->gam; e=pc->e;

if(a>c){ axl.x=alp.x+4.*e.x; axl.y=alp.y+4.*e.y; bxl.x=bet.x-4.*e.x; bxl.y=bet.y-4.*e.y; cxl.x=gam.x; cxl.y=gam.y; ax2.x=axl.x; ax2.y=axl.y; bx2.x=bet.x; bx2.y=bet.y; cx2.x=gam.x-4.*e.x; cx2.y=gam.y-4.*e.y;

elsef ax 1.x=alp.x+4.*e.x; ax 1 ,y=alp.y+4. *e.y ; bxl.x=bet.x; bxl.y=bet.y; ex 1.x=gam.x+4.*e.x; cxl.y=gam.y+4.*e.y; ax2.x=axl.x; ax2.y=axl.y; bx2.x=bet.x+4.*e.x; bx2.y=bet.y+4.*e.y; cx2.x=gam.x; cx2.y=gam.y;

} pn->e=eu=cogi(pn); pnp->e=et=cogi(pnp); pn->e=eu; pnp->e=et; alp.x=axl.x-2.*et.x; bet.x=bxl .x; gam.x=cx 1.x-2.*et.x; alp.y=axl.y-2.*et.y; bet.y=bxl.y; gam.y=cxl .y-2.*t , , pn->alp.x=(2.*bet.x+alp.x+gam.x)/2.; ρn->bet.x=(bet.x-alp.x+gam.x); pn->gam.x=(2.*bet.x-alρ.x+3.*gam.x)/2.; ρn->alp.y=(2.*bet.y+alp.y+gam.y)/2.; pn->bet.y=(bet.y-alp.y+gam.y); pn->gam.y=(2.*bet.y-alp.y+3.*gam.y)/2.;

alp.x=ax2.x-2.*eu.x; bet.x=bx2.x; gam.x=cx2.x+2.*eu.x; alp.y=ax2.y-2.*eu.y; bet.y=bx2.y; gam.y=cx2.y+2.*eu.y; pnp->alp.x=(-2.*bet.x+alp.x-gam.x)/2.; pnp->bet.x=alp.x+bet.x+gam.x; pnp->gam.x=(2.*bet.x+alp.x+3.*gam.x)/2.; pnp->alp.y=(-2.*bet.y+aIp.y-gam.y)/2.; pnp->bet.y=alp.y+bet.y+gam.y; pnp->gam.y=(2.*bet.y+alp.y+3.*gam.y)/2.;

Claims

Methods of D al Filtering and Multi-Dimensional Datε ompression

Using the Farey Quadrature and Arithmetic. Fan. and Modular Wavelets

Claims

I claim:

1. A method of compressing and reconstructing input data comprising the steps of: providing input data in an initial form; sampling the input data at the points of the Farey quadrature; transforming the Farey-sampled data into an intermediate form; quantizing the Farey-sampled data while in the intermediate form; using the quantized intermediate form to store or transmit the Farey-sampled data; de-quantizing the quantized intermediate form of the Farey-sampled data prior to reconstructing the data; and reconstructing the intermediate form of the Farey-sampled data in order to obtain reconstructed data for use.

2. A method as recited in claim 1 wherein: the step of transforming the Farey-sampled data into an intermediate form is achieved via a wavelet transformation; and the intermediate form of the Farey-sampled data takes the form of wavelet coefficients.

3. A method as recited in claim 2 wherein the wavelet transformation is based on arithmetic wavelets.

4. A method as recited in claim 3 wherein the arithmetic wavelets are defined in accordance with the following expressions:

+2 cos θ + 2 sin θ - 2 if 0 < 6> < f ,

+2 cos θ - 2 sin θ + 2 if f < θ ≤ π,

MW -2 cos θ - 2 sin θ - 2. ii π ≤ θ ≤ ,

-2 cos θ + 2 sin θ + 2 if ^~ < θ < 2 r; a b 1 0 and if A = ≠ / labels an arithmetic arrow; c dj ^ 0 1 and if A labels a top arithmetic arrow and a > c, then

( +4(d² - c²)cos θ + 8cd sin θ - 4(d² + c²); i :rt + ta- -n-132 l^d < <- V a ^ < + ta - -n-1 ₍₆2₊(b_d+₎ic-_)(₍o_α+₊c_c)₎ϋ ,

+2(α² + d² - b² - c² + 2ac - 2bd) cos θ +4(cd — ab — be — ad) sin θ + 2(α² + b² - c² - d² + 2ac + 2bd); if tnn -¹ 2(b+d)(α+c) _< β _{< tar}-l 26α A(Θ) ^{11 tan} (fc-f--^-(α-fc)² - - b^a* ^>

+2(b² + d² - a² - c² + 2αc - 2bd) cos θ +4(ab + cd — ad — be) sin θ + 2(-α² - b² - c² - d? + 2ac + 2bd);

" ;f t ^tainⁿ-¹ -^2h°^ < ft 0 < < t taorn. -* _(d 2_(d_b-₎i6_)(_{e_c-_α_a)₎2 , . 0; otherwise; and if A labels a top ; hmetic arrow and c> a, then

( +2(a² + c²-b²-d² - 2ac + 2bd) cos θ +4(ad + bc — ab — cd) sin θ +2(a² + b² + c² + d²- 2ac - 2bd); if _f,nn-l 2(-.-b)(c-«) < a < _tflT1-l l d ^{11 tan} (d-b) -{c-a) - ^υ ^ ^tan d^- -

+2(α² + d² - b² - c² + 2bd - 2ac) cos 0 +4(ad + bc + cd — ab) sin 0 ϋ_A(θ) = { + 2(α² + b²-c² -d²- 2bd - 2ac); ii tan ^.---- - ^ c ^ tan _(6+d)i---₍-._+c)2 ,

+4(α² - b²) cos θ - 8ab sin t9 + 4(α² + b²);

¹ i¹f t tnann-¹ ₍₆2₊(-_d.+₎i-._)(j_;α_α+₊e_c)₎j i -.f «' <-- tan -1₆s 2_6α_αs,

10; otherwise; and if J labels a bottom arithmetic arrow and a > — e, then

+2( ² + c² - ό² - d² + 2ac - 2bd) cos θ

+4(—ad — bc — ab cd) sin θ +2(a² + b² + c² + d² + 2ac + 2bd); if + _T1-l 2(6+--)(o+e) _{< ft <} +_ατ,-l 26α

" ^tan (6+<-)2-( +c)2 y tan _Γ-_?

+2(b² + c²-a² - d² + 2ac - 2bd) cos θ +4(ab - ad — be — cd) sin θ

#A(Θ) + 2(c² + d²- a² -b² + 2ac + 2bd); i "f ttnann-^{1 2}^^baμ.. _ <. θ v < t tann-¹ _(d 2_(c₆.-₎S6_)(j_:c_c-__α)₎-. ,

+4 (c² - d²) cos £ -

I 0; otherwise; and if A labels a bottom arithmetic arrow and — c > a, then

( +4(b² - a²) cos θ + 8ab sin θ - 4(a² + b²); if tan^"1 ^ ≤ θ ≤ tan^"1 _{ %_ ,

+2(b² + c² -a² -d² + 2bd - 2ac) cos θ +4(ad + bc + ab — cd) sin θ + 2(c² + d² -a² -b²- 2ac - 2bd);

#A Θ) if an-¹ _{^_ Z^ ≤ θ < an-¹^,

+2(6² + d² - a² - c² + 2bd - 2αc) cos θ +4(ad + bc + cd + ab) sin θ + 2(-α² - b² - c² - d² - 2ac - 26--); if tan^'1 A < θ < tan- Α^ ¹ , .

10; otherwise.

5. A method as reci in claim 2 wherein the wavelet transfc ation is based on modular wavelets.

6. A method as recited in claim 5 wherein the modular wavelets are defined in accordance with the following expressions:

where

{ A(Θ)\ if A ≠ U- T- A(Θ) = < ϋu-ι(θ)-28i_θ; XA = U-¹,

[tf_τ-ι(0)-f-2sint?; iϊA = T-\

0 -1 1 0 1 1 with S = U = and T ; and

1 0 1 1 0 1

+2 cos θ + 2 sin θ - 2 if 0<0< f,

+2 cos θ - 2 sin θ + 2 if f < θ < TV, (θ) -2 cos θ - 2 sin θ - 2 if π<θ < ^~f,

-2 cos θ + 2 sin θ + 2 if ^ < θ < 2π;

and if J = I ) T n _i ) ⁼ ^ labels an arithmetric arrow; and if A labels a top arithmetic arrow and a > c, then

' +4(d² - c²)cos θ + 8cd sin θ - 4(d² + c²); i ifi tnann-l < Λ

v < _. t tnan_ti"¹ ₍₆2₊(h_d+₎-!c._)(₍α_α+₊c_c)₎2,

+2(α² + d²-b²-c² + 2ac - 2bd) cos θ +4(cd — ab — bc — ad) sin θ + 2(a² + b²-c²-d² + 2ac + 26--);

;f _t-n-l 2(6+-.)(α+e) ^ n ^ _t»---l 2bα 0,4(0) = " ^tan (b+d)i-(α+c)2 < 0 < tan JJ 3T,

+2(6² + ci² - α² +4(α64- cd — + 2(-a² -b²

otherwise;

and if A labels a top arithmetic arrow and c> α, then { +2(o? + c² -b² -d² - 2αc + 2bd) cos θ +4(ad + bc — ab — cd) sin θ +2(a² + b² + c² + d²- 2ac - 2bd); i iif taηn-¹- _(d-S_₆i₎2_H₍-_c-_--_α-L₎2 S < ft t < tnarn.-¹ ^j ^2cd^-.r,

+2(α² + d² - 6² - c² + 26ci - 2αc) cos θ +4(ad + bc + cd- ab) sin θ

*A(Θ) = I + 2(a² + b²-c² -d²- 2bd - 2ac); iftaπ-ⁱ^≤fl≤tan-^^^,

+4(α² - 6²) cos 0 - 8ab sin 0 + 4( ² + 6²); n if t tnann"¹ ₍₆2₊(-_d.+₎2--_)(₍α_Q+₊c_c)_{)2 < β} f

I 0; otherwise; and if A labels a bottom arithmetic arrow and α > — c, then

( +2(a² + c²-b²-d? + 2ac - 2bd) cos 0 +4(—ad — be — ab — cd) sin 0 +2(α² + 6² + c² + d² + 2ac + 2bd); i ifi t tnann"¹ _(b ² ₊(^fe _d+₎2^rf_)(₍ ^α _α+₊ ^c _c)₎₂ < VQ < t ta_a„n-l

+2(6² + c²-a²-d² + 2ac - 2bd) cos 0 +4(α6 — ad — be — cd) sin 0

#A(Θ) + 2(c² + d² - a² - 6² + 2αc + 2bd); if tan-¹^ ≤ 0 ≤ tan-¹ > ,

+4 (c² - d²) cos 0 -

I 0; otherwise; and if J labels a bottom arithmetic arrow and — c > α, then

( +4(6² - a²) cos 0 + 8ab sin 0 - 4(α² + 6²); if tan-¹^, ≤ 0 ≤ tan-¹ %>_£_%_>

+2(6² + c² d² + 2bd - 2ac) cos 0

+4(ad + bc + ab- cd) sin 0

+ 2(c² + d² - a² - b² - 2ac - 26ci); if tnn-¹ 2(d-b)(c-a) _<

*Λ(Θ) = π tan _(d_₆₎..__(c__α)_ < ≤tan^"1^,

+2(6² + d² - a² - c²

+4(αc. + 6c + cd +

+ 2(-α² - 6² - c²

10; otherwise.

7. A method as recite .1 claim 2 wherein the wavelet transform. ->n is based on fan wavelets.

8. A method as recited in claim 7 wherein the fan wavelets are defined in accordance with the following expressions:

φ_A(θ) = ^~Tϋu _A(θ), n 0

where A(Θ)\ HA≠U-^λ,T-

#A(Θ) ϋu-ι(θ)-2s 0; HA = U^~ tf_τ-ι(0) + 2sin0; _Α = T^~1,

0 -1 with S and

1 0 u 1 0 1 1 and r = 1 1 0 1

+2 cos 0 + 2 sin 0 - 2; if 0<0< f,

+2 cos 0 - 2 sin 0 + 2; if f < 0 < π, θ) -2 cos 0 - 2 sin 0 - 2 if TΓ ≤Ø≤ ^~f,

-2 cos 0 + 2 sin 0 + 2 iϊ <θ≤2π;

and if A = / labels an arithmetric arrow;

and if A labels a top arithmetic arrow and a > c, then

( +4(d² - c²)cos 0 + 8cd sin 0 - 4(ci² + c²);

1 i1f t_fainn-¹ JT ^2c-^d^„ < _.0 v < _, t tnann-¹ ₍ ²(b b₊+_d+^d)(a+c)

(_h c.)γz_-₍(_n<H_-_rcγ)² ,

+2(α² + d? - b² - c² + 2ac - 2bd) cos 0

+4(cd — ab - 6c — ad) sin 0 + 2(α² + 6² -c²---² + 2αc + 26c.);

:r ₊__ — 1 2(6+t.)(o+e) -- Λ _<-- <. -1 2bα

0A(Θ) = { 11 tan _(6+d)11__(β+c)s < 0 < tan -p_^?_>

+2(6² + d² - α² - c² + 2αc - 2bd) cos +4(α6 + cd — ad — be) sin 0 + 2(-α² -b²-c² -d² + 2ac + 2bd); if tan^"1^ 6^— α-¹ < — 0 < — tan-¹ ( Jdⁱ—⁽l 67)_-'» —⁾⁽, (^cc^~—^o α⁾),_² ,'

10; otherwise;

and if A labels a top arithmetic arrow and c> a, then ( +2(a² + c² -b² -d² - 2ac + 2bd) cos 0 +4(ad + bc — ab — cd) sin θ +2(a² + b² + c² + d² - 2ac - 2bd); i ^- ! ^b-⁾{- ≤ θ ≤ tan-¹^

+2(α² + d² - b² - c² + 2bd - 2ac) cos 0

+4(ad + bc + cd — ab) sin 0

#A(Θ) = { + 2( ² + b²-c² -d²- 2bd - 2αc); i ifi ta_nn.-l ^^cd7_„ <- 0 v < _, t tnann-¹ ₍₆ ² ₊(^b _d+₎d_)(₍a_o+₊c_c)₎2 ,

+4(α² - 6²) cos 0 - 8α6 sin 0 + 4( ² + 6²);

1 if1 t t_na_nn-¹ ₍₆ ² ₊ ⁽⁶ _d+₎2^rf_⁾⁽ ₍ ^α _α+₊ ^c _c ⁾ ₎s S < ft P S < _t t_nan_n-l 5^ 2fea-?'

I 0; otherwise; and if A labels a bottom arithmetic arrow and α > — c, then

( +2(a² + c²-b²-d² + 2ac - 2bd) cos 0 +4(—ad — be — ab — cd) sin 0 +2(α² + b² + c² + d² + 2ac + 2bd); if t -¹ 2(t-+t-)(-.+c) < Λ < _{t n}-l 2bo

^{11 tan} ₍(,_+d) __(α+c)2 v S tan 52^7 _:

+2(6² + c² -a² -d² + 2ac - 2bd) cos 0 +4(α6 — ad — be — cd) sin 0

^■&A(Θ) = { + 2(c² + d² - a² - b² + 2ac + 2bd); if tan-¹ ^ ≤ 0 ≤ tan-¹ JH^_^l^_^~)_a ,

+4 (c² - ci²) cos 0 -

.0; otherwise; and if A labels a bottom arithmetic arrow and — c > α, then

J. A method as recitec claim 2 wherein the steps of sampling tl nput data at the points of the Farey quadrature and transforming the Farey-sampled data into wavelet coefficients is characterized by a binary cascade of arrow structures arising from the calculation of wavelet coefficients by sampling the input data in its initial form at the points of the Farey quadrature as defined along a circle in the complex plane.

10. A method as recited in claim 2 wherein the wavelet transformation is based on arithmetic wavelets and the step of transforming the Farey-sampled data into arithmetic wavelet coefficients includes the step of regularization.

11. A method as recited in claim 2 wherein the λvavelet transformation and the reconstruction are based on a finite number of wavelets.

12. A method as recited in claim 2 wherein: wavelet transformation is achieved via a plurality of wavelet transformation operations each for a specified arc along a circle in the complex plane; and the plurality of wavelet transformation operations are calculated in parallel with each other.

13. A method as recited in claim 1 wherein the quantized intermediate form of the data is stored and retrieved prior to de-quantizing and reconstructing the data for use.

14. A method as recited in claim 1 wherein: the steps of transforming the Farey- sampled data into the intermediate form and quantizing the intermediate form of the Farey-sampled data occur at a first physical location; the quantized intermediate form of the data is transmitted to a second physical location; and the steps of de-quantizing the quantized intermediate form of the Farey-sampled data and reconstructing the intermediate form of the Farey-sampled data in order to obtain reconstructed data occurs at the second physical location.

15. A method as recited in claim 14 further comprising the step of transmitting a reconstruction algorithm to the second physical location in order to accomplish reconstruction at the second physical location.

16. A method as recited in claim 1 wherein the input data is digital data.

17. A method as recited in claim 1 wherein the input data is analog data and the analog input data is sampled by placing analog sensors at locations corresponding to the points of the Farey quadrature.

18. A method as recited in claim 1 wherein the input data and the reconstructed data are multi-dimensional.

19. A method as recited in claim 2 wherein the step of reconstructing the intermediate form of the Farey-sampled data in order to obtain the reconstructed data is achieved via an inverse wavelet transformation which is characterized by a binary cascade.

20. A method as recited in claim 2 wherein the calculation of the wavelet transformation cludes the step of rei malization.

21. A method as recited in claim 1 wherein the sampling of data at the points of the Farey quadrature includes temporary storage of sequences of data in a buffer.

22. A digital signal processing method for obtaining data in a transform domain comprising the steps of: providing input data in an initial form; sampling the input data at the points of the Farey quadrature; using an intermediate transform to transform the Farey-sampled data into an intermediate form; and using a primary transform to convert the data in intermediate form of the Farey-sampled data into transform coefficients representative of the data in a transform domain.

23. A method as recited in claim 22 wherein the primary transformation consists of one of the groups of following transformations: Fourier, Hubert, Haar, Laplace, Bessel, Laguerre, Hermite, Chebyshev, Hotelling, Mersenne and Fermat.

24. A method as recited in claim 22 wherein the primary transformation consists of the Fourier transform, and the intermediate transform is a wavelet transformation based on arithmetic wavelets.

25. A method as recited in claim 24 wherein the intermediate transform is an arithmetic wavelet transform and the primary transform is a Fourier transform, and the conversion of the data from the intermediate form to the Fourier coefficients, c^A, is given by the following expressions:

a 6 1 0

A = ≠ = I labels an arithmetic arrow; and c d 0 1

(6 — d) — i(a — c)

- i (n³ — n) c^A — -[(c - a)² + (6 - d)²] (6 - d) + i(a - c) d — ic n ιa

+ 2(c² + d²) + 2(α + 6²) d + ic b + ia

and the coefficients c^A, c^A, c^A _x are given by ±ι = ^θh(c² - d² ± 2icd) + θ_r[(b - d)² (a - c)² + 2i(b - d)(a - c)}/2 + θ_t(a² - b² ± 2iab) + θ_e[(b + d)² (a + c)² + 2i(b + d)(a + c)]/2,

πc^A = 2θ_h(c² + d²) θ_r[(b - d)² + (a - c)²} + 2θ_t(a² + b²) θ_e[(b + d)² + (a + c)²}, where the angles are

-ι _/ 2c& x _Λ _ι _/ 2(6 + d)(a + c) -. ^tan ( ^>d_Λ ² T^)- fc = an ( ^■ ₍(b_h _ +_ d^)²_ _- _f(a„ _ +_ c)₂ ^< ) ,

26. A method as recited in claim 22 wherein the intermediate transform is a wavelet transform based on arithmetic wavelets, and the calculation of the intermediate transform includes the step of renormalization.

27. A method as recited in claim 22 wherein the primary transformation consists of the Fourier transform, and the intermediate transform is a wavelet transformation based on modular wavelets.

28. A method as recited in claim 22 wherein the intermediate transform is a wavelet transform based on modular wavelets, and the calculation of the intermediate transform includes the step of renormalization.

29. A method as recited in claim 22 wherein the primary transformation consists of the Fourier transform, and the intermediate transform is a wavelet transformation based on fan wavelets.

30. A method as recited in claim 22 wherein the intermediate transform is a wavelet transform based on fan wavelets, and the calculation of the intermediate transform includes the step of renormalization.

31. A method as recited in claim 22 wherein the sampling of data at the points of the Farey quadrature includes temporary storage of sequences of data in a buffer.

32. A method as recited in claim 22 wherein the input data is multi-dimensional.

33. A method as recited in claim 22 wherein the constants and coefficients required in the calculation of the intermediate transform and in the calculation of the primary transform are archived and recovered from memory at run time.

34. A method of data processing comprising the steps of: providing input data in an initial form; sampling the input data at the points of the Farey quadrature; transforming the Farey-sampled data to a second form; and processing the data in the second form.

35. A method as recited in claim 34 wherein: the step of transforming the Farey- sampled data into the second form is achieved via a wavelet transformation; and the second form of the Farey-sampled data takes the form of wavelet coefficients.

36. A method as recited in claim 35 wherein the wavelet transformation is based on arithmetic wavelets.

37. A method as recited in claim 36 wherein the arithmetic wavelets are defined in accordance with the following expressions: +2 cos θ - 2 sin 0 - 2 ifO≤ø≤ f,

+2 cos 0 - 2 sin 0 + 2 if f < 0 ≤ 7r, θ) = -2 cos 0 - 2 sin 0 — 2 if π<0< ^,

-2 cos 0 + 2 sin 0 + 2 if < 0 < 2π

and if -Λ --= ) ^ ( n _i ) ⁼ ^ ^^els an arithmetric arrow;

and if A labels a top arithmetic arrow and a > c, then

and if A labels a top arithmetic arrow and c> a, then

( +2(a² + c²-b²-d² - 2αc + 26ci) cos 0 +4(ad + bc — ab — cd) sin 0 +2(α² + b² + c² + d² - 2ac - 26d); if tan^"1 (c *--^idbrf) -⁾⁽;(c-^a 2cc! rαλ)i - <--.θ< tan^"1 d^

+2(α² + d² - 6² - c² + 2bd - 2αc) cos 0

+4(ad + bc + cd — ab) sin 0

^rø = + 2(α² + b² -c² -d² - 2bd - 2αc); i iif tnann-¹ ^ ²^^d„ <.0 c < ^ t tnann-¹ ₍₆ ² ₊ ^(b _d+₎2^d_⁾⁽ ₍ ^a _α+₊ ^c _c ⁾ ₎ „,

+4(α² - 6²) cos 0 -

.0; otherwise;

and if J labels a bottom arithmetic arrow and α > — c, then '

+2(6² + c² - a² - d² + 2ac - 2bd) cos 0 +4( 6 — ad — be — cd) sin 0

#A(Θ) = + 2(c² + d² - a² - 6² + 2αc + 2bd); i iif t tann_n-^{1 2}^^b y- < ft f < tann- _(d 2_(<₆--₎₂6_)(₍e_c-_α_α)₎2 , +4 (c² - d²) cos 0 - 8cci sin 0 + 4(c² + ci²); i :fi t ta--n--l _(d 2_(d₆-₎₄b_)(₍c_c-_α_α)₎2 <■ » a <■ tan -1 ^ 2ccj c.^-c²'

I 0; otherwise; and if A labels a bottom arithmetic arrow and —c> a, then

+4(6² a²) cos 0 +

+2(6² + c² - α² - d² + 26d - 2αc) cos 0 +4(0--" + 6c + ab — cd) sin 0

+ 2(c² + d² 2ac-2bd); if tan ._^" 1 2(d-b)(c- )

$A(Θ) = < (-.-b)²-(c-α)² <0<tan-^x d²²-^cdc²'

+2(6² + d² -a? -c² + 2bd - 2ac) cos 0 +4(αc- + 6c + cd + ab) sin 0 + 2(-α² - 6² - c² - ci² - 2αc - 26c.); i 1f1 <

< ft 0 < <

10; otherwise.

38. A method as recited in claim 35 wherein the wavelet transformation is based on modular wavelets.

39. A method as recited in claim 38 wherein the modular wavelets are defined in accordance with the following expressions:

ΦA(Θ) = ^~~~ n un_A(θ), n>ϋ

where

0 -1 1 0 1 1 with S — and T = ana

1 0 ' 1 1 0 1

+2 cos 0 + 2 sin 0 - 2 if O<0 — < £ 2'

+2 cos 0 - 2 sin 0 + 2 if f < θ) = -2 cos 0 - 2 sin 0 - 2 if π < 0 < 3τr 2

I -2 cos 0 + 2 sin 0 + 2 if r < θ < 2τr;

and ifJ= .

I.J I = / labels an arithmetric arrow;

and if A labels a top arithmetic arrow and a > c, then

( +4(d² - c)cos 0 + 8cci sin 0 4(ci²+c²); if tan-¹^ < 0 ≤ tan-¹ _{^_y ,

+2(a² + d²-b²-c² + 2ac - 26ci) cos 0 +4(cc- — ab — be — ci) sin 0 + 2(α² + 6² - c² - ci² + 2 c + 26i);

^ an^"1 ₍ ffl_tt+$ ≤ θ ≤ tan-¹^,

^•&A(Θ)

+2(6² + d² - a² +4(α6 + cd — + 2(-α² - 6²

otherwise;

and if J labels a top arithmetic arrow and c> a, then

( +2(a² + c²-b²-d²- 2ac + 2bd) cos 0 +4(ad + bc — ab — cd) sin 0 +2(α² + 6² + c² + d² - 2αc - 26c-);

^ ^-'₍S^ < θ < tan-¹^,

+2(α² + d -b²-c² + 2bd - 2ac) cos 0 +4(αc- + 6c + cd — ab) sin 0 A(Θ) = { + 2(α² + b²-c² -d²- 2bd - 2ac); if tan-¹ ^ ≤ 0 ≤ tan-¹ ₍₆ ²|g£^_a ,

+4( ² - 6²) cos 0 - 8α6 sin 0 + 4(α² + 6²); iiff taann^-"11 ₍ff_ζ ffffl_ttt < g ≤ tan^"1^

.0; otherwise;

and if A labels a bottom arithmetic arrow and α > — c, then ( +2(a² + c² - b² - d² + 2ac - 2bd) cos 0 +4(—ad — be — ab — cd) sin 0 +2(α² + b² + c² + d² + 2ac + 26--); i if nnn- l ₍₆2₄(-b_d+₎id_)(_(α+_-t-c_c)₎i _< ---- c a < ---:

+2(6² + c² - a² - d² + 2ac - 2bd) cos 0 +4(α6 — ad — be — cd) sin 0 A (Θ) = < + 2(c² + d² - a² - b² + 2ac + 2bd); i ift t tnann" ^{1 2}^^6α^„- < 0 v < t tnann^{" 1} _(d ²_(^rf ₆-₎2^fc_)(₍ ^c _c-_^Q _α)₎2 , +4 (c² - d²) cos 0 - 8cci sin 0 + 4(c² + d²); i :ri tan-l _(d 2_(-₆--₎j6_)(₍c_c-_α_α)_)i <- f) S<•

. 0; otherwise; and if J labels a bottom arithmetic arrow and — c > α, then

f +4(6² - a²) cos 0 + 8α6 sin 0 - 4(α² + 6²); i •i_f tann-l p 2ibα? S --- f o _< t _f n -1 _(d 2-(-.d₆-₎S6-l)(₍c_c-_α_α)'₎. ,

+2(6² + c² - a² - d² + 26ci - 2αc) cos 0 +4( c- + bc + ab — cd) sin 0 + 2(c² + d² - a? - b² - 2ac - 26c-);

:r _t --- -1 2(rf-b)(c-α) ^ n ^ + ---. -1 2ccj ϋ_A(θ) = { ^{ύ tan} (ci-b)²-(c-α)^ ≤ ^u ≤ tan 5—7 ,

+2(6² + d² - α² - c² + 26cZ - 2αc) cos θ +4(ad + bc + cd + ab) sin 0 + 2(-α² - 6² - c² - d² - 2αc - 26ci);

1 0; otherwise.

40. A method as recited in claim 35 wherein the wavelet transformation is based on fan wavelets.

41. A method as recited in claim 40 wherein the fan wavelets are defined in accordance with the following expressions:

where #Λ(0)

0 1 0 1 1 with S , and T and

1 1 1 0 1

+2 cos 0 + 2 sin 0 - 2 if0<0< f,

+2 cos 0 - 2 sin 0 + 2 if f < ≤ π,

^■& = -2 cos 0 - 2 sin 0 - 2 if π< < 3-1

-2 cos 0 + 2 sin 0 + 2 if 4f < 0 < 2π;

and if -A = f , ) = / labels an arithmetric arrow;

and if A labels a top arithmetic arrow and a> c, then

and if A labels a top arithmetic arrow and c> a, then

( +2(a² + c² b² -d² - 2ac + 26c.) ^{cos θ} +4(ad + bc — ab — cd) sin 0 +2(α² + 6² + c² + ci² - 2αc - 2bd) if tan^" 2(d-b)(c- 2cd

(<--b) -(c ^-α) - <-i 0 < tan^" -1 d J5^ — :

+2(α² + d²-b²-c² + 2bd - 2ac) cos 0 +4(αc- + 6c + cd — ab) sin 0 A(Θ) + 2(α² + 6² - c² - d² - 26ci - 2αc); if nn-l 2cd _{< β < t} -i 2(h+d)(o+c) ii tan 52 ^ -- " tan _(b+d)2__(α+c)j ,

+4(α² - 6²) cos 0 -

10; otherwise;

and if A labels a bottom arithmetic arrow and a > —c, then ( +2(a² + c² - b² - d² + 2ac - 2bd) cos 0 +4(—ad — be — ab — cd) sin 0 +2( ² + b² + c² + d² + 2ac + 2bd); i iif tnann"¹ ₍₆2₊(b_d+₎2d_)(₍α_α+₊c_c)₎i _< -. » β - _<. t , an - 1 ₆42_b_αs ,

+2(6² + c² - a² - ci² + 2αc - 26ci) cos 0 +4(α6 — ad — be — cd) sin 0

#A (Θ) = < + 2(c² + d² - a² - b² + 2ac + 2bd); if tan-¹ ^, ≤ 0 ≤ tan^"1 _{^__ Z ,

+4 (c² - d²) cos 0 - 8cd sin 0 + 4(c² + ci²); i ifi t tnann"¹ _d 2_(-₆.-₎ι6_)(c_c-_α_α)₎4 < -. ff a <-i _t ta„_nn-l 352 2cc<d

3- ,

I. 0; otherwise; and if A labels a bottom arithmetic arrow and — c > α, then

( +4(b² - a²) cos 0 + 8α6 sin 0 - 4( ² + 6²);

:f i._{n n} -l 2bα <- a <• t-.--. -l 2(d-b)(c-α) it tan ^J^^Γ σ -S tan _(d_₆₎j_ _c__o)4 ,

+2(6² + c² - α² - d² + 2bd - 2ac) cos 0 +4(α i + 6c + ab — cd) sin 0 + 2(c² + d² - a² - b² - 2ac - 2bd); δ_A ) if tan^"1 ff ff-g₈ ≤ 0 ≤ tan-¹ ^,

+2(6² + d² - α² - c² + 26d - 2 c) cos 0 +4(αc- + 6c + cd + ab) sin 0 + 2(-α² - 6² - c² - d² - 2αc - 26d); if tan-¹ ^ ≤ 0 ≤ tan-¹ ₍%£% , ; otherwise.

42. A method as recited in claim 35 wherein the steps of sampling the input data at the points of the Farey quadrature and transforming the Farey-sampled into wavelet coefficients is characterized by a binary cascade of arrow structures arising from the calculation of wavelet coefficients by sampling the input data in its initial form at the points of the Farey quadrature as defined along a circle in the complex plane.

43. A method as recited in claim 42 wherein the wavelet transformation is based on arithmetic wavelets and the step of transforming the Farey-sampled data into arithmetic wavelet coefficients includes the step of regularization.

44. A method as recited in claim 35 wherein the wavelet transformation is based on a finite number of wavelets.

45. A method as recited in claim 35 wherein: wavelet transformation is achieved via a plurality of wavelet transformation operations each for a specified arc along a circle in the complex plane; anc e plurality of wavelet transformation opei αns are calculated in parallel with each other.

46. A method as recited in claim 34 wherein the input data is digital data.

47. A method as recited in claim 34 wherein the input data is analog date and the analog input data is sampled by placing analog sensors at locations corresponding to the points of the Farey quadrature.

48. A method as recited in claim 34 wherein the input data is multi-dimensional.

49. A method as recited in claim 34 wherein the sampling of data at the points of the Farey quadrature includes temporary storage of sequences of data in a buffer.

50. A digital processing system comprising a digital filter that receives an input signal and outputs a transformed signal, the digital filter being an arithmetic wavelet filter in which the arithmetic wavelets are defined in accordance with the following expressions:

+2 cos 0 + 2 sin 0 - 2 if 0 < 0 < f ,

+2 cos 0 - 2 sin 0 + 2 if f < 0 < -, ϋι(θ) = -2 cos 0 - 2 sin 0 - 2 if π < 0 < ,

-2 cos 0 + 2 sin 0 + 2 if < θ < 2τr;

and if J = labels an arithmetric arrow;

and if A labels a top arithmetic arrow and a > c, then

( +4(d² - c²)cos 0 + 8cd sin 0 - 4(d² + c²); i :tf i t.a„„n- l <- «

0 a < t.an - 1 _(b 2₊(b_d+₎4c-_)(₍α_α+₊c_c)₎ ,

+2

ϋ_A{θ) = {

+2(6² + d²

+4(α6 + cd -

+ 2(-α² - b²

1 0; otherwise:

and if A labels a top arithmetic arrow and c > a, then { +2(a² + c² -b² -d² - 2ac + 2bd) cos 0 +4(αci + 6c — ab — cd) sin 0 +2( ² + 6² + c² + ci² - 2αc - 26ci); if n -l 2(d-b)(_C-a) _< „ _< , -l 2cd ii tan _(d_₆)4__(c__α)4 if i tan _d2__c2 ,

+2( ² + d² - b² - c² + 2bd - 2ac) cos 0 +4( ci + bc + cd — ab) sin 0

^■d_A(θ) = { + 2(α² + 6² - c² - d² - 26ci - 2αc); if tn -l 2cd <- a < to "^{1 2}(fe+rf)("+c) ⁱⁱ tan -rr^r » tan (_.+--)* -(α-t-c)⁴ >

+4(α² - 6²) cos 0 - 8 6 sin 0 + 4(α² + 6²); if tan-¹ _{b +f \_⁺ ₊ ≤ 0 ≤ tan-¹,^, l 0; otherwise; and if A labels a bottom arithmetic arrow and α > — c, then

f+2( ² +c² 6² - ci² + 2αc - 26ci) cos 0 +4(— ci — be — ab — cd) sin 0 +2( ² + 6² + c² + ci² + 2αc + 26ci); if tnn-l 2(b+d)(α+c) _< Λ <^• _fnr.-l 2b ii tan _(ό+d)2_(_Q+c)2 f tan j S^*r^

+2(6² + c² - α +4(α6 — ad ϋ_A(θ) = { + 2(c² + d²

+4 (c² - ci²) cos 0 - 8cci sin 0 + 4(c² + d²); if tan^"1 (ff g_* ≤ <? ≤ ten-¹^ r-r.

.0; otherwise; and if J labels a bottom arithmetic arrow and — c > a, then

( +4(b² - a²) cos 0 + 8 6 sin 0 - 4(α² + 6²); if tnn-l 2bα _{< ff < tnT1}-l 2(d~b)(c-α) ii tan -^ϊ f S tan _(d-₆₎ ²-_(c-_α) ² '

+2(6² + c² - α² - ci² + 26d - 2αc) cos 0 +4( ci + bc + ab — cd) sin 0 + 2(c² + d² -a² -b² - 2ac - 2bd); A(Θ) = { i tan^"1 _{!iy_ Z ≤ θ ≤ tan-¹^,

+2(6² + d² - α² - c² + 26d - 2αc) cos 0 +4(αci + 6c + cd + ab) sin 0 + 2(-α² - 6² - c² - d² - 2αc - 26ci); i ifi tnann-l _< _, v a < _, _f ta--n-l ₍₆2₊(b_d+₎ad_)(₍a_α+₊c_c)₎₂ ,

10; otherwise.

51. A digital processir. ystem comprising a digital filter that rec as an input signal and outputs a transformed signal, the digital filter being a modular wavelet filter in which the modular wavelets are defined in accordance with the following expressions:

ΦA(Θ) = ∑ n ursA(θ), n>0

where

( &); i(AφU- T- A(Θ) < ϋu-ι(θ)-2sinθ; HA = U-¹, l_{> r}-ι(0) + 2sin0; ifA = T^~1,

1 0 1 1 with 5 = U = andT: ; and

1 0 1 1 0 1

+2 cos 0 + 2 sin 0 - 2 if0<0< f,

+2 cos 0 - 2 sin 0 + 2 if f < 0 < π, θ) = -2 cos 0 - 2 sin 0 - 2 if π<0< ,

-2 cos 0 + 2 sin 0 + 2 if ≤ < 2τr;

and if J = { ) ^ ( n ι ) ⁼^ labels an arithmetric arrow;

and if A labels a top arithmetic arrow and a> c, then

and if A labels a top arithmetic arrow and c> a, then f +2(α + c² d² — 2αc + 26d) cos

+4( d + bc- ab- cd) sin 0

+2(a² + b² + c² + d² 2 c-26d); if tan-¹ * !_ -_ < θ ≤ tan-¹^ d -c²^

*A(Θ) =

+4(α² - 6²) cos 0 - 8α6 sin 0 + 4(α² + 6 ,2) if tm- W°⁺A> < θ < tan-¹ 2b

(b+d)²-{a+c)² — ^w — ^uσ-" W-

K.0; otherwise; and if A labels a bottom arithmetic arrow and a > — c, then

( +2(a² + c² -b²-d² + 2ac - 26d) cos 0 +4(— ad — be — ab — cd) sin 0 +2(α² + 6² + c² + d² + 2αc + 26d); i ifi t tna_nn"¹ ₍₆(l^b _d+₎S^d_H₍°_α+₊ ^c _c)_{) <} i-f Λ 1 _< t ^,—an -i ι-∑ 2_b_7,

+2(6² + c² - a² - d² + 2ac - 2bd) cos 0 +4(α6 — ad — bc — cd) sin 0 A(Θ) = < + 2(c² + d² - a² - b² + 2ac + 26d); i ifi tnann-l 2_bα_ι < ^ v? < i t ta_nn-l _(d 2_(d₆-₎b_)(₍c_c-_α_α)₎2 ,

+4 (c² - d²) cos 0 - 8cd sin 0 + 4(c² + d²); if an^"1 ffff^ < θ < tan^"1^, l 0; otherwise; and if A labels a bottom arithmetic arrow and — c > a, then

( +4(b² - a²) cos 0 + 86 sin 0 - 4( ² + 6²); if tan-¹^^ <0< tan^" -1 2(d-b)(c-α) ⁱ _d-b) -{c-ay>

+2(6² + c²-a²-d² + 2bd - 2ac) cos 0 +4( d + bc + ab — cd) sin 0 + 2(c² + d² - α² - 6² - 2αc - 26d); i ,„n-l ;j- 2

#A(Θ) = { iif t tannn-¹ _(d 2_(d₆-₎4b_)(₍e_c-__α)₎4 _< β V _< S _t ta_! —c ?,

+2(6² + d² -a² -c² + 2bd - 2ac) cos 0 +4( d + 6c + cd + ab) sin 0 + 2(-α² - 6² - c² - d² - 2αc - 26d); iftan^^≤ø≤tan-¹^^^, ; otherwise.

02. A digital processing /stem comprising a digital filter that recc :s an input signal and outputs a transformed signal, the digital filter being a fan wavelet filter in which the fan wavelets are defined in accordance with the following expressions:

where

HA≠U- T- A(Θ) = { ϋ_-ι (θ) - 2 sin 0; XA = U-¹, tf_τ-ι(0) + 2sin0; i£A = T-¹,

1 0 1 1 with S = U = and T ; and 1 1 0 1

+2 cos 0 + 2 sin 0 - 2 if O<0< ,

+2 cos 0 - 2 sin 0 + 2 if f < 0 < π,

-2 cos 0 - 2 sin 0 - 2 ifτr<0< ^~f,

-2 cos 0 + 2 sin 0 + 2 if ≤ θ < 2π;

and if J = ( j φ I J = / labels an arithmetric arrow;

and if A labels a top arithmetic arrow and a > c, then

' +4(d² - c²)cos 0 + 8cd sin 0 - 4(d² + c²); iftan^^≤ø≤tan-¹,^^^,

+2(α² + d² - 6² - c² + 2αc - 26d) cos 0 +4(cd — ab — bc — ad) sin 0 + 2(α² + 6² - c² - d² + 2αc + 26d); ι ₇?._Λ( _/0„-) = { " _jf t tannn-l ₍₆2₊(b_d+₎2d_)(₍α_α+₊c_c)₎2 _< f a _< < t .an

+2(6² + d²-a²-c² + 2ac - 26d) cos 0 +4(α6 + cd — ad — be) sin 0 + 2(- ² - 6² - c² - d² + 2 c + 26d); if tan-¹^ ≤ θ ≤ t_Αn- (d-^db^~)Α-(^cc-^aλa)² ; otherwise;

and if A labels a top arithmetic arrow and c> a, then +2(α² b2² - d² - 2ac + 2bd) cos θ

+4(αd + 6c - ab — cd) sin 0 +2(α² + 6² + c + d² 2αc-26d); i iif t tnann"¹ _(d ²_(^d ₆-₎2^fc--)-(₍c_c-_a_α)₎2 - <--- P a - _{< f} t_na_nn-l _d22_ed_c4,

+2(α² + d² - b² - c² + 2bd - 2ac) cos 0 +4(αd + bc + cd — ab) sin 0

^•9A(Θ) + 2( ² + 6² - c² - d² - 26d - 2αc); iftan-^≤ø≤tan-¹,^!^,

+4(α² - 6²) cos 0 - 8α6 sin 0 + 4(α² + 6²); if tnn-l 2(6+^t.)(α+c) _< a _{< trιrι}-l 2ba

^{11 tan} (b+d)²-( +c)² - ^υ ≤= tan ^- v.0; otherwise; and if A labels a bottom arithmetic arrow and a > — c, then

' +2(α² + c² - 6² - d² + 2αc - 26d) cos 0 +4(— ad — bc — ab — cd) sin 0 +2(α² + b² + c² + d² + 2ac + 26d); if tnn-l 2(b+d)(α+c) _{< β <} t _n-l 2b ^{11 tan} (b+d) -(_α+c)² - ° ^ ^tan ^ 7)

+2(b² + c² -a² - d² + 2αc - 26d) cos 0

+4( 6 — ad — bc- cd) sin 0

*A(Θ) = { + 2(c² + d²-a² -6² + 2αc + 26d); if tan^"1 ^ τ < - ^tCl11 (d-b)²-(c-θ)²

+4 (c² - d²) cos 0 - 8cd sin 0 + 4(c² + d²); iftan-¹^-^. ≤ø≤tan-¹ 2cd d^-c¹ l 0; otherwise; and if A labels a bottom arithmetic arrow and — c > α, then

' +4(6² - a²) cos 0 + 8α6 sin 0 - 4(α² + 6²); i ifi t tnann-l < -ϊ- β ^y S <^■. t ta_nnn- _(d 2

_(d₆-₎5b_)(₍-_c ^,-_a_α)₎s ,

+2

A(Θ) =

+2(6² + d² - a² - c² + 2bd - 2αc) cos 0 +4(αd + bc + cd + ab) sin 0 + 2(-α² - 6² - c² - d² - 2 if an-¹

otherwise.

53. A method of obtai -g an inverse transform of input digital a comprising the steps of: providing input digital data in an initial form; transforming the data in the initial form to an intermediate form; and transforming the intermediate form of the data using an intermediate inverse transform to obtain output values at the points of the Farey quadrature.

54. A method as recited in claim 53 wherein the input digital data is in the form of Fourier coefficients and the combination of transforming the data from the Fourier coefficients to the intermediate form and from the intermediate form to obtain output values at the points of the Farey quadrature constitutes an inverse Fourier transform.

55. A method as recited in claim 54 wherein the intermediate form of the data is in the form of wavelet coefficients.

56. A method as recited in claim 55 wherein the wavelet coefficients are coefficients for arithmetic wavelets,

with U -1 and conversion of the data from the initial form to the intermediate form is based on the following expressions:

Λnθ = cos nθ + i sin nθ

= - tø + c^~e» + c__ιe-^iθ] + ~_ά ∑{ ( + η") + ^ V +± ξ ( - _η ⁿ) }ϋ_A(θ), v - ζ where the sum is over all arithmetic arrows, the underlying chord of which has complex endpoints ξ, η, and

-1, n = 0(4) 0, n ≡ 0(4) 0, n ≡ 0(4)

0, n ≡ 1(4) -1, n ≡ l(4) 0, n ≡ 1(4)

-1, 71 -= 2(4) i, n ≡ 2(4) -1 ". _ i, n ≡ 2(4)

0, n ≡ 3(4) 0, n ≡ 3(4) -1, n = 3(4).

57. A method as recited in claim 55 wherein the wavelet coefficients are coefficients for modular wavelets, A , and conversion of the data from the initial form to the intermediate form is based on the following expressions:

e Λ"n"θ" — cos nθ + i sin nθ = - [c^ + c^e^iθ + cϋ₁

+ ϊ Ε{ⁿ^ + Η + ij-zf ( " - *?")} tΨuA ^θ) - 2φ_A(θ) + Φu-^s _A{θ)) , where the sum is over < -rithmetic arrows, the underlying chord c 'hich has complex

1 0 endpoints ξ, η, where U t±i - , and where ±1 1

-1, n = 0(4) 0, n ≡ 0(4) 0, n ≡ 0(4)

0, n ≡ 1(4) -1, 7i = l(4) 0, n ≡ 1(4)

^Co = -1, n = 2(4) i, n ≡ 2(4) -i, n ≡ 2(4)

0, n ≡ 3(4) 0, n ≡ 3(4) -1, n = 3(4).

58. A method as recited in claim 55 wherein the wavelet coefficients are coefficients for fan wavelets, φ_A, and conversion of the data from the initial form to the intermediate form is based on the following expressions: e^inθ = cos nθ + i sin nθ

= - [c% + c⁷{e^iθ + c_₁e-^iθ]

+ I EI Γ + ⁿ) + ~ (T - V^~)) [ΦA(Θ) - ΦUA(Θ)], η - ξ where the sum is over all arithmetic arrows, the underlying chord of which has complex endpoints ξ, η, where N = I I , and where

-1, 7i = 0(4); ( 0, n = 0(4); ( 0, n ≡ 0(4) . 0, n ≡ 1(4); _ -1, n = l(4); „ I 0, n = l(4)

^C° ^~ -1, n ≡ 2(4); ^{Cl ~} i, n = 2(4); -¹ ] -*, n = 2(4)

0, n = 3(4); 0, n ≡ 3(4); I -1, n ≡ 3(4).

59. A method as recited in claim 53 wherein the intermediate form of the data is in the form of wavelet coefficients.

60. A method as recited in claim 53 wherein the initial form of the input digital data is Fourier coefficients, the intermediate form of the data is in the form of wavelet coefficients, and the output values at the points of the Farey quadrature are obtained via the use of a binary cascade of arrow structures arising from the calculation of wavelet coefficients as defined along a circle in the complex plane.

61. A method as recited in claim 60 further comprising the step of storing terminal arrow structures in order for restart or iterative refinement.

62. A method as recited in claim 59 wherein the wavelet coefficients are coefficients for arithmetic wavelets.

63. A method as recited in claim 59 wherein the wavelet coefficients are coefficients for modular wavelets.

64. A method as recited in claim 59 wherein the wavelet coefficients are coefficients for fan wavelets.

65. A method as recited in claim 59 wherein the input data is multi-dimensional.