# Gamma Splines and Wavelets

Authors: Hannu Olkkonen; Juuso T. Olkkonen
Journal: Journal of Engineering
 (2013)

---

## Abstract

In this work we introduce a new family of splines termed as gamma splines for continuous signal approximation and multiresolution analysis. The gamma splinesgm,b(t)=tme-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.


---

## Body

## 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 signalx(t) is based on the convolution [3]
(1)x(t)=∑kc[k]βK(t-k),
where c[k] is the scale sequence. We use brackets [] to emphasize the discrete-time signal. The integer k denotes the discrete-time t=kT, where the sampling interval is normalized to T=1.The B-splineβK(t) of the order K is constructed by convolving an indicator function β(t) by itself
(2)βK(t)=β(t)*β(t)*⋯β(t)︸Ktimes,
where
(3)β(t)={10≤t≤1,0elsewhere.
The analytic form of the B-spline is
(4)βK(t)=1(K-1)!∑k=0K(-1)k(Kk)(t-k)K-1u(t-k),
where the unit step function is defined as
(5)u(t)={1fort≥0,0fort<0.
An interesting property of the B-splines is that they obey the two-scale equation [3] as follows:
(6)βK(t2)=2∑kBK[k]βK(t-k),BK[k]=12K(Kk)=K!2Kk!(K-k)!,
where BK[k] is the binomial kernel.


### 2.2. Discrete B-Splines

The Laplace transform of the B-spline comes from(7)L{β(t)}=1s(1-e-s)⟹βK(s)=1sK(1-e-s)K.
The discrete B-spline βp[k] equals to the continuous B-spline at integer values of time. Hence, the Laplace transform (7) and the z-transform of the discrete B-spline should have inverse transforms which coincide at integer values in the time domain. Using the relation
(8)L-1(1sK)=tK-1(K-1)!u(t)
we obtain the z-transform of the discrete B-spline as follows:
(9)βK(z)=Z{βK[k]}=Z{L-1(1sK(1-e-s)K)}=NK(z)(1-z-1)K,
where
(10)NK(z)=Z{L-1(1sK)}=∑k=0∞kK-1(K-1)!z-k.
We have N1(z)=1/(1-z-1). By differentiating with respect to z, we obtain a recursion
(11)NK+1(z)=-zKdNK(z)dz.
Table 1 gives the NK(z) and the discrete B-splines βK(z) for K=1 to 6.Table 1
N
K
(
z
) and the discrete B-spline βK(z) for K=1 to 6.
1
1
1
-
z
-
1
1
2
z
-
1
(
1
-
z
-
1
)
2
z
-
1
3
z
-
1
+
z
-
2
2
(
1
-
z
-
1
)
3
z
-
1
+
z
-
2
2
4
z
-
1
+
4
z
-
2
+
z
-
3
6
(
1
-
z
-
1
)
4
z
-
1
+
4
z
-
2
+
z
-
3
6
5
z
-
1
+
11
z
-
2
+
11
z
-
3
+
z
-
4
24
(
1
-
z
-
1
)
4
z
-
1
+
11
z
-
2
+
11
z
-
3
+
z
-
4
24
6
z
-
1
+
26
z
-
2
+
66
z
-
3
+
26
z
-
4
+
z
-
5
120
(
1
-
z
-
1
)
4
z
-
1
+
26
z
-
2
+
66
z
-
3
+
26
z
-
4
+
z
-
5
120


### 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 form,n∈Z is generated from one-single function by dilation and translation operations as follows [9]
(12)ψm,n(t)=2-m/2ψ(2-mt-n)
which constitute an orthonormal basis. The wavelet transform of a signal x(t) is defined by the analysis/synthesis equations
(13)xnm=∫-∞+∞x(t)ψm,n(t)dt,x(t)=∑m∑nxnmψm,n(t),
where xnm are the wavelet coefficients in the subscale m. 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 ψm,n(t) behaves like the B-spline via the dilation and translation equation (12). If we keep the B-spline as a wavelet, the binomial kernel BK[n] works as a scaling filter in the wavelet analysis.


## 2.1. B-Splines

The B-spline approximation (transform) of the signalx(t) is based on the convolution [3]
(1)x(t)=∑kc[k]βK(t-k),
where c[k] is the scale sequence. We use brackets [] to emphasize the discrete-time signal. The integer k denotes the discrete-time t=kT, where the sampling interval is normalized to T=1.The B-splineβK(t) of the order K is constructed by convolving an indicator function β(t) by itself
(2)βK(t)=β(t)*β(t)*⋯β(t)︸Ktimes,
where
(3)β(t)={10≤t≤1,0elsewhere.
The analytic form of the B-spline is
(4)βK(t)=1(K-1)!∑k=0K(-1)k(Kk)(t-k)K-1u(t-k),
where the unit step function is defined as
(5)u(t)={1fort≥0,0fort<0.
An interesting property of the B-splines is that they obey the two-scale equation [3] as follows:
(6)βK(t2)=2∑kBK[k]βK(t-k),BK[k]=12K(Kk)=K!2Kk!(K-k)!,
where BK[k] is the binomial kernel.


## 2.2. Discrete B-Splines

The Laplace transform of the B-spline comes from(7)L{β(t)}=1s(1-e-s)⟹βK(s)=1sK(1-e-s)K.
The discrete B-spline βp[k] equals to the continuous B-spline at integer values of time. Hence, the Laplace transform (7) and the z-transform of the discrete B-spline should have inverse transforms which coincide at integer values in the time domain. Using the relation
(8)L-1(1sK)=tK-1(K-1)!u(t)
we obtain the z-transform of the discrete B-spline as follows:
(9)βK(z)=Z{βK[k]}=Z{L-1(1sK(1-e-s)K)}=NK(z)(1-z-1)K,
where
(10)NK(z)=Z{L-1(1sK)}=∑k=0∞kK-1(K-1)!z-k.
We have N1(z)=1/(1-z-1). By differentiating with respect to z, we obtain a recursion
(11)NK+1(z)=-zKdNK(z)dz.
Table 1 gives the NK(z) and the discrete B-splines βK(z) for K=1 to 6.Table 1
N
K
(
z
) and the discrete B-spline βK(z) for K=1 to 6.
1
1
1
-
z
-
1
1
2
z
-
1
(
1
-
z
-
1
)
2
z
-
1
3
z
-
1
+
z
-
2
2
(
1
-
z
-
1
)
3
z
-
1
+
z
-
2
2
4
z
-
1
+
4
z
-
2
+
z
-
3
6
(
1
-
z
-
1
)
4
z
-
1
+
4
z
-
2
+
z
-
3
6
5
z
-
1
+
11
z
-
2
+
11
z
-
3
+
z
-
4
24
(
1
-
z
-
1
)
4
z
-
1
+
11
z
-
2
+
11
z
-
3
+
z
-
4
24
6
z
-
1
+
26
z
-
2
+
66
z
-
3
+
26
z
-
4
+
z
-
5
120
(
1
-
z
-
1
)
4
z
-
1
+
26
z
-
2
+
66
z
-
3
+
26
z
-
4
+
z
-
5
120


## 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 form,n∈Z is generated from one-single function by dilation and translation operations as follows [9]
(12)ψm,n(t)=2-m/2ψ(2-mt-n)
which constitute an orthonormal basis. The wavelet transform of a signal x(t) is defined by the analysis/synthesis equations
(13)xnm=∫-∞+∞x(t)ψm,n(t)dt,x(t)=∑m∑nxnmψm,n(t),
where xnm are the wavelet coefficients in the subscale m. 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 ψm,n(t) behaves like the B-spline via the dilation and translation equation (12). If we keep the B-spline as a wavelet, the binomial kernel BK[n] 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)gm,b(t)=tme-btu(t)
is born by (m+1)-times convolution of the exponential e-btu(t) by itself. The term gamma spline originates from the gamma variate defined in the literature as tαe-t/β. For clarity, with the change of the variables to m=αandb=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
(15)gm,b(t-d)=(t-d)me-b(t-d)u(t-d)=∑k=0mak,dgk,b(t)u(t-d),
where
(16)ak,d=(mk)(-d)m-kebd.
The maximum value of the delayed gamma spline is obtained at time tmax=d+m/b and the area (1/b)m+1Γ(m+1). The continuous-time signal x(t) can be approximated by a weighted sum of delayed gamma splines as
(17)x(t)=∑kc[k]gm,b(t-k),
where c[k],k=0,1,2,… is the scale sequence. We may compute the scale sequence directly in time domain from the recursion
(18)x[1]=c[0]gm,b[1]⟹c[0]=x[1]gm,b[1],x[2]=c[0]gm,b[2]+c[1]gm,b[1]⟹c[1]=(x[2]-c[0]gm,b[2])gm,b[1],[2]⋮x[n]=∑k=0n-2c[k]gm,b[n-k]+c[n-1]gm,b[1]⟹c[n-1]=(x[n]-∑k=0n-2c[k]gm,b[n-k])gm,b[1].
In the following we introduce shortly the most essential gamma spline signal processing tools.Figure 1
Two typical gamma spline waveforms.Interpolation:
(19)x(t2)=2-m∑kc[k]gm,b/2(t-2k).Decimation:
(20)x(2t)=2m∑kc[k]gm,2b(t-k2).Derivative:
(21)ddtx(t)=∑kc[k](mgm-1,b(t-k)-bgm,b(t-k))=m∑kc[k]gm-1,b(t-k)-bx(t).Integral:
(22)∫0tx(t)dt=∑kc[k]∫0tgm,b(t-k)dt=b-1m∑kc[k]∫0tgm-1,b(t-k)dt-b-1x(t)
which follows from the integration by parts
(23)∫0tgm,b(t-k)dt=b-1m∫0tgm-1(t-k)dt-b-1gm,b(t-k).
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

Thez-transform of the gamma spline (14) may be conducted by differentiating
(24)Gm(z)=∑k=0∞kme-bkz-k⟹dGm(z)dz=-z-1∑k=0∞km+1e-bkz-k
which yields the recursion
(25)G0(z)=11-e-bz-1⟹Gm+1(z)=-zdGm(z)dz.
This yields a general form
(26)Gm(z)=Pm(z)(1-e-bz-1)m+1,
where Pm(z) is a polynomial in z-1, which is given in Table 2 for m=0to5. Gm(z) can be written in cascade realization as follows:
(27)Gm(z)=Pm(z)∏i=1m+111-e-bz-1.
The implementation of the gamma spline is possible as a discrete-time IIR filter in the case the pole e-b<1. If the pole lies outside the unit circle, the inverse filtering procedure (Appendix) has to be used.Table 2
P
m
(
z
) polynomial for m=0 to 5.
m
P
m
(
z
)
0
1
1
e
-
b
z
-
1
2
e
-
b
z
-
1
+
e
-
2
b
z
-
2
3
e
-
b
z
-
1
+
4
e
-
2
b
z
-
2
+
e
-
3
b
z
-
3
4
e
-
b
z
-
1
+
11
e
-
2
b
z
-
2
+
11
e
-
3
b
z
-
3
+
e
-
4
b
z
-
4
5
e
-
b
z
-
1
+
26
e
-
2
b
z
-
2
+
66
e
-
3
b
z
-
3
+
26
e
-
4
b
z
-
4
+
e
-
5
b
z
-
5


### 3.3. Scale Sequence

When we apply (17), for example, in interpolation or decimation we have to determine the scale sequence c[k]. In the following we present a feasible discrete-time filter for the computation of the scale coefficient sequence. The z-transform of (17) gives
(28)X(z)=C(z)Gm(z)⟹C(z)=Gm-1(z)X(z)
which yields
(29)C(z)=Pm-1(z)∏i=1m+1(1-e-bz-1)X(z).
The denominator polynomial Pm(z) 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 Pm-1(z) 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 thez-transform of the discrete signal x[n] is X(z), the interpolation of the discrete-time signal yields X(z2), which lacks the odd sequence. When we consider interpolation of the continuous signal x(t), the procedure yields also the odd (approximated intermediate) signal values. Based on (28) we have the z-transform of the interpolated signal as follows:
(30)Z{x[n2]}=C(z2)Z{gm,b[n2]},
where
(31)Z{gm,b[n2]}=∑k=0∞(k2)me-bk/2z-k=∑k=0∞kme-bkz-2k+z-1×∑k=0∞(k+12)me-b(k+1/2)z-2k.
Due to (24) we have
(32)Z{gm,b[n2]}=Gm(z2)+z-1Fm(z2),
where
(33)Fm(z2)=∑k=0∞(k+12)me-b(k+1/2)z-2k
which can be computed from the recursion
(34)F0(z)=e-b/21-e-bz-1⟹Fm+1(z)=12Fm(z)-zdFm(z)dz.
From (30) we obtain
(35)Z{x[n2]}=C(z2)[Gm(z2)+z-1Fm(z2)]=C(z2)Im(z),
where
(36)Im(z)=Gm(z2)+z-1Fm(z2).
Finally we have
(37)Z{x[n2]}=X(z2)+z-1C(z2)Fm(z2).
We may interpret the Im(z) as the gamma spline interpolation filter, which has the polyphase components Gm(z2) and Fm(z2). The application of Im(z) requires the computation of the scaling sequence C(z) via (29). We may observe that the polyphase component Fm(z2) 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 z transform of the shifted gamma spline comes from
(38)Z{gm,b[k+Δ]}=Gm,b(Δ,z)=∑k=0∞(k+Δ)me-b(k+Δ)z-k.
Using the binomial expansion
(39)(k+Δ)m=∑n=0m(mn)km-nΔn
we have
(40)Gm,b(Δ,z)=e-bΔ∑n=0m(mn)ΔnGm-n,b(z).
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)Z{x[k+Δ]}=C(z)Gm,b(Δ,z).
For example, the derivative of the signal can be calculated by differentiating with respect to Δ as follows:
(42)Z{ddtx[k+Δ]}=C(z)e-bΔ∑n=0m(mn)(nΔn-1-b)Gm-n,b(z).
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:
(43)Z{x[n-d]}=Z{x[n-1+1-d]}=z-1Z{x[n+1-d]},
where the delay d∈[0,1]. Let us now define Δ=1-d, which yields
(44)Z{x[n-d]}=z-1Z{x[n+Δ]}
which can be computed by applying the original time-shift algorithm (41).


## 3.1. Gamma Splines

In this work we introduce that the gamma spline representation(14)gm,b(t)=tme-btu(t)
is born by (m+1)-times convolution of the exponential e-btu(t) by itself. The term gamma spline originates from the gamma variate defined in the literature as tαe-t/β. For clarity, with the change of the variables to m=αandb=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
(15)gm,b(t-d)=(t-d)me-b(t-d)u(t-d)=∑k=0mak,dgk,b(t)u(t-d),
where
(16)ak,d=(mk)(-d)m-kebd.
The maximum value of the delayed gamma spline is obtained at time tmax=d+m/b and the area (1/b)m+1Γ(m+1). The continuous-time signal x(t) can be approximated by a weighted sum of delayed gamma splines as
(17)x(t)=∑kc[k]gm,b(t-k),
where c[k],k=0,1,2,… is the scale sequence. We may compute the scale sequence directly in time domain from the recursion
(18)x[1]=c[0]gm,b[1]⟹c[0]=x[1]gm,b[1],x[2]=c[0]gm,b[2]+c[1]gm,b[1]⟹c[1]=(x[2]-c[0]gm,b[2])gm,b[1],[2]⋮x[n]=∑k=0n-2c[k]gm,b[n-k]+c[n-1]gm,b[1]⟹c[n-1]=(x[n]-∑k=0n-2c[k]gm,b[n-k])gm,b[1].
In the following we introduce shortly the most essential gamma spline signal processing tools.Figure 1
Two typical gamma spline waveforms.Interpolation:
(19)x(t2)=2-m∑kc[k]gm,b/2(t-2k).Decimation:
(20)x(2t)=2m∑kc[k]gm,2b(t-k2).Derivative:
(21)ddtx(t)=∑kc[k](mgm-1,b(t-k)-bgm,b(t-k))=m∑kc[k]gm-1,b(t-k)-bx(t).Integral:
(22)∫0tx(t)dt=∑kc[k]∫0tgm,b(t-k)dt=b-1m∑kc[k]∫0tgm-1,b(t-k)dt-b-1x(t)
which follows from the integration by parts
(23)∫0tgm,b(t-k)dt=b-1m∫0tgm-1(t-k)dt-b-1gm,b(t-k).
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

Thez-transform of the gamma spline (14) may be conducted by differentiating
(24)Gm(z)=∑k=0∞kme-bkz-k⟹dGm(z)dz=-z-1∑k=0∞km+1e-bkz-k
which yields the recursion
(25)G0(z)=11-e-bz-1⟹Gm+1(z)=-zdGm(z)dz.
This yields a general form
(26)Gm(z)=Pm(z)(1-e-bz-1)m+1,
where Pm(z) is a polynomial in z-1, which is given in Table 2 for m=0to5. Gm(z) can be written in cascade realization as follows:
(27)Gm(z)=Pm(z)∏i=1m+111-e-bz-1.
The implementation of the gamma spline is possible as a discrete-time IIR filter in the case the pole e-b<1. If the pole lies outside the unit circle, the inverse filtering procedure (Appendix) has to be used.Table 2
P
m
(
z
) polynomial for m=0 to 5.
m
P
m
(
z
)
0
1
1
e
-
b
z
-
1
2
e
-
b
z
-
1
+
e
-
2
b
z
-
2
3
e
-
b
z
-
1
+
4
e
-
2
b
z
-
2
+
e
-
3
b
z
-
3
4
e
-
b
z
-
1
+
11
e
-
2
b
z
-
2
+
11
e
-
3
b
z
-
3
+
e
-
4
b
z
-
4
5
e
-
b
z
-
1
+
26
e
-
2
b
z
-
2
+
66
e
-
3
b
z
-
3
+
26
e
-
4
b
z
-
4
+
e
-
5
b
z
-
5


## 3.3. Scale Sequence

When we apply (17), for example, in interpolation or decimation we have to determine the scale sequence c[k]. In the following we present a feasible discrete-time filter for the computation of the scale coefficient sequence. The z-transform of (17) gives
(28)X(z)=C(z)Gm(z)⟹C(z)=Gm-1(z)X(z)
which yields
(29)C(z)=Pm-1(z)∏i=1m+1(1-e-bz-1)X(z).
The denominator polynomial Pm(z) 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 Pm-1(z) 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 thez-transform of the discrete signal x[n] is X(z), the interpolation of the discrete-time signal yields X(z2), which lacks the odd sequence. When we consider interpolation of the continuous signal x(t), the procedure yields also the odd (approximated intermediate) signal values. Based on (28) we have the z-transform of the interpolated signal as follows:
(30)Z{x[n2]}=C(z2)Z{gm,b[n2]},
where
(31)Z{gm,b[n2]}=∑k=0∞(k2)me-bk/2z-k=∑k=0∞kme-bkz-2k+z-1×∑k=0∞(k+12)me-b(k+1/2)z-2k.
Due to (24) we have
(32)Z{gm,b[n2]}=Gm(z2)+z-1Fm(z2),
where
(33)Fm(z2)=∑k=0∞(k+12)me-b(k+1/2)z-2k
which can be computed from the recursion
(34)F0(z)=e-b/21-e-bz-1⟹Fm+1(z)=12Fm(z)-zdFm(z)dz.
From (30) we obtain
(35)Z{x[n2]}=C(z2)[Gm(z2)+z-1Fm(z2)]=C(z2)Im(z),
where
(36)Im(z)=Gm(z2)+z-1Fm(z2).
Finally we have
(37)Z{x[n2]}=X(z2)+z-1C(z2)Fm(z2).
We may interpret the Im(z) as the gamma spline interpolation filter, which has the polyphase components Gm(z2) and Fm(z2). The application of Im(z) requires the computation of the scaling sequence C(z) via (29). We may observe that the polyphase component Fm(z2) 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 z transform of the shifted gamma spline comes from
(38)Z{gm,b[k+Δ]}=Gm,b(Δ,z)=∑k=0∞(k+Δ)me-b(k+Δ)z-k.
Using the binomial expansion
(39)(k+Δ)m=∑n=0m(mn)km-nΔn
we have
(40)Gm,b(Δ,z)=e-bΔ∑n=0m(mn)ΔnGm-n,b(z).
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)Z{x[k+Δ]}=C(z)Gm,b(Δ,z).
For example, the derivative of the signal can be calculated by differentiating with respect to Δ as follows:
(42)Z{ddtx[k+Δ]}=C(z)e-bΔ∑n=0m(mn)(nΔn-1-b)Gm-n,b(z).
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:
(43)Z{x[n-d]}=Z{x[n-1+1-d]}=z-1Z{x[n+1-d]},
where the delay d∈[0,1]. Let us now define Δ=1-d, which yields
(44)Z{x[n-d]}=z-1Z{x[n+Δ]}
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 Figure2. The analysis part consists of the low-pass scaling H0(z) and the high-pass wavelet filter H1(z). In the synthesis part the reconstruction filters G0(z) and G1(z) are related to the analysis filters by the perfect reconstruction (PR) condition as follows:
(45)[H0(z)H1(z)H0(-z)H1(-z)][G0(z)G1(z)]=[2z-k0].
The G0(z) and G1(z) are defined as
(46)G0(z)=H1(-z),G1(z)=-H0(-z).
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 sequenceswi[n],i=1,2,…,n and one scaling sequence sn[n].Lemma 1.
If the scaling filterH0(z), the wavelet filter H1(z), and the reconstruction filters G0(z)G1(z) are related according to the PR condition (45), the following modified analysis and synthesis filters obey the PR condition as follows:
(47)H0(z)=F(z)H0(z),H1(z)=F-1(-z)H1(z),G0(z)=F-1(z)G0(z),G1(z)=F(-z)G1(z),
where F(z) is any polynomial in z-1 and F-1(z) 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)H0(z)=Gm-1(z),H1(z)=z-1Gm(-z),G0(z)=z-1Gm(z),G1(z)=Gm-1(-z).
We may verify that by fixing H0(z)=1 and H1(z)=z-1 and following the result of the Lemma 1 the filter bank (48) obeys the PR condition (45). The gamma spline Gm(z) can be implemented by cascade of  m+1 IIR filters having exponential impulse response (27). The inverse filter
(49)Gm-1(z)=Pm-1(z)∏i=1m+1(1-e-bz-1)=Pm-1(z)∏i=1m+1R(z),
where Gm-1(z) is realizable by a cascade of m FIR filters having the impulse response R(z)=1-e-bz-1.
The implementation of the inverse filterPm-1(z) is described in the Appendix. The gamma spline scaling filter H0(z)=Gm-1(z) in (48) produces the scale sequence C(z) via (28). In multiscale analysis the scaling coefficients S(z)=C(z)↓2 (Figure 1) are fed to the following stage. The gamma spline wavelet H1(z)=z-1Gm(-z) 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 signalx[n] and the corresponding time shifted signal x[n-τ]; 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: h0[n] and h0[n-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ℋ(z), which has the frequency response
(50)ℋ(ω)=e-jπ/2sgn(ω),
where sgn(ω)=1 for ω≥0 and sgn(ω)=-1 for ω<0. The construction of ℋ(z) is based on the half-sample delay filter
(51)D(z)=z-0.5=A(z)B(z)
whose frequency response is D(ω)=e-jω/2. The quadrature mirror filter D(-z) has the frequency response D(ω-π)=e-j(ω-π)/2, correspondingly. Now the frequency response of the filter D(z)/D(-z) is
(52)D(ω)D(ω-π)=e-jω/2ej(ω-π)/2=e-jπ/2.
Comparing (50) and using the IIR filter notation (51) we obtain the Hilbert transform filter as
(53)ℋ(z)=A(z)B(-z)B(z)A(-z).
The half-delay filter D(z) in (51) can be realized by the polyphase components of the gamma spline interpolation filter Im(z), which yields
(54)D(z)=z-1Fm(z)Gm(z)⟹ℋ(z)=-Fm(z)Gm(-z)Gm(z)Fm(-z).


### 5.2. Shift Invariant Wavelet Transform

The Hilbert transform filter is inserted in the BF bank using the result of Lemma1 (47). The modified gamma spline wavelet transform filter bank is
(55)H-0(z)=ℋ(z)Gm-1(z),H-1(z)=ℋ-1(-z)z-1Gm(-z),G-0(z)=ℋ-1(z)z-1Gm(z),G-1(z)=ℋ(-z)Gm-1(-z).
The filter bank (33) can be highly simplified by noting the following equivalents:
(56)ℋ-1(-z)=ℋ(z),ℋ-1(z)=ℋ(-z).
By inserting (35) in (33) we obtain a highly simplified filter bank
(57)H-0(z)=ℋ(z)Gm-1(z),H-1(z)=ℋ(z)z-1Gm(-z),G-0(z)=ℋ(-z)z-1Gm(z),G-1(z)=ℋ(-z)Gm-1(-z).
The modified filter bank (57) can be realized by the Hilbert transform filter ℋ(z), which works as a prefilter for the analysed signal. On the other hand, the Hilbert transform filter ℋ(-z) 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(58)ℋa(z)=1+jℋ(z).
By filtering the real-valued signal x[n] by the Hilbert transform operator results in an analytic signal
(59)xa[n]=x[n]+jℋ{x[n]}
whose magnitude response is zero at negative side of the frequency spectrum
(60)Xa(ω)={2X(ω)0≤ω<π,0-π≤ω<0.
The wavelet sequence is obtained by decimation of the high-pass filtered analytic signal
(61)W(ω)=[Xa(ω)H1(ω)]↓2=Wa(ω)↓2=12Xa(ω2)H1(ω2).
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 Wa(ω) 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 w[n-τ/2] of the fractionally delayed signal x[n-τ] is e-jωτ/2Wa(ω/2)/2. The energy of the decimated wavelet coefficients is |W(ω/2)|/2, which does not depend on the fractional delay τ. Hence, the wavelet coefficients are shift invariant in respect to their energy content.


## 5.1. Hilbert Transform Filter

We define the Hilbert transform filterℋ(z), which has the frequency response
(50)ℋ(ω)=e-jπ/2sgn(ω),
where sgn(ω)=1 for ω≥0 and sgn(ω)=-1 for ω<0. The construction of ℋ(z) is based on the half-sample delay filter
(51)D(z)=z-0.5=A(z)B(z)
whose frequency response is D(ω)=e-jω/2. The quadrature mirror filter D(-z) has the frequency response D(ω-π)=e-j(ω-π)/2, correspondingly. Now the frequency response of the filter D(z)/D(-z) is
(52)D(ω)D(ω-π)=e-jω/2ej(ω-π)/2=e-jπ/2.
Comparing (50) and using the IIR filter notation (51) we obtain the Hilbert transform filter as
(53)ℋ(z)=A(z)B(-z)B(z)A(-z).
The half-delay filter D(z) in (51) can be realized by the polyphase components of the gamma spline interpolation filter Im(z), which yields
(54)D(z)=z-1Fm(z)Gm(z)⟹ℋ(z)=-Fm(z)Gm(-z)Gm(z)Fm(-z).


## 5.2. Shift Invariant Wavelet Transform

The Hilbert transform filter is inserted in the BF bank using the result of Lemma1 (47). The modified gamma spline wavelet transform filter bank is
(55)H-0(z)=ℋ(z)Gm-1(z),H-1(z)=ℋ-1(-z)z-1Gm(-z),G-0(z)=ℋ-1(z)z-1Gm(z),G-1(z)=ℋ(-z)Gm-1(-z).
The filter bank (33) can be highly simplified by noting the following equivalents:
(56)ℋ-1(-z)=ℋ(z),ℋ-1(z)=ℋ(-z).
By inserting (35) in (33) we obtain a highly simplified filter bank
(57)H-0(z)=ℋ(z)Gm-1(z),H-1(z)=ℋ(z)z-1Gm(-z),G-0(z)=ℋ(-z)z-1Gm(z),G-1(z)=ℋ(-z)Gm-1(-z).
The modified filter bank (57) can be realized by the Hilbert transform filter ℋ(z), which works as a prefilter for the analysed signal. On the other hand, the Hilbert transform filter ℋ(-z) 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(58)ℋa(z)=1+jℋ(z).
By filtering the real-valued signal x[n] by the Hilbert transform operator results in an analytic signal
(59)xa[n]=x[n]+jℋ{x[n]}
whose magnitude response is zero at negative side of the frequency spectrum
(60)Xa(ω)={2X(ω)0≤ω<π,0-π≤ω<0.
The wavelet sequence is obtained by decimation of the high-pass filtered analytic signal
(61)W(ω)=[Xa(ω)H1(ω)]↓2=Wa(ω)↓2=12Xa(ω2)H1(ω2).
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 Wa(ω) 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 w[n-τ/2] of the fractionally delayed signal x[n-τ] is e-jωτ/2Wa(ω/2)/2. The energy of the decimated wavelet coefficients is |W(ω/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 bym-times convolution of the exponential e-bt 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/(s+b) and 1/(1+e-bz-1). Depending on an implementation the application of the gamma spline parameters m and b can be optimally adjusted. In many applications we have selected m=b, which means that the maximum is attained at t=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 e-b 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 parameterm the shape of the gamma spline becomes more symmetrical, and the time support decreases. The m parameter affects the initial slope of the gamma spline and b the exponential damping rate. At very high values of m 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 m and b parameters. Hence, we may easily construct adaptive gamma spline wavelets, where m and b 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 m and b 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 N different excitation sequences are used and the excitation rate is f, the sampling rate can be reduced to f/N. 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.



---

*Source: 625364-2013-06-10.xml*