Neville Dias, Subramanyam Ganapathi, and Sudhir Ramanna
TechOnline

The multi-resolution decomposition of a function is similar to passing the signal through a bank of orthogonal filters. In signal-detection applications, a signal in the presence of noise is decomposed with the wavelet transform. A wavelet matched to the signal would produce a sharper peak in the time-scale space as compared to standard non-matched wavelets. The sharper peak is the result of the similarity between the shape of the wavelet and the signal to be detected. However, this requires the design of a suitable wavelet for a given application. The design techniques developed to date do not specifically address this need for maximizing correlation of the designed wavelet with a specific signal.

The technique of deriving orthonormal wavelets directly from a signal of interest has specific applications in signal detection, image enhancement, and target detection. This article describes a technique of designing wavelets that almost match a desired signal and also satisfy the necessary conditions for an MRA (multiresolution analysis). The matching algorithm is performed on the wavelet spectrum magnitude only and, hence, produces symmetric wavelets. This is sufficient for target detection.

The Radon transform involves taking projections of the object at various angles. For target detection, it reduces a 2D-detection problem into a set of 1D cases. We will use the projections of the training object for designing matched wavelets, corresponding to different angles. In the detection phase, the designed wavelets detect the rotated, shifted, and scaled version of the object.

<y_{j, k}, y_{l, m}> = d_{j, l} * d_{k, m}      (1)

Equation 1 indicates that the family of wavelets y is orthogonal in two dimensions. Integer translations of the wavelet at a given scale are orthogonal to one another (d_{j, l}), as are wavelets at two different scales (d_{k, m}). Orthogonality across scales indicates that W_j W_l for all j l, where W_j is a subspace of L²(R). In an orthonormal MRA, a signal f(x), is decomposed into an infinite series of detail functions, g_j(x), such that:

f(x) = S_{j = -} g_j(x)      (2)

The decomposition of the function involves the projection of the function onto two orthogonal subspaces, V and W. Further decomposition is carried on the subspace V. The detail captured by the wavelet, at a particular level is unique.

<f_{j, k}, y_{j, m}> = 0      (3)

For the scaling function f to satisfy the orthonormality condition

<f_{j, k}, f_{j, m}> = d_{k, m}      (4)

The Fourier Transform of the Equation 4 gives the Poisson's equation given by:

S_{m = -}⁺ |f(w+2pm)|² = 1      (5)

The scaling function feV_oV_-1 and wavelet yeW_oV_-1, where V_-1 = V_o + W_o. Representing this in the frequency domain:

f(w) = f(w/2) * H(w/2)
y(w) = y(w/2) * G(w/2)      (6)

H and G are the frequency responses of quadrature mirror filters satisfying the equation,

|H(w)|² + |G(w)|² = 1      (7)

From the filter-bank interpretation of wavelet transforms the wavelet should "match" the signal of interest such that the output of the matched filter bank is maximized. The wavelet transform is a correlation transform and the correlation is maximized when the spectrum of the wavelet (for the corresponding value of dilation and translation) is matched to the spectrum of the signal of interest. The wavelet is assumed to have zero phase.

F(w) = K * y_{j, k}(w)      (8)

Therefore, we develop a method for matching the amplitude spectrum of the wavelet to that of the desired signal. You cannot design a wavelet matched to the given signal by merely satisfying the OMRA conditions on the wavelet alone, as this would not necessarily yield a scaling function that generates an OMRA. Therefore conditions for an orthonormal MRA on the scaling function are propagated onto the wavelet. Wavelet matching under these conditions is performed, and the wavelet determines the scaling function, satisfying the conditions for an orthonormal MRA.

The first step in deriving the OMRA conditions for the wavelet spectrum amplitude is providing a means of deriving the scaling function from the mother wavelet. From the conditions provided earlier in this article and from Equations 6 and 7, we get the following expression (Refer to Chapa for proof):

|f(w)|² = |f(2w)|² + |y(2w)|²      (9)

Repeated substitution of |f(2^jw)|² into Equation 9 gives the following closed-form solution:

|f(w)|² = S_j=1|y(2^jw)|²;    w 0      (10)

To guarantee a closed-form solution, the scaling function is assumed to be band-limited with support over the interval:

we[-w_{m, Wm}]    w_m < p + a.0 < a < p/3      (11)

The proof of the above statement can be found in Chapa and Rao. Thus a band-limited orthonormal scaling function spectrum has support that can vary from [-p,p] to [-4p/3, 4p/3] with the structure given as:

|f(w)| = 1,    for |w| < p-a
|f(w)|² + |f(2p-w)|² = 1    for p-a < |w| < p+a      (12)

Because the orthonormal conditions are in terms of the magnitude of f only and since is f is band limited, the phase of f can assume any function. The corresponding orthonormal wavelet has support |w| e [p-a, 2p+2a] where a p/3. This is proved in Chapa and Rao and the results are summarized in Equation 13.

	{0,	0 < |w| < p-a
|y(w)| =	|f(2p-w)|	p-a < |w| < p+a
	1	p+a < |w| < 2p-2a	(13)
	|f(w/2)|	2p-2a < |w| < 2p+2a
	0	2p+2a < |w|
	}

The construction of |y(w)| exhibits a symmetry about 4p/3 (refer to Chapa). This leads to the following sufficient conditions on |y|:

|y(w)| = |y(4p-2w)|	for p-a < |w| < 4p/3
|y(w)|² + |y(2p-w)|² = 1	for p-a < |w| < 4p/3	(14)

The spectrum amplitude of an orthonormal wavelet has support that can vary from |w| e [p, 2p] as the lower limit to |w| e [2p/3, 8p/3] as the upper limit. We assume a = p/3 for the best possible match.

The next step is to develop the necessary and sufficient condition on y that guarantees orthonormality. Let

Y(k) = |y(kDw)|² k Z      (15)

Where Dw = 2p/2^L. The necessary and sufficient condition on Y to guarantee that |f(k)| generates an orthonormal MRA is found by substituting Equation 15 into the Poisson summation, Equation 5.

SLp=0 Sm=-	Y (2L/2P (k+2L+1m) = 1
	2L-1/3 < |2L/2P (k+2L+1m)| < 2L+2/3	(16)

Dw = 2p/2^L is the sample spacing of f(kDw) and y(kDw). In our implementation, we have assumed L=4 and a 64 point magnitude spectrum. This gives us 16 points in the wavelet passband between 2p/3 and 8p/3. By substituting different values of k into Equation 16, we obtain a set of constraint equations. Eliminating any equations that might repeat, we get a set of L unique constraint equations. We can reduce these equations to a constraint matrix A such that:

AY = 1      (17)

A is a 16 * L matrix. Elements of A may take values of 0, 1, or 2 only since a particular frequency may appear in Equation 16 only once and Y(k) is an even function.

A = {a_{i k} {0, 1, 2}; i = 1... L; j = 1 ... 2^L}      (18)

with 1 being a L x 1 vector given by 1 = {1 1 ...1}^T

Figure 1 shows the constraint matrix for 16 points in the wavelet passband.

Figure 1:  Constraint Matrix for 16 points in the wavelet passband

Let W and Y be vectors containing the samples of |F(kDw)|² and |y(kDw)|² respectively, in the passband.

W = {|F(kDw)|²;  k = [2^L/3]... [2^L+2/3]}
Y = {|y(kDw)|²;  k = [2^L/3]... [2^L+2/3]}      (19)

Where F(w) is the spectrum of the signal for which we desire a matched wavelet and y(w) is the matched wavelet spectrum. The error we need to minimize is given by:

E = (W - aY)^T(W - aY) / W^TW      (20)

This equation reduces to a constrained optimization problem with the error given by Equation 20 and the constraint equation by Equation 21 optimal wavelet spectrum is then given as,

Y = 1/a * W + A^T(AA^T)^-1 (1 - 1/a * AW)      (21)

Where

a = 1^T(AA^T)^-1AW / 1^T(AA^T)^-11      (22)

Figure 2:  Matched wavelet spectrum for D4 wavelet

Figure 2 shows the matched wavelet amplitude spectrum for the Daubechies D4 wavelet. The D4 wavelet satisfies the conditions for an MRA, hence the obtained match is quite accurate. However, since the D4 is not a band-limited wavelet, the match is not perfect. Figure 3 shows the result we obtain for the Meyer wavelet, which is both band-limited and satisfies the conditions for an MRA. As seen, we get perfect matching.

Figure 3:  Power spectrum for Meyer wavelet

Figure 4:  Radon transform

The next step is to take the Fourier transform of each of the projections and dilate the signals so that maximum of their energy comes into the passband. Figure 5 shows the results of amplitude matching after applying the wavelet constraints and Poisson equation discussed earlier in this paper. The shape of the wavelet is the best possible match of the wavelet spectrum to that of the projection. Since the projection signal does not satisfy the conditions for an OMRA, the match is not as good as in the case of the D4 wavelet.

Figure 5:  Matched spectrum for 45º projection

We then take the inverse transform of the designed wavelet.

Figure 6:  Matched wavelet for 45º projection

Next, the CWT (Continuous Wavelet Transform) of each projection is calculated using the Matched Wavelet for that particular projection. This is referred to as the Wavelet Radon Transform. The scale corresponding to the location of maximum power for each of the eight projections is stored in an 8X1 vector S_max. This gives an indication of the orientation of the object in the training image. The direction with maximum scale corresponds to the longest dimension of the object, the direction with minimum scale corresponds to the shortest dimension of the object, and so on. If we are now given a new image containing the same object but with different size and rotation, the S_max vector is a scaled and circularly shifted version of the vector obtained in the training image.

Figure 7:  Reconstruction by back projection

At the point of intersection of the peaks the target is present. Thus, the accurate size and orientation of the object is seen in the back projected image.

S = (S_max * S'_max)/max(S_max * S'_max)

Here S is normalized by the maximum value of the circular correlation. The direction of the object is determined from the location of the peak in S.

Figure 8:  Training image

Figures 9 through 11 show images containing the same object with different rotation/zoom and added noise. The ellipse indicates the size and rotation of the object. The angles indicate clockwise rotation.

Figure 9:  Detected object angle = 67.5º

Figure 10:  Detected object angle = 157.5º

Figure 11:  Detected object angle = 67.5º

This algorithm gives only the orientation; it cannot distinguish between, say, 0º and 180º.