Gamma Splines and Wavelets

In this workwe introduce a new family of splines termed as gamma splines for continuous signal approximation andmultiresolution analysis.The gamma splines g m,b (t) = t m e −bt , t ≥ 0 are born by (m+1)-times convolution of the exponential e−bt 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.


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][4][5][6][7][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.

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  = 1.
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.

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  1 () = 1/(1 −  −1 ).By differentiating with respect to , we obtain a recursion Table 1 gives the   () and the discrete B-splines   () for  = 1 to 6.
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.

Gamma Splines
3.1.Gamma Splines.In this work we introduce that the gamma spline representation is born by ( + 1)-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  =  and  = 1/ 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  max =  + / and the area (1/) +1 Γ( + 1).
The continuous-time signal () can be approximated by a weighted sum of delayed gamma splines as where [],  = 0, 1, 2, . . . is the scale sequence.We may compute the scale sequence directly in time domain from the recursion [1] , . . .
In the following we introduce shortly the most essential gamma spline signal processing tools. 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.

Discrete Gamma Splines.
The -transform of the gamma spline ( 14) may be conducted by differentiating which yields the recursion Table 2:   () polynomial for  = 0 to 5.
This yields a general form where   () is a polynomial in  −1 , which is given in Table 2 for  = 0 to 5.   () 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  − < 1.If the pole lies outside the unit circle, the inverse filtering procedure (Appendix) has to be used.
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  −1  () polynomial.

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 ( 2 ), 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   ( 2 ) and   ( 2 ).The application of   () requires the computation of the scaling sequence () via (29).We may observe that the polyphase component   ( 2 ) is responsible for the approximation of the intermediate points between integer values of time.

Fractional Time-Shift Gamma Spline
Filter.Next we introduce the gamma spline filter, which produces the fractional time-shift Δ ∈ [0, 1] 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 timeshift algorithm (41) in the following way: where the delay  ∈ [0, 1].Let us now define Δ = 1−, which yields which can be computed by applying the original time-shift algorithm (41).

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  0 () and the high-pass wavelet filter  1 ().In the synthesis part the reconstruction filters  0 () and  1 () are related to the analysis filters by the perfect reconstruction (PR) condition as follows: The  0 () and  1 () 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).

Lemma 1.
If the scaling filter  0 (), the wavelet filter  1 (), and the reconstruction filters  0 () 1 () 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  −1 and  −1 () 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  0 () = 1 and  1 () =  −1 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  + 1 IIR filters having exponential impulse response (27).The inverse filter where  −1  () is realizable by a cascade of  FIR filters having the impulse response () = 1 −  −  −1 .

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  ∈ [0,1], 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: ℎ 0 [] and ℎ 0 [ − 0.5], 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.

Hilbert Transform Filter.
We define the Hilbert transform filter H(), which has the frequency response where sgn() = 1 for  ≥ 0 and sgn() = −1 for  < 0. The construction of H() is based on the half-sample delay filter whose frequency response is () =  −/2 .The quadrature mirror filter (−) has the frequency response ( − ) =  −(−)/2 , 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

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 H(), which works as a prefilter for the analysed signal.On the other hand, the Hilbert transform filter H(−) 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 highpass filtered analytic signal ) . (61) 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 0 ≤  < , but the frequency spectrum of the decimated analytic signal has the frequency band 0 ≤  < 2.Hence, the decimation does not produce overlapping and leakage to the negative frequency range.

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 Bspline 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 singlepole functions 1/( + ) and 1/(1 +  −  −1 ).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  = 1.However, only the discretetime 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 Bspline 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.

Figure 2 :
Figure 2: The analysis and synthesis parts of the tree structured discrete wavelet transform.The decomposition consists of the wavelet sequences   [],  = 1, 2, . . .,  and one scaling sequence   [].
The Fourier transform of the decimated wavelet sequence [ − /2] of the fractionally delayed signal [ − ] is  −/2   (/2)/2.The energy of the decimated wavelet coefficients is |(/2)|/2, which does not depend on the fractional delay .Hence, the wavelet coefficients are shift invariant in respect to their energy content.