Gamma Splines and Wavelets

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 𝑔_𝑚,𝑏(𝑡)=𝑡^𝑚𝑒^−𝑏𝑡,𝑡≥0 are born by (𝑚+1)-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] 𝑥(𝑡)=_𝑘𝑐[𝑘]𝛽_𝐾(𝑡−𝑘), (1) 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 𝛽_𝐾(𝑡)=𝛽(𝑡)∗𝛽(𝑡)∗⋯𝛽_𝐾times(𝑡), (2) where 𝛽(𝑡)=10≤𝑡≤1,0elsewhere. (3) The analytic form of the B-spline is 𝛽_𝐾(𝑡)=1(𝐾−1)!^𝐾_𝑘=0(−1)^𝑘𝐾𝑘(𝑡−𝑘)^𝐾−1𝑢(𝑡−𝑘), (4) where the unit step function is defined as 𝑢(𝑡)=1for𝑡≥0,0for𝑡<0. (5) An interesting property of the B-splines is that they obey the two-scale equation [3] as follows: 𝛽_𝐾𝑡2=2_𝑘𝐵_𝐾[𝑘]𝛽_𝐾(),𝐵𝑡−𝑘_𝐾[𝑘]=12^𝐾𝐾𝑘=𝐾!2^𝐾()!,𝑘!𝐾−𝑘 (6) where 𝐵_𝐾[𝑘] is the binomial kernel.

2.2. Discrete B-Splines

The Laplace transform of the B-spline comes from 𝐿𝛽(𝑡=1)𝑠1−𝑒^−𝑠⟹𝛽_𝐾(𝑠)=1𝑠^𝐾1−𝑒^−𝑠𝐾. (7) 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 𝐿⁻¹1𝑠^𝐾=𝑡^𝐾−1()!𝑢(𝑡)𝐾−1 (8) we obtain the 𝑧-transform of the discrete B-spline as follows: 𝛽_𝐾(𝑧)=𝑍𝛽_𝐾[𝑘]=𝑍𝐿⁻¹1𝑠^𝐾1−𝑒^−𝑠𝐾=𝑁_𝐾(𝑧)1−𝑧^−1𝐾, (9) where 𝑁_𝐾(𝑧)=𝑍𝐿⁻¹1𝑠^𝐾=^∞_𝑘=0𝑘^𝐾−1()!𝑧𝐾−1^−𝑘. (10) We have 𝑁₁(𝑧)=1/(1−𝑧⁻¹). By differentiating with respect to 𝑧, we obtain a recursion 𝑁_𝐾+1(𝑧)=−𝑧𝐾𝑑𝑁_𝐾(𝑧).𝑑𝑧 (11) Table 1 gives the 𝑁_𝐾(𝑧) and the discrete B-splines 𝛽_𝐾(𝑧) for 𝐾=1 to 6.

Table 1 

𝑁_𝐾(𝑧) and the discrete B-spline 𝛽_𝐾(𝑧) for 𝐾=1 to 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] 𝜓_𝑚,𝑛(𝑡)=2^−𝑚/2𝜓2^−𝑚𝑡−𝑛 (12) which constitute an orthonormal basis. The wavelet transform of a signal 𝑥(𝑡) is defined by the analysis/synthesis equations 𝑥^𝑚_𝑛=^+∞_−∞𝑥(𝑡)𝜓_𝑚,𝑛(𝑡)𝑥(𝑡)𝑑𝑡,=_𝑚𝑛𝑥^𝑚_𝑛𝜓_𝑚,𝑛(𝑡), (13) 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 𝑔_𝑚,𝑏(𝑡)=𝑡^𝑚𝑒^−𝑏𝑡𝑢(𝑡) (14) 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 𝑔_𝑚,𝑏(𝑡−𝑑)=(𝑡−𝑑)^𝑚𝑒^{−𝑏(𝑡−𝑑)}𝑢()=𝑡−𝑑_𝑚𝑘=0𝑎_𝑘,𝑑𝑔_𝑘,𝑏(𝑡)𝑢(𝑡−𝑑), (15) where 𝑎_𝑘,𝑑=𝑚𝑘(−𝑑)^𝑚−𝑘𝑒^𝑏𝑑. (16) 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 𝑥(𝑡)=_𝑘𝑐[𝑘]𝑔_𝑚,𝑏(𝑡−𝑘), (17) where 𝑐[𝑘],𝑘=0,1,2,… is the scale sequence. We may compute the scale sequence directly in time domain from the recursion 𝑥[1]=𝑐[0]𝑔_𝑚,𝑏[1][0]=𝑥[1]⟹𝑐𝑔_𝑚,𝑏[1],𝑥[2]=𝑐[0]𝑔_𝑚,𝑏[2]+𝑐[1]𝑔_𝑚,𝑏[1]⟹𝑐[1]=𝑥[2]−𝑐[0]𝑔_𝑚,𝑏[2]𝑔_𝑚,𝑏[1],⋮𝑥[𝑛]=^𝑛−2_𝑘=0𝑐[𝑘]𝑔_𝑚,𝑏[𝑛−𝑘]+𝑐[𝑛−1]𝑔_𝑚,𝑏[1]⟹𝑐[]=𝑥[𝑛]𝑛−1−∑^𝑛−2_𝑘=0𝑐[𝑘]𝑔_𝑚,𝑏[𝑛−𝑘]𝑔_𝑚,𝑏[1]. (18) In the following we introduce shortly the most essential gamma spline signal processing tools.

Figure 1 

Two typical gamma spline waveforms.

Interpolation: 𝑥𝑡2=2^−𝑚_𝑘𝑐[𝑘]𝑔_𝑚,𝑏/2(𝑡−2𝑘). (19)

Decimation: 𝑥(2𝑡)=2^𝑚_𝑘𝑐[𝑘]𝑔_𝑚,2𝑏𝑘𝑡−2. (20)

Derivative: 𝑑𝑥(𝑡)𝑑𝑡=_𝑘𝑐[𝑘]𝑚𝑔_{𝑚−1,𝑏}(𝑡−𝑘)−𝑏𝑔_𝑚,𝑏(𝑡−𝑘)=𝑚_𝑘𝑐[𝑘]𝑔_{𝑚−1,𝑏}(𝑡−𝑘)−𝑏𝑥(𝑡). (21)

Integral: ^𝑡₀𝑥(𝑡)𝑑𝑡=_𝑘𝑐[𝑘]^𝑡₀𝑔_𝑚,𝑏(𝑡−𝑘)𝑑𝑡=𝑏⁻¹𝑚_𝑘𝑐[𝑘]^𝑡₀𝑔_{𝑚−1,𝑏}(𝑡−𝑘)𝑑𝑡−𝑏⁻¹𝑥(𝑡) (22) which follows from the integration by parts ^𝑡₀𝑔_𝑚,𝑏(𝑡−𝑘)𝑑𝑡=𝑏⁻¹𝑚^𝑡₀𝑔_𝑚−1(𝑡−𝑘)𝑑𝑡−𝑏⁻¹𝑔_𝑚,𝑏(𝑡−𝑘). (23) 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 𝐺_𝑚(𝑧)=^∞_𝑘=0𝑘^𝑚𝑒^−𝑏𝑘𝑧^−𝑘⟹𝑑𝐺_𝑚(𝑧)𝑑𝑧=−𝑧^∞−1_𝑘=0𝑘^𝑚+1𝑒^−𝑏𝑘𝑧^−𝑘 (24) which yields the recursion 𝐺₀(𝑧)=11−𝑒^−𝑏𝑧⁻¹⟹𝐺_𝑚+1(𝑧)=−𝑧𝑑𝐺_𝑚(𝑧).𝑑𝑧 (25) This yields a general form 𝐺_𝑚(𝑧)=𝑃_𝑚(𝑧)1−𝑒^−𝑏𝑧^−1𝑚+1, (26) where 𝑃_𝑚(𝑧) is a polynomial in 𝑧⁻¹, which is given in Table 2 for 𝑚=0to5. 𝐺_𝑚(𝑧) can be written in cascade realization as follows: 𝐺_𝑚(𝑧)=𝑃_𝑚(𝑧)^𝑚+1_𝑖=111−𝑒^−𝑏𝑧⁻¹. (27) 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.

Table 2 

𝑃_𝑚(𝑧) polynomial for 𝑚=0 to 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 𝑋(𝑧)=𝐶(𝑧)𝐺_𝑚(𝑧)⟹𝐶(𝑧)=𝐺_𝑚⁻¹(𝑧)𝑋(𝑧) (28) which yields 𝐶(𝑧)=𝑃_𝑚⁻¹(𝑧)^𝑚+1_𝑖=11−𝑒^−𝑏𝑧⁻¹𝑋(𝑧). (29) 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: 𝑍𝑥𝑛2=𝐶𝑧²𝑍𝑔_𝑚,𝑏𝑛2, (30) where 𝑍𝑔_𝑚,𝑏𝑛2=^∞_𝑘=0𝑘2^𝑚𝑒^{−𝑏𝑘/2}𝑧^−𝑘=^∞_𝑘=0𝑘^𝑚𝑒^−𝑏𝑘𝑧^−2𝑘+𝑧⁻¹×^∞_𝑘=01𝑘+2^𝑚𝑒^{−𝑏(𝑘+1/2)}𝑧^−2𝑘. (31) Due to (24) we have 𝑍𝑔_𝑚,𝑏𝑛2=𝐺_𝑚𝑧²+𝑧⁻¹𝐹_𝑚𝑧², (32) where 𝐹_𝑚𝑧²=^∞_𝑘=01𝑘+2^𝑚𝑒^{−𝑏(𝑘+1/2)}𝑧^−2𝑘 (33) which can be computed from the recursion 𝐹₀(𝑧)=𝑒^−𝑏/21−𝑒^−𝑏𝑧⁻¹⟹𝐹_𝑚+1(𝑧)=12𝐹_𝑚(𝑧)−𝑧𝑑𝐹_𝑚(𝑧).𝑑𝑧 (34) From (30) we obtain 𝑍𝑥𝑛2=𝐶𝑧²𝐺_𝑚𝑧²+𝑧⁻¹𝐹_𝑚𝑧²=𝐶𝑧²𝐼_𝑚(𝑧), (35) where 𝐼_𝑚(𝑧)=𝐺_𝑚𝑧²+𝑧⁻¹𝐹_𝑚𝑧². (36) Finally we have 𝑍𝑥𝑛2=𝑋𝑧²+𝑧⁻¹𝐶𝑧²𝐹_𝑚𝑧². (37) 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 Δ∈[0,1] for the analyzed signal. The 𝑧 transform of the shifted gamma spline comes from 𝑍𝑔_𝑚,𝑏[𝑘+Δ]=𝐺_𝑚,𝑏(Δ,𝑧)=^∞_𝑘=0(𝑘+Δ)^𝑚𝑒^{−𝑏(𝑘+Δ)}𝑧^−𝑘. (38) Using the binomial expansion (𝑘+Δ)^𝑚=^𝑚_𝑛=0𝑚𝑛𝑘^𝑚−𝑛Δ^𝑛 (39) we have 𝐺_𝑚,𝑏(Δ,𝑧)=𝑒^{𝑚−𝑏Δ}_𝑛=0𝑚𝑛Δ^𝑛𝐺_{𝑚−𝑛,𝑏}(𝑧). (40) 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 𝑍{𝑥[𝑘+Δ]}=𝐶(𝑧)𝐺_𝑚,𝑏(Δ,𝑧). (41) For example, the derivative of the signal can be calculated by differentiating with respect to Δ as follows: 𝑍𝑑𝑥[(𝑧)𝑒𝑑𝑡𝑘+Δ]=𝐶^{𝑚−𝑏Δ}_𝑛=0𝑚𝑛𝑛Δ^𝑛−1−𝑏𝐺_{𝑚−𝑛,𝑏}(𝑧). (42) 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: 𝑍{𝑥[𝑛−𝑑]}=𝑍{𝑥[𝑛−1+1−𝑑]}=𝑧⁻¹𝑍{𝑥[𝑛+1−𝑑]}, (43) where the delay 𝑑∈[0,1]. Let us now define Δ=1−𝑑, which yields 𝑍{𝑥[𝑛−𝑑]}=𝑧⁻¹𝑍{𝑥[𝑛+Δ]} (44) 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: 𝐻₀(𝑧)𝐻₁(𝑧)𝐻₀(−𝑧)𝐻₁()𝐺−𝑧⁰(𝑧)𝐺₁(𝑧)=2𝑧^−𝑘0. (45) The 𝐺₀(𝑧) and 𝐺₁(𝑧) are defined as 𝐺₀(𝑧)=𝐻₁(),𝐺−𝑧₁(𝑧)=−𝐻₀(−𝑧). (46) 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 𝑤𝑖[𝑛],𝑖=1,2,…,𝑛 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: 𝐇₀(𝑧)=𝐹(𝑧)𝐻₀(𝑧),𝐇₁(𝑧)=𝐹⁻¹(−𝑧)𝐻₁(𝑧),𝐆₀(𝑧)=𝐹⁻¹(𝑧)𝐺₀(𝑧),𝐆₁(𝑧)=𝐹(−𝑧)𝐺₁(𝑧), (47) 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 𝐻₀(𝑧)=𝐺_𝑚⁻¹(𝑧),𝐻₁(𝑧)=𝑧⁻¹𝐺_𝑚(),𝐺−𝑧₀(𝑧)=𝑧⁻¹𝐺_𝑚(𝑧),𝐺₁(𝑧)=𝐺_𝑚⁻¹(−𝑧). (48) We may verify that by fixing 𝐻₀(𝑧)=1 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  𝑚+1 IIR filters having exponential impulse response (27). The inverse filter 𝐺_𝑚⁻¹(𝑧)=𝑃_𝑚⁻¹(𝑧)^𝑚+1_𝑖=11−𝑒^−𝑏𝑧⁻¹=𝑃_𝑚⁻¹(𝑧)^𝑚+1_𝑖=1𝑅(𝑧), (49) where 𝐺_𝑚⁻¹(𝑧) is realizable by a cascade of 𝑚 FIR filters having the impulse response 𝑅(𝑧)=1−𝑒^−𝑏𝑧⁻¹.
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 𝑆(𝑧)=𝐶(𝑧)↓2 (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 𝜏∈[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: ℎ₀[𝑛] and ℎ₀[𝑛−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.

5.1. Hilbert Transform Filter

We define the Hilbert transform filter ℋ(𝑧), which has the frequency response ℋ(𝜔)=𝑒^{−𝑗𝜋/2}sgn(𝜔), (50) where sgn(𝜔)=1 for 𝜔≥0 and sgn(𝜔)=−1 for 𝜔<0. The construction of ℋ(𝑧) is based on the half-sample delay filter 𝐷(𝑧)=𝑧^−0.5=𝐴(𝑧)𝐵(𝑧) (51) whose frequency response is 𝐷(𝜔)=𝑒^{−𝑗𝜔/2}. The quadrature mirror filter 𝐷(−𝑧) has the frequency response 𝐷(𝜔−𝜋)=𝑒^{−𝑗(𝜔−𝜋)/2}, correspondingly. Now the frequency response of the filter 𝐷(𝑧)/𝐷(−𝑧) is 𝐷(𝜔)𝐷(𝜔−𝜋)=𝑒^{−𝑗𝜔/2}𝑒^{𝑗(𝜔−𝜋)/2}=𝑒^{−𝑗𝜋/2}. (52) Comparing (50) and using the IIR filter notation (51) we obtain the Hilbert transform filter as ℋ(𝑧)=𝐴(𝑧)𝐵(−𝑧)𝐵(𝑧)𝐴().−𝑧 (53) The half-delay filter 𝐷(𝑧) in (51) can be realized by the polyphase components of the gamma spline interpolation filter 𝐼_𝑚(𝑧), which yields 𝐷(𝑧)=𝑧⁻¹𝐹_𝑚(𝑧)𝐺_𝑚(𝑧)(𝑧)𝐹⟹ℋ=−^𝑚(𝑧)𝐺_𝑚(−𝑧)𝐺_𝑚(𝑧)𝐹_𝑚().−𝑧 (54)

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 𝐻₀(𝑧)=ℋ(𝑧)𝐺_𝑚⁻¹(𝑧),𝐻₁(𝑧)=ℋ⁻¹(−𝑧)𝑧⁻¹𝐺_𝑚(−𝑧),𝐺₀(𝑧)=ℋ⁻¹(𝑧)𝑧⁻¹𝐺_𝑚(𝑧),𝐺₁(𝑧)=ℋ(−𝑧)𝐺_𝑚⁻¹(−𝑧). (55) The filter bank (33) can be highly simplified by noting the following equivalents: ℋ⁻¹()(𝑧),ℋ−𝑧=ℋ₋₁(𝑧)=ℋ(−𝑧). (56) By inserting (35) in (33) we obtain a highly simplified filter bank 𝐻₀(𝑧)=ℋ(𝑧)𝐺_𝑚⁻¹(𝑧),𝐻₁(𝑧)=ℋ(𝑧)𝑧⁻¹𝐺_𝑚(−𝑧),𝐺₀(𝑧)=ℋ(−𝑧)𝑧⁻¹𝐺_𝑚(𝑧),𝐺₁(𝑧)=ℋ(−𝑧)𝐺_𝑚⁻¹(−𝑧). (57) 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 ℋ_𝑎(𝑧)=1+𝑗ℋ(𝑧). (58) By filtering the real-valued signal 𝑥[𝑛] by the Hilbert transform operator results in an analytic signal 𝑥_𝑎[𝑛]=𝑥[𝑛]+𝑗ℋ{𝑥[𝑛]} (59) whose magnitude response is zero at negative side of the frequency spectrum 𝑋_𝑎(𝜔)=2𝑋(𝜔)0≤𝜔<𝜋,0−𝜋≤𝜔<0. (60) The wavelet sequence is obtained by decimation of the high-pass filtered analytic signal 𝑊(𝜔)=𝑋_𝑎(𝜔)𝐻₁(𝜔)_↓2=𝑊_𝑎(𝜔)_↓2=12𝑋_𝑎𝜔2𝐻₁𝜔2. (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. 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.

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 1/(𝑠+𝑏) and 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 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: 𝑃_𝑛⁻¹(𝑧)=𝑐^𝑟_𝑖=111−𝑧_𝑖𝑧⁻¹=𝑐^𝑟_𝑖=1𝐻_𝑖(𝑧), (A.1) 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 𝐻_𝑖𝑧⁻¹=−𝑧_𝑖⁻¹𝑧⁻¹1−𝑧_𝑖⁻¹𝑧⁻¹ (A.2) which has the impulse response ℎ_𝑖[−𝑛]=0−𝑧_𝑖⁻¹−𝑧_𝑖⁻²⋯−𝑧_𝑖^−𝑚. (A.3) The length of filter has been selected, so that the coefficient −𝑧_𝑖^−𝑚−1 is negligible. Now 𝐻_𝑖(𝑧) can be implemented by the FIR filter having the impulse response ℎ_𝑖[𝑛]=−𝑧_𝑖^−𝑚−𝑧_𝑖^−𝑚+1−𝑧_𝑖^−𝑚+2⋯−𝑧_𝑖⁻¹. (A.4) 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 𝑌𝑧⁻¹=𝐻_𝑖𝑧⁻¹𝑈𝑧⁻¹=−𝑧_𝑖⁻¹𝑧⁻¹1−𝑧_𝑖⁻¹𝑧⁻¹𝑈𝑧⁻¹. (A.5) The final output is obtained by reversing the resulting sequence. The following MATLAB program rfilter.m demonstrates the method 𝑦=rfilter𝑢,𝑧_𝑖𝑢=𝑢(end∶−1∶1);𝑎=0−𝑧_𝑖⁻¹;𝑏=1−𝑧_𝑖⁻¹;();();𝑦=filter𝑎,𝑏,𝑢𝑦=𝑦end∶−1∶1 (A.6)

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).