Journal of Engineering

Volume 2013 (2013), Article ID 625364, 8 pages

http://dx.doi.org/10.1155/2013/625364

Research Article

Gamma Splines and Wavelets

Hannu Olkkonen¹ and Juuso T. Olkkonen²

¹Department of Applied Physics, University of Eastern Finland, P.O. Box 1627, 70211 Kuopio, Finland

²VTT Technical Research Centre of Finland, 02044 Espoo, Finland

Received 28 January 2013; Revised 15 April 2013; Accepted 15 April 2013

Academic Editor: Nizar Bouguila

Copyright © 2013 Hannu Olkkonen and Juuso T. Olkkonen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In this work we introduce a new family of splines termed as gamma splines for continuous signal approximation and multiresolution analysis. The gamma splines are born by -times convolution of the exponential by itself. We study the properties of the discrete gamma splines in signal interpolation and approximation. We prove that the gamma splines obey the two-scale equation based on the polyphase decomposition. to introduce the shift invariant gamma spline wavelet transform for tree structured subscale analysis of asymmetric signal waveforms and for systems with asymmetric impulse response. Especially we consider the applications in biomedical signal analysis (EEG, ECG, and EMG). Finally, we discuss the suitability of the gamma spline signal processing in embedded VLSI environment.

1. Introduction

Processing discrete-time signals and images requires the interpolation, decimation, and approximation procedures. Signal approximation via the B-spline transform is an established tool in a variety of signal management procedures based on discrete-time filtering, which are summarized in two excellent review articles [1, 2]. Compared to the conventional sinc interpolation the B-spline algorithms are extremely fast and suitable for microprocessor and VLSI environment [3–8]. The key feature in B-spline signal processing is a link between the continuous and discrete time domains. The approximation of signals by B-splines is based on the assumption of the underlying continuous signal, though we process discrete-time signal sequences.

In many applications of the B-spline signal processing the disadvantage comes from the symmetrical shape of the B-spline. The impulse response of the system is not usually symmetric but owns an exponential tail or resembles a damping sine wave. In this case the B-spline fit is not perfect.

In this work we introduce a new family of splines termed as gamma splines for signal approximation. The shape of the gamma spline is nonsymmetric having an exponential tail. We apply gamma splines to the numerical signal processing, shift invariant wavelet analysis, and multiresolution representation of the continuous-time signals and images.

2. Theoretical Considerations

2.1. B-Splines

The B-spline approximation (transform) of the signal is based on the convolution [3]

where is the scale sequence. We use brackets to emphasize the discrete-time signal. The integer denotes the discrete-time , where the sampling interval is normalized to .

The B-spline of the order is constructed by convolving an indicator function by itself

where

The analytic form of the B-spline is

where the unit step function is defined as

An interesting property of the B-splines is that they obey the two-scale equation [3] as follows:

where is the binomial kernel.

2.2. Discrete B-Splines

The Laplace transform of the B-spline comes from

The discrete B-spline equals to the continuous B-spline at integer values of time. Hence, the Laplace transform (7) and the -transform of the discrete B-spline should have inverse transforms which coincide at integer values in the time domain. Using the relation

we obtain the -transform of the discrete B-spline as follows:

where

We have . By differentiating with respect to , we obtain a recursion

Table 1 gives the and the discrete B-splines for to .

Table 1: and the discrete B-spline for to .

1

1

2

3

4

5

6

2.3. Wavelet Transform

Wavelet transform has gained a general acceptance in the multiscale analysis and synthesis of signals and images in various areas of science and technology. The family of wavelets for is generated from one-single function by dilation and translation operations as follows [9]

which constitute an orthonormal basis. The wavelet transform of a signal is defined by the analysis/synthesis equations

where are the wavelet coefficients in the subscale . The connection of the B-spline approximation (1) and the wavelet analysis/synthesis equations (13) is obvious. The B-spline approximation transfers the information to the next subscale via the two-scale equation (6). The wavelet behaves like the B-spline via the dilation and translation equation (12). If we keep the B-spline as a wavelet, the binomial kernel works as a scaling filter in the wavelet analysis.

3. Gamma Splines

3.1. Gamma Splines

In this work we introduce that the gamma spline representation

is born by -times convolution of the exponential by itself. The term gamma spline originates from the gamma variate defined in the literature as . For clarity, with the change of the variables to we obtain the gamma spline representation (14). Two typical gamma spline waveforms are described in Figure 1. The delayed version of the gamma spline is

where

The maximum value of the delayed gamma spline is obtained at time and the area . The continuous-time signal can be approximated by a weighted sum of delayed gamma splines as

where is the scale sequence. We may compute the scale sequence directly in time domain from the recursion

In the following we introduce shortly the most essential gamma spline signal processing tools.

Figure 1: Two typical gamma spline waveforms.

Interpolation:

Decimation:

Derivative:

Integral:

which follows from the integration by parts

The previous computational aids are based on the convolution in the time domain. The drawback comes from the time consuming computations for long data sequences, since the number of multiplications increases linearly as a function of the signal length. The replacement of the convolution by discrete-time filtering yields fast algorithms. In the following we construct the discrete gamma spline filter and apply it to the signal interpolation and approximation.

3.2. Discrete Gamma Splines

The -transform of the gamma spline (14) may be conducted by differentiating

which yields the recursion

This yields a general form

where is a polynomial in , which is given in Table 2 for . can be written in cascade realization as follows:

The implementation of the gamma spline is possible as a discrete-time IIR filter in the case the pole . If the pole lies outside the unit circle, the inverse filtering procedure (Appendix) has to be used.

Table 2: polynomial for to .

0
1

1

2

3

4

5

3.3. Scale Sequence

When we apply (17), for example, in interpolation or decimation we have to determine the scale sequence . In the following we present a feasible discrete-time filter for the computation of the scale coefficient sequence. The -transform of (17) gives

which yields

The denominator polynomial requires special treatment, since it is possible that all the roots are not inside the unit circle. However, we show in the Appendix that both the FIR and IIR filtering methods can be used for implementation of the polynomial.

3.4. Gamma Spline Interpolation Filter

The basic idea in gamma spline signal processing is the assumption of the underlying continuous signal, though we process discrete signal sequences. If the -transform of the discrete signal is , the interpolation of the discrete-time signal yields , which lacks the odd sequence. When we consider interpolation of the continuous signal , the procedure yields also the odd (approximated intermediate) signal values. Based on (28) we have the -transform of the interpolated signal as follows:

where

Due to (24) we have

where

which can be computed from the recursion

From (30) we obtain

where

Finally we have

We may interpret the as the gamma spline interpolation filter, which has the polyphase components and . The application of requires the computation of the scaling sequence via (29). We may observe that the polyphase component is responsible for the approximation of the intermediate points between integer values of time.

3.5. Fractional Time-Shift Gamma Spline Filter

Next we introduce the gamma spline filter, which produces the fractional time-shift for the analyzed signal. The transform of the shifted gamma spline comes from

Using the binomial expansion

we have

The previous result indicates that the time-shifted gamma spline can be represented by a linear combination of the gamma spines. The time-shifted signal is now obtained as

For example, the derivative of the signal can be calculated by differentiating with respect to as follows:

Many systems utilize the fractional delay filters, where the time-shift is negative. However, we may modify the time-shift algorithm (41) in the following way:

where the delay . Let us now define , which yields

which can be computed by applying the original time-shift algorithm (41).

4. Gamma Spline Wavelet Transform

The general tree structured dual-channel discrete wavelet transform is illustrated in Figure 2. The analysis part consists of the low-pass scaling and the high-pass wavelet filter . In the synthesis part the reconstruction filters and are related to the analysis filters by the perfect reconstruction (PR) condition as follows:

The and are defined as

In this work we apply the following essential result obtained in our previous work [4], which concerns on the PR condition (45).

Figure 2: The analysis and synthesis parts of the tree structured discrete wavelet transform. The decomposition consists of the wavelet sequences and one scaling sequence .

Lemma 1. If the scaling filter , the wavelet filter , and the reconstruction filters are related according to the PR condition (45), the following modified analysis and synthesis filters obey the PR condition as follows:

where is any polynomial in and is its inverse. The result can be proved by direct insertion of (47) into (45).

Definition 2. The gamma spline wavelet transform consists of the filter bank

We may verify that by fixing and and following the result of the Lemma 1 the filter bank (48) obeys the PR condition (45). The gamma spline can be implemented by cascade of   IIR filters having exponential impulse response (27). The inverse filter

where is realizable by a cascade of FIR filters having the impulse response .

The implementation of the inverse filter is described in the Appendix. The gamma spline scaling filter in (48) produces the scale sequence via (28). In multiscale analysis the scaling coefficients (Figure 1) are fed to the following stage. The gamma spline wavelet works as a high-pass filter.

5. Shift Invariant Gamma Spline Wavelet Transform

Main drawback in wavelet analysis is the dependence of the total energy of the wavelet coefficients on the fractional time shift of the analysed signal. Let us suppose that we have a discrete-time signal and the corresponding time shifted signal ; where , there occurs a significant difference in energy of the wavelet coefficients as a function of the time shift. In a nearly shift invariant method the real and imaginary parts of the complex wavelet coefficients are approximately a Hilbert transform pair [10]. The energy of the wavelet coefficients equals to the envelope, which provides smoothness and approximate shift invariance. In two parallel conjugate quadrature filter banks, where the impulse responses of the scaling filters are half-sample delayed versions of each other: and , the wavelet bases are a Hilbert transform pair [11]. The corresponding relation in biorthogonal filter banks is obtained using B-spline half-delay filters [4, 8]. In this work we introduce the Hilbert transform filter based on the gamma spline representation and construct the shift invariant gamma spline wavelet transform.

5.1. Hilbert Transform Filter

We define the Hilbert transform filter , which has the frequency response

where for and for . The construction of is based on the half-sample delay filter

whose frequency response is . The quadrature mirror filter has the frequency response , correspondingly. Now the frequency response of the filter is

Comparing (50) and using the IIR filter notation (51) we obtain the Hilbert transform filter as

The half-delay filter in (51) can be realized by the polyphase components of the gamma spline interpolation filter , which yields

5.2. Shift Invariant Wavelet Transform

The Hilbert transform filter is inserted in the BF bank using the result of Lemma 1 (47). The modified gamma spline wavelet transform filter bank is

The filter bank (33) can be highly simplified by noting the following equivalents:

By inserting (35) in (33) we obtain a highly simplified filter bank

The modified filter bank (57) can be realized by the Hilbert transform filter , which works as a prefilter for the analysed signal. On the other hand, the Hilbert transform filter works as a postfilter in the reconstruction stage.

Two parallel wavelet transforms can be realized by a single tree by defining the complex Hilbert transform operator

By filtering the real-valued signal by the Hilbert transform operator results in an analytic signal

whose magnitude response is zero at negative side of the frequency spectrum

The wavelet sequence is obtained by decimation of the high-pass filtered analytic signal

The result (61) means that the decimation does not produce aliasing, but the frequency spectrum is dilated by two. The frequency spectrum of the undecimated wavelet sequence contains frequency components in the range , but the frequency spectrum of the decimated analytic signal has the frequency band . Hence, the decimation does not produce overlapping and leakage to the negative frequency range. The Fourier transform of the decimated wavelet sequence of the fractionally delayed signal is . The energy of the decimated wavelet coefficients is , which does not depend on the fractional delay . Hence, the wavelet coefficients are shift invariant in respect to their energy content.

6. Discussion

The symmetric property of the compactly supported that B-splines is advantageous in many real-time applications, such as video and image compression and enhancement, as it reduces the complexity and machine run time. The B-spline signal processing is used vastly in CAD and other computer graphics for reproducing functions from only a few control points and/or boundary conditions. B-splines involve the parametric and geometric continuities, which, for example, improve the reconstruction performance of the 3D tomography.

In this framework we introduce the gamma splines, which are born by -times convolution of the exponential by itself. The main reason for the selection of the exponential kernel function is that it can be described uniquely both in continuous and discrete-time domains by the single-pole functions and . Depending on an implementation the application of the gamma spline parameters and can be optimally adjusted. In many applications we have selected , which means that the maximum is attained at . However, only the discrete-time domain solution of the scale sequence (18) requires this condition. The cascade realization (27) permits an extremely fast implementation in embedded VLSI environment, even in the case the single-pole lies outside the unit circle. The Appendix describes two variants of the time inversed filtering procedure.

It should be pointed out that the gamma splines are not compactly supported. By increasing the value of the parameter the shape of the gamma spline becomes more symmetrical, and the time support decreases. The parameter affects the initial slope of the gamma spline and the exponential damping rate. At very high values of the gamma spline approaches the delta function, and (17) is close to the conventional Shannon’s sampling theorem. The PR property of the shift invariant gamma spline wavelet bank (57) is not dependent on the selection of the and parameters. Hence, we may easily construct adaptive gamma spline wavelets, where and parameters are self-adjusting according to some cost criteria or error function. Recently our research group introduced an adaptive matrix-vector gradient algorithm [12], which can be used for optimization of the and parameters to match for special signal features. Our preliminary results in wireless transmission of the ambulatory ECG indicate that the adaptive gamma spline wavelets improved compression performance compared with the static B-spline wavelets. The clinical results will be published elsewhere [13].

In wavelet analysis the smoothness of the scaling filter is an important feature, which would allow shift invariance [10]. In this work we introduced the shift invariant gamma spline wavelet transform, which is based on a Hilbert transform filter (54). It appeared that the modified filter bank (36) can be realized by the Hilbert transform filter, which works as a prefilter for the analysed signal and as a postfilter in the reconstruction stage, respectively. By defining the complex Hilbert transform operator (58), two parallel wavelet transforms can be realized by a single-tree structured transform (Figure 2).

The present gamma spline formulation offers a plenty of new tools for signal and image processing. For example, in biomedical signal analysis EEG, ECG, and EMG signals are not symmetric but owe an exponential tail or resemble damping sinusoidal waveforms. We have previously shown that the shift invariant B-spline wavelets are useful in the subscale analysis of EEG signals [8]. According to our experience the shift invariant gamma spline wavelets fit somewhat better to the neuroelectric spikes compared with the biorthogonal B-spline wavelets. However, a larger medical study is required to warrant this preliminary observation.

In our previous work we introduced a novel wavelet excitation method for measurement of the system transfer function [14]. The method employs system excitation by wavelet shaped waveforms instead of commonly used impulse or sinusoidal excitation. The system transfer function can be reconstructed from the output measurements. Data acquisition can be designed, so that if different excitation sequences are used and the excitation rate is , the sampling rate can be reduced to . Gamma splines can be realized by a cascade of identical RC filters [15]. Hence, the system can be excited by analog waveforms instead of the zero-order hold signal produced by the digital-to-analog converter. The gamma spline excitation method permits high speed applications, where the sampling rate may be considerably lower compared with system bandwidth. The method is especially advantageous in testing systems, where impulse or sinusoidal excitation cannot be applied.

Appendix

Implementation of the Inverse Filter

The inverse filter can be generally factored into the product of the filter as follows:

where is a constant and are the roots of the polynomial. Let us consider a filter , whose root lies outside the unit circle. By replacing by we have

which has the impulse response

The length of filter has been selected, so that the coefficient is negligible. Now can be implemented by the FIR filter having the impulse response

The disadvantage of the FIR implementation is that if the root is close to the unit circle, the sequence (A.3) vanishes relatively slowly. A competitive method is the inverse IIR filtering procedure. If the input signal is and the output signal is we have

The final output is obtained by reversing the resulting sequence. The following MATLAB program rfilter.m demonstrates the method

Acknowledgements

The authors are greatly indebted to the two anonymous reviewers, whose comments improved the original paper significantly. This work was supported by the National Technology Agency of Finland (TEKES).

References

M. Unser, A. Aldroubi, and M. Eden, “B-spline signal processing—part I: theory,” IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821–833, 1993.

M. Unser, A. Aldroubi, and M. Eden, “B-spline signal processing—part II: efficient design and applications,” IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 834–848, 1993.

J. T. Olkkonen and H. Olkkonen, “Fractional delay filter based on the B-spline transform,” IEEE Signal Processing Letters, vol. 14, no. 2, pp. 97–100, 2007.

H. Olkkonen and J. T. Olkkonen, “Half-delay B-spline filter for construction of shift-invariant wavelet transform,” IEEE Transactions on Circuits and Systems II, vol. 54, no. 7, pp. 611–615, 2007.

J. T. Olkkonen and H. Olkkonen, “Fractional time-shift B-spline filter,” IEEE Signal Processing Letters, vol. 14, no. 10, pp. 688–691, 2007.

H. Olkkonen and J. T. Olkkonen, “Shifted linear interpolation filter,” Journal of Signal and Information Processing, vol. 1, no. 1, pp. 44–49, 2010.

H. Olkkonen and J. T. Olkkonen, “Shifted B-spline interpolation filter,” Current Development on Theory and Applications of Wavelets, vol. 4, no. 3, pp. 289–305, 2010.

H. Olkkonen and J. T. Olkkonen, “Shift-invariant B-spline wavelet transform for multi-scale analysis of neuroelectric signals,” IET Signal Processing, vol. 4, no. 6, pp. 603–609, 2010.

I. Daubechies, “Orthonormal bases of compactly supported wavelets,” Communications on Pure and Applied Mathematics, vol. 41, pp. 909–996, 1988.

N. Kingsbury, “Complex wavelets for shift invariant analysis and filtering of signals,” Journal of Applied and Computational Harmonic Analysis, vol. 10, no. 3, pp. 234–253, 2001.

I. W. Selesnick, “The design of approximate Hilbert transform pairs of wavelet bases,” IEEE Transactions on Signal Processing, vol. 50, no. 5, pp. 1144–1152, 2002.

J. T. Olkkonen, S. Ahtiainen, K. Jarvinen, and H. Olkkonen, “Adaptive matrix/vector gradient algorithm for design of IIR filters and ARMA models,” Journal of Signal and Information Processing, vol. 4, no. 2, 2013.

K. Jarvinen, [Ph.D. dissertation], unpublished.

H. Olkkonen and J. T. Olkkonen, “Wavelet excited measurement of system transfer function,” Review of Scientific Instruments, vol. 78, no. 2, Article ID 025104, 5 pages, 2007.

H. Olkkonen, S. Ahtiainen, K. Jarvinen, and J. T. Olkkonen, “Reconstruction of impulse train by parallel sampling of cascaded identical RC filters,” Wireless Engineering and Technology, vol. 4, no. 2, pp. 71–76, 2013.