AN IMPROVED ENERGY DEMODULATION ALGORITHM USING SPLINES

Dimitrios Dimitriadis and Petros Maragos

Department of Electrical & Computer Engineering,

National Technical University of Athens,
Zografou, 15773 Athens, Greece
ddim@cs.ntua.gr, maragos@cs.ntua.gr

ABSTRACT The instantaneous energy separation methodology has led to

several classes of algorithms for demodulating discrete-time AM-
A new algorithm is proposed for demodulating discrete-time FM signals
AM-FM signals, which first interpolates the signals with smooth
n
splines and then uses the continuous-time energy separation al-
gorithm (ESA) based on the Teager-Kaiser energy operator. This x[n] = x(nT ) = A[n] cos( Ω[k]dk) (3)
0
spline-based ESA retains the excellent time resolution of the ESAs
based on discrete energy operators but performs better in the pres- where A[n] = a(nT ) and Ω[n] = T ω(nT ). A direct approach is
ence of noise. Further, its dependence on smooth splines allows to apply the discrete-time Teager-Kaiser operator Ψd [xn ] ≡ x2n −
some optimal trade-off between data fitting versus smoothing. xn−1 xn+1 , where xn = x[n], to the discrete AM-FM signal (3)
and derive discrete energy equations of the form
1. INTRODUCTION Ψd [xn ] ≈ A2 [n] sin2 (Ω[n]) (4)

Demodulating AM-FM signals, i.e., nonstationary sines This yields the following algorithm, called Discrete ESA [3]:

 
t Ψd [xn −xn−1 ]+Ψd [xn+1 −xn ]
arccos 1− ≈ Ω[n] (5)
x(t) = a(t) cos( ω(τ )dτ ) (1) 4Ψd [xn ]

0

Ψd [xn ]
sin2 (Ω[n])
≈ |A[n]| (6)
that have a combined amplitude modulation (AM) and frequency
modulation (FM), has been a significant research problem with Another approach involves estimating the instantaneous frequency
many applications in communications systems, speech processing, by modeling the discrete-time signal xn via the exact Prony, as
and general nonstationary signal analysis. To solve it, a new ap- shown in [1, 5]. This yields algorithms that also contain the discrete
proach was developed in the 1990's based on nonlinear differential energy operator as their main ingredient. The best such algorithm
operators that can track the instantaneous energy or its derivatives is [1]  
of a source producing an oscillation [2, 4]. The main such rep- Υ3 [xn ] + Υ3 [xn−1 ]
resentative is the continuous-time Teager-Kaiser energy operator arccos ≈ Ω[n] (7)
4Ψd [xn ]
Ψ[x(t)] ≡ [ẋ(t)]2 − x(t)ẍ(t), where ẋ(t) = dx(t)/dt. Applied
to the AM-FM signal (1), Ψ yields the instantaneous source en- where Υ3 [xn ] ≡ xn xn+1 −xn−1 xn+2 is a third-order energy oper-
ergy, i.e. Ψ[x(t)] ≈ a2 (t)ω 2 (t), where the approximation error ator measuring the energy derivative [4]. The instantaneous ampli-
becomes negligible [3] if the instantaneous amplitude a(t) and in- tude is obtained from Eq. (4). We call the combination of Eq. (6),(7)
stantaneous frequency ω(t) do not vary too fast or too much with a Prony ESA.
respect to the average value of ω(t). Thus, AM-FM demodulation The advantages of the ESAs are efficiency, low computational
can be achieved by separating the instantaneous energy into its am- complexity and excellent time resolution (5-sample window) [3].
plitude and frequency components. Ψ is the main ingredient of the The main disadvantage is a moderate sensitivity to noise. A more
first energy separation algorithm (ESA) systematic approach is, hereby, developed where we first interpo-
late the discrete-time signal using smooth splines [6, 7], and then
 apply the continuous-time ESA of Eq. (2). We begin with a brief
Ψ[x(t)] Ψ[x(t)]
≈ ω(t) ,  ≈ |a(t)| (2) background on splines and then develop and test the new algorithm.
Ψ[ẋ(t)] Ψ[ẋ(t)]

developed by Maragos et al. [3] and used for signal and speech 2. SPLINE SIGNAL INTERPOLATION
AM–FM demodulation.
2.1. Exact Splines
This research work was supported by the Greek Secretariat for Research
and Technology and by the European Union under the program EΠET-98
Spline functions are piecewise polynomials constructed as a lin-
with Grant # 98ΓT26 and the program ΠENE∆-99 with Grant # 99E∆164. ear combination of B-Splines. A spline function of order n has
It was also supported by the basic research program ARCHIMEDES of the continuous (smooth) derivatives up to order n − 1, which is very
NTUA Institute of Communication and Computer Systems. important when using Ψ.
 Given the initial signal samples x[n], n = 1, . . . , N , the inter- curve. The positive parameter λ quantifies how smooth the interpo-
polating spline function is given by lating curve will be and how close to the data points the interpolant
will pass. For λ = 0 there is no smoothing and the interpolation

+∞
curve fits exactly the signal samples. If λ = 0, the deviation from
gν (t) = c[n]βν (t − n) (8) the data samples increases with the parameter λ.
n=−∞ As shown in [6, 7], the interpolating function sν (t) minimizing
the mean square error  is a linear combination of splines βν , as in
where βν (t) is the B-spline of order ν, and the coefficients c[n] Eq. (8), but the coefficients c[n] are computed as the output of an
depend only on the data x[n] and the analytic expression of the IIR filter that is different from the filter in (10):
B-spline. The B-spline can be formed as the the (ν + 1) th–fold
convolution of the zeroth–order B-spline with itself: X(z)
C(z) = Hνλ (z)X(z) = (11)
Pνλ (z)
βν (t) ≡ β0 (t) ∗ β0 (t) ∗ · · · ∗ β0 (t),


(ν+1) times where Pνλ (z) is equal to
ν+1
where the zeroth–order B-spline is defined by Pνλ (z) = Bν (z) + λ(−z + 2 − z −1 ) 2 (12)

1 if −1/2 < t < 1/2 The IIR filter Hν (z) has a symmetric impulse response and all
β0 (t) = 1/2 if |t| = 1/2 its poles are inside the unit circle. Thus, the spline coefficients
0 otherwise c[n] can be determined stably via a few recursive equations [6, 7].
Henceforth, smooth splines with λ = 0 are applied unless stated
Using the discrete B-spline bν (n) ≡ βν (n), Eq. (8) becomes otherwise.
gν (n) = (c ∗ bν )(n) (9)
3. SPLINE ESA
Note that for the exact interpolation problem is gν (n) = x[n]. By
transforming Eq. (9) in the Z domain, Our previous discussion has established that, by using the theory of
smooth splines, we can interpolate the original discrete-time signal
X(z) x[n] using smooth splines of νth order and create a continuous-time
C(z) = (10)
Bν (z) signal

+∞
Thus, the spline coefficients c[n] can be determined recursively sν (t) = c[n]βν (t − n) (13)
from the above equation. Note that each original sample gν (n) = n=−∞
x[n] is resynthesized by the contributions of ν + 1 neighbor spline
coefficients. Obviously, the signal sν (t) is a continuous-time expansion of the
original discrete signal x[n]. Now, the basic idea of the new ap-
proach for ESA-based demodulation is to apply the continuous-time
2.2. Smooth Splines energy operator Ψ and the continuous ESA to the continuous-time
We have used splines to improve the performance of the ESA algo- signal sν (t), instead of using the discrete energy operator Ψd and
rithm. Several experiments were conducted concerning different the DESA on the discrete signal x[n]:
kinds of splines. In these experiments we used noisy AM-FM sig-  2
nals with different levels of SNR. The results were disappointing ∂sν (t) ∂ 2 sν (t)
Ψ[sν (t)] = − sν (t) (14)
as the exact fitting of the curve, due to the presence of noise, was ∂t ∂t2
creating large estimation errors. This led to the search for a solution
dealing with the problem of noise. The next step was to find a way In order to use the continuous ESA we have to compute the first-,
to optimally interpolate signal samples passing closely but not ex- second- and third-order derivatives of the signal. Given the coef-
actly through them. Smooth splines are ideal for this purpose. The ficients c[n] of the spline interpolation (13), we can derive the fol-
main advantage of smooth splines is that the interpolating polyno- lowing closed-form expressions for these derivatives that involve
mial does not pass through the signal samples but close enough. only the coefficients c[n] and the B-spline functions:
This concept is similar to the idea of polynomial fitting with least
squares. ∂sν (t)
= (c[n] − c[n − 1]) βν−1 (t − n + 1/2) (15)
The smooth spline interpolating function is the function sν of ∂t
n
order ν = 2r − 1 that minimizes the mean square error criterion
∂ 2 sν (t)

+∞
+∞  r
2 = (c[n + 1] − 2c[n] + c[n − 1]) βν−2 (t − n)
∂ sν (x) ∂t2
= (x[n] − sν (n))2 +λ dx n
−∞
∂xr (16)
n=−∞


 ∂ 3 sν (t)

d

s = (c[n + 1] − 3c[n] + 3c[n − 1] − c[n − 2])
∂t3
n
where d is the mean square error of the interpolation function and ·βν−3 (t − n + 1/2)
s is the mean square error introduced by the need for a smoothed (17)
curve. This criterion is a compromise between the need for close- By using these signal derivatives in the continuous ESA equations
to-the-data points interpolation curve and the need for a smoothed (2), we can estimate the instantaneous amplitude a(t) and frequency
ω(t) of the continuous signal sν (t). Finally, by sampling these 4. EXPERIMENTS AND DISCUSSION
information-bearing signals, we obtain estimates of the instanta-
neous amplitude A[n] = a(nT ) and frequency Ω[n] = T ω(nT ) For the experiments conducted we used parametric AM-FM signals,
of the original discrete signal x[n]. This whole approach above is [1, 3] and different levels of white Gaussian noise:
called the Spline ESA. πn πn πn
An important part of the Spline ESA is the computation of the x[n] = (1 + 0.05k cos( )) cos( + m sin( )) + e[n],
100 5 100
spline coefficients c[n]. Next we discuss the details of this pro-
cedure. First, we find the zeros of the denominator polynomial where e[n] is the added noise and m = 1, . . . , 10, k = 1, . . . , 10.
Pνλ (z) in Eq. (11). Due to the symmetric form of this polynomial, The first experiments concerned Spline ESA using exact splines,
the zeros come in pairs (zi , zi−1 ), i = 1, . . . , r. Thus, the transfer from fifth to seventh order. The experimental results were not so
function in Eq. (11) can be written as promising as the frequency mean estimation error of the spline ESA
was always worse than the corresponding one of the Prony ESA

r
−zi and the DESA as shown in Fig. 1.
Hν (z) = c0 (18)
(1 − zi z −1 )(1 − zi z)
i=1 10
Prony ESA
DESA
From Eq. (11), (18) and [6, 7] the recursive equations are:

Frequency Estimation Mean Abs. Error %

9 Exact Splines, Order=5
Exact Splines, Order=6
Exact Splines, Order=7
yi+ [n] = yi−1 [n] + zi y + [n − 1] , n = 2, . . . , N 8

yi [N ] = ai (2yi+ [N ] − yi−1 [N ]) (19) 7

6
yi [n] = zi (yi [n + 1] − yi+ [n]) , n = N − 1, . . . , 1
5

where ai = −zi /(1 − zi2 ), yi−1 [n] is the input and yi [n] is the
4
output of a digital filter with transfer function
3
−zi
Ti (z) = ,
(1 − zi z −1 )(1 − zi z) 2

and y0 [n] = x[n]. If the above step is repeated as many times as 1

the number of pole pairs (zi , zi−1 ), the final output sequence yr [n]
0
will equal c[n]. As boundary condition we set, 40 45 50 55 60 65 70 75 80
SNR (dB)

k0
|k−1|
yi+ [1] = zi yi−1 [k] Fig. 1. Comparison of Spline ESA with exact splines with DESA
i=1 and Prony ESA.

where k0 is an integer that ensures a certain level of precision. The exact spline fitting is responsible for the large mean square
We present an example to clarify some of the steps in determining error of the algorithm, as the noisy samples insert a significant er-
c[n] for splines of order ν = 5(r = 3) and for two different values ror. Since noisy signals have large variations and sharp edges, the
of λ, λ = 0 and λ = 0.5. First, if λ = 0, we have interpolation need of a smoothing factor is apparent. So, the next series of exper-
with exact splines of order ν = 5. The denominator of the transfer iments was conducted using smooth splines with different values of
function H5 (z) will be P5 (z) ≡ B5 (z) and the poles will be, λ through a wide range of SNR values, which are shown in Fig. 2.
Spline ESA using smooth splines performs better than the DESA
z1 = −0.04309, z2 = 0.43057,
and the Prony ESA in the presence of noise. The choise of value of
z3 = z1−1 , z4 = z2−1
λ is not completely arbitrary. We have attempted to find experimen-
Now, for λ = 0.5, by setting r = 3 in Eq. (12) we find the poles tally a good value for λ for different SNRs. Fig. 3 shows the mean
of H5 (z) = 1/P5 (z): error of the cubic and fifth–order smooth splines for various values
of λ when SNR=35 dB. In these experiments the corresponding
z1 = 0.32548, mean square error curves had minima for particular values of λ.
z2 = 0.32154 − 0.47128i, z3 = 0.32154 + 0.47128i, More specifically, the minima occurred when λ ∈ [0.1, 1] inde-
z4 = z1−1 , z5 = z2−1 , z6 = z3−1 pendently of the SNR's values. Note that, the mean square error
of the fifth–order smooth spline is always smaller that the corre-
In both cases we find c[n] by using the algorithm of Eq. (19). The sponding one of the cubic smooth spline independently of the value
only difference between the case with λ = 0 and the case with of λ or SNR. The optimal value of λ is not known and can be de-
λ = 0 is the number and the values of the poles. Finally, having termined only through experimentation. This happens because the
computed c[n], the coefficient signal is convolved with the B-spline errors depend on the SNR, the signal, and the application.
βν to yield the inerpolating signal sν . Now, we compare the Spline ESA with DESA and Prony ESA,
In general, the evaluation of the spline coefficients by the fil- fixing λ equal to a constant value. For SNR= ∞, λ could be
tering approach presented above is less computationally complex set to zero. Otherwise for SNR= ∞, λ takes values in the order
than the standard numerical analysis approach using sparse Toeplitz of 0.25. The order of the smooth spline will be ν = 5 (r = 3)
matrices. because the mean square error of the Spline ESA with fifth-order
 25 14
Prony ESA Smooth Splines, Order=3
DESA Smooth Splines, Order=5

Frequency Estimation Mean Abs. Error %

Frequency Estimation Mean Abs. Error %

Smooth Splines, Order=3 SNR=35 db

Smooth Splines, Order=5 12

20

10

15
8

6
10

5
2

0 0
30 35 40 45 50 55 60 65 70 0.0075 0.01 0.025 0.05 0.075 0.1 0.25 0.5 0.75 1 2.5 5
SNR (dB) Values of λ

Fig. 2. Comparison of Spline ESA with smooth splines (λ = 0.25) Fig. 3. Mean Square Frequency Error for Spline ESA with cubic
with DESA and Prony ESA. and fifth-order smooth splines when SNR=35 dB.

smooth splines is always smaller than the error of the Spline ESA In conclusion the proposed new algorithm provides an efficient
with cubic splines, as it is shown in Fig. 3. In addition, by setting way to represent a discrete signal in continuous-time domain, which
ν = 5, the time–window (i.e., the number of coefficients required is very important for applying the Teager-Kaiser energy operator
to produce one output sample) of Spline ESA is equal to that of the where the signal's derivatives must be determined with high preci-
DESA and the Prony ESA. sion. Further, using smooth splines gives the new algorithm some
additional robustness in the presence of noise.
We tested all three algorithms for SNR= +∞ (no noise is added
to the signal) and for SNR=45 dB. In Table 1 we compare the DESA,
Prony ESA and Spline ESA with different values of λ. When 5. REFERENCES
the SNR is small, smooth splines have better performance and the
estimation error is much smaller than the corresponding one for [1] L. B. Fertig and J. H. McClellan, "Instantaneous Frequency
the exact splines. This is the main reason for using smooth splines Estimation Using Linear Prediction With Comparisons to the
instead of exact ones. The Prony ESA algorithm, even though it DESAs", IEEE Signal Process. Lett., vol. 3, pp. 54-56, Feb.
yields very good results (mean square error 0.14%) when there is 1996.
no noise, is not robust for noisy signals and the corresponding error [2] J. F. Kaiser, "On a simple algorithm to calculate the ènergy' of
increases dramatically as the signal's SNR is decreasing. On the a signal", in Proc. IEEE ICASSP-90, Albuquerque, NM, Apr.
contrary, the DESA and especially, the Spline ESA, when λ = 0, 1990, pp. 381–384.
are more robust for signals with low SNRs. This fact is clearly
[3] P. Maragos, J. F. Kaiser, and T. F. Quatieri, "Energy Separation
represented in Fig. 2 and Table 1. Also, it must be noticed that in
in Signal Modulations with Application to Speech Analysis",
Fig. 2 the value of λ does not change with the SNR values (its value
IEEE Trans. Signal Processing, vol. 41, pp. 3024–3051, Oct.
is constant). This is why for SNR values bigger than 70 dB, DESA
1993.
seems to perform better than Spline ESA. Changing λ in inverse
proportion to the SNR value, the Spline ESA yields a mean square [4] P. Maragos and A. Potamianos, "Higher Order Differential En-
error that is smaller than the corresponding one of DESA. ergy Operators", IEEE Signal Process. Lett., vol. 2, p. 152-154,
Aug. 1995.
[5] C. S. Ramalingam, "On the Equivalence of DESA-1a and
Prony's Method When the Signal is a Sinusoid", IEEE Signal
Table 1. Comparison of ESA Demodulation Algorithms. Process. Lett., vol. 3, pp. 141-143, May 1996.
Frequency Mean Absolute Error (%) [6] M. Unser, "Splines A Perfect Fit for Signal and Image Process-
SNR= +∞ dB ing", IEEE Signal Process. Magaz., pp. 22-38, Nov. 1999.
DESA Prony ESA Spline ESA Spline ESA [7] M.Unser, A. Aldroubi and M. Eden "B-Spline Signal Process-
ν = 5, λ = 0 ν = 5, λ = 0.25 ing: Part I–Theory. Part II–Efficient Design and Applications",
0.38767 0.14774 0.37205 0.79749 IEEE Trans. Signal Processing, vol. 41, pp. 821–848, Feb.
SNR=45 dB 1993.
DESA Prony ESA Spline ESA Spline ESA
ν = 5, λ = 0 ν = 5, λ = 0.25
1.32737 6.02938 4.56928 0.95701