The received Doppler signal of a stationary sensor, as emitted by a transiting acoustic source, is used to estimate source motion parameters, including speed, closest distance, rest frequency, and closest point of approach (CPA) time. First, the instantaneous frequency, amplitude, and CPA time are accurately estimated by the polynomial chirplet transform of the Doppler signal. Thereafter, the three other source motion parameters are obtained with a simplified amplitude-weighted nonlinear least squares method. The proposed scheme is successfully applied to the analysis of the characteristics of a moving truck.

The frequency of a pure tone emitted by a transiting source, when received by a stationary observer, changes with time because of the acoustical Doppler effect. This variation in instantaneous frequency (IF) can be used for source motion parameter estimation, which is a significant process in acoustics and has attracted much research attention.1–5 The basic principles involved are measuring the temporal variation in the nonlinear IF of the Doppler signal and then minimizing the sum of the squared deviations of the IF estimates from their predicted values. Therefore, central to the success of this frequency-based motion parameter estimation scheme is the need for an accurate IF estimate. Various IF estimation methods for the Doppler signal have been used. A common approach is to locate the peak of the time–frequency distribution (TFD), such as short-time Fourier transform (STFT), Wigner–Ville distribution (WVD), and polynomial WVD.6–9 However, these methods have some drawbacks, such as their low energy concentration on the IF law of the Doppler signal, appearance of inner artifacts, and intensive computation.

In this study, we consider the polynomial chirplet transform (PCT), which is a parameterized time–frequency analysis method developed by extension of the conventional chirplet transform.10 By replacing the chirplet kernel that has a linear IF law in the chirplet transform with a new kernel that has a polynomial nonlinear IF law, the PCT can generate a TFD that satisfies the energy concentration for the Doppler signal without the inner artifacts. The accurate estimation of the IF and its corresponding amplitude can be effectively achieved through a signal-dependent iterative procedure, the first iteration of which is the STFT. Meanwhile, the closest point of approach (CPA) time is directly estimated through a search of the minimum of the instantaneous frequency rate (IFR), which is obtained with the polynomial kernel coefficients of the PCT. As the polynomial kernel coefficients are estimated by the iterative procedure based on the STFT, the PCT is more robust than the polynomial phase estimation method that obtains the phase through Hilbert transform.6,11 It is also different from the multilinear function class methods, which are need to design a complicated signal kernel and are suitable for the polynomial phase signal of low order (usually less than 5).12–14 With the obtained IF, amplitude, and CPA time, the other motion parameters, such as the source speed, closest distance, and rest frequency, can be estimated with a simplified amplitude-weighted nonlinear least squares method. This proposed scheme is applied to the analysis of the characteristics of a moving truck. These features include the four motion parameters, the cylinder firing rate, engine firing rate, number of cylinders, and tachometer reading.

The geometry of the relative motion between a stationary sensor and a moving source is shown in Fig. 1. The source is traveling in a straight line at a constant speed v, so at time tc, the source is at the CPA to the stationary sensor at a separation distance R0. The source emits an acoustic tone of constant frequency f0 (i.e., the rest frequency), and the speed of sound in air is c. Because of Doppler effect, the IF of the signal received by the sensor at time t is given by2 

(1)

where α=f0c2/(c2v2), β=f0cv/(c2v2), y(t;tc,s)=(ttc)/s2+(ttc)2, and s=R0c2v2/cv. The IF of the Doppler signal is a function with respect to the four parameters { tc,v,R0,f0 }, which completely specify the source trajectory. Therefore, these motion parameters, or equivalently { tc,α,β,s }, can be estimated by measurement of the temporal variation in the IF and then the minimization of the sum of the squared deviations of the IF estimates from their predicted values over a period of time.

Fig. 1.

Geometry of a stationary sensor and a source traveling in a straight line at a constant speed.

Fig. 1.

Geometry of a stationary sensor and a source traveling in a straight line at a constant speed.

Close modal

As mentioned earlier, the IF law of the Doppler signal is a nonlinear function of time. The conventional chirplet transform, which is a parameterized time–frequency analysis method particularly designed for the analysis of chirp-like signals with the linear IF law, cannot track the evolution versus time of the IF of the Doppler signal. To improve the efficacy of the chirplet transform in analyzing signals with a nonlinear IF law, a modified version known as PCT is defined as15 

(2)

where z(t) is the analytical signal; hσ(t) stands for a nonnegative, symmetric, and normalized real window, usually taken as the Gaussian function e(t/σ)2/2/2πσ; and ΦaR(t) and ΦaM(t,t0) denote the nonlinear frequency rotating operator and the frequency shift operator, whose respective expressions are

(3)

Here, a=(a1,...,aK) are the polynomial kernel parameters.

Weierstrass approximation theorem16 states that any continuous function on a closed and bounded interval can be uniformly approximated by a polynomial of a finite-order K. Therefore, the Doppler signal can be expressed in the form of a polynomial phase signal, i.e., z(t)=A(t)exp(jk=1K+1bk1tk/k), where A(t) is the slowly varying amplitude, and (b1,,bK) denote the polynomial phase coefficients. Substituting the polynomial form of the Doppler signal into Eq. (2), we can obtain

(4)

When the kernel parameter vector α exactly matches the polynomial phase coefficients (b1,,bK), i.e., (a1,,aK)=(b1,,bK) of the Doppler signal, Eq. (4) is reduced to

(5)

where fr(t)=k=0Kbktk/2π is the polynomial form of the IF of the Doppler signal, hσ(f) is the Fourier transform of hσ(t), and f denotes the convolution operator with respect to f. Obviously, the frequency resolution of the PCT at all specific moments t0 is equal to the frequency bandwidth of the Gaussian window 1/σ. Therefore, the PCT can produce a high-quality TFD with an energy concentration on the IF law of the Doppler signal by giving a proper polynomial kernel parameter vector a, which can be determined by a signal-dependent iterative procedure with the termination condition

(6)

where η is a predetermined threshold, and IF(i)(t) denotes the estimated IF in each iteration. The iterative procedure for the IF estimates is described as follows:

  1. Initialization: Set η as a specific value. Set the polynomial order K. Set α=(a1,...,aK)=0, and set the window size.

  2. Estimation: Calculate the PCT with (a1,...,aK), and locate the peak in the obtained TFD of the PCT. Then, approximate the TFD peak data with a polynomial function of order K with the least square method and denote the coefficients as (a¯0,a¯1,...,a¯K).

  3. Recursion: Calculate the termination criterion ξ; if ξ>η, update α=(a1,...,aK) with (a¯1,...,a¯K) and repeat step (2); otherwise, go to step (4).

  4. Quit: Take the IF trajectory as the obtained K-order polynomial function of time.

The predetermination for the order K has been discussed in Ref. 14. Following the IF estimates, the peak values on the estimated IF law of the TFD are the estimated amplitude Â(t). Because the initialized kernel parameter vector α=(a1,...,aK)=0, the first iteration of the PCT is the STFT. The PCT can be viewed as a generalized STFT. Therefore, taking the Chirp z-transform17 instead of the Fourier transform in the PCT of the narrow-band Doppler signal can speed up the entire procedure.

The IFR is the derivative of the IF, so it can be also expressed in the following polynomial form:

(7)

The IFR has a minimum at the CPA time, i.e., when t=tc, IFR(tc)=IFRmin(t). Therefore, searching the minimum of the IFR estimation, which is obtained as k=1K1kaktk/2π, can facilitate the direct estimation of the CPA time.

The performance of the PCT-based IF and CPA time estimation method is assessed by a test on the numerically generated Doppler signal, whose IF has the following parameters: tc=15s, f0=118Hz, v=6m/s, R0=35m, and c=347m/s. The sampling frequency is 600 Hz and the duration of the Doppler signal is 30 s. The signal-to-noise ratio (SNR) is −5 dB, and it is defined as the ratio of average power of the Doppler signal and white Gaussian noise. When the PCT-based method is used to estimate the Doppler shift, the size of the Gaussian window is predetermined as 4500, the order of polynomial kernel K = 12, and the threshold δ=0.1%. For the purpose of comparison, Figs. 2(a)–2(c) depict the TFDs of the Doppler signal, which are generated by the STFT (first iteration of the PCT), PCT, and WVD. All three TFDs could render the inherent time-frequency pattern of the Doppler signal. However, the TFD generated by the PCT has the best concentration, the TFD by the STFT has a spreading concentration, and the TFD by the WVD has inner artifacts that disperse the concentration. Figure 2(d) shows the estimated IF curves of the three methods, as well as the real Doppler shift (denoted by red line). Corresponding to the respective TFDs, the IF estimated by the PCT (denoted by green line) overlaps the real Doppler shift, the IF by the STFT (denoted by blue line) cannot follow the nonlinear variation of the Doppler shift specially when the source approaches to the CPA, and the IF by the WVD (denoted by black line) is biased and sensitive to the noise. Figure 2(e) clearly shows that the IFR estimated by the PCT (denoted by green dot line) is almost the same as the real IFR (denoted by red line), and the estimated CPA time is 14.96 s that is approximately equal to the actual value 15 s. In conclusion, the PCT has been demonstrated to be effective for the analysis of Doppler signal, even when the SNR is low.

Fig. 2.

(Color online) The analysis of the numerical Doppler signal by using the PCT, STFT, and WVD: (a) TFD generated by the STFT, (b) TFD generated by the PCT, (c) TFD generated by the WVD, (d) the IF estimation curves and real IF, and (e) the IFR curve and CPA time estimated by the PCT.

Fig. 2.

(Color online) The analysis of the numerical Doppler signal by using the PCT, STFT, and WVD: (a) TFD generated by the STFT, (b) TFD generated by the PCT, (c) TFD generated by the WVD, (d) the IF estimation curves and real IF, and (e) the IFR curve and CPA time estimated by the PCT.

Close modal

After the IF estimates are obtained with the amplitude, the motion parameters { tc,v,R0,f0 }, or equivalently { tc,α,β,s }, can be estimated with a nonlinear least squares method, i.e., by minimization of the sum of the squared deviations of the IF estimates from their predicted values. The CPA time tc has been accurately estimated, so the amplitude-weighted NLS estimates of { v,R0,f0 } are specifically given by

(8)

where f̂r(tn) is the IF estimate at time tn, and N is the number of IF estimates. The weighted factor Â(tn) is used to highlight the accurate IF estimates in the duration when the transiting source approaches the CPA, and the SNR is relatively high. According to the linear form of the IF in Eq. (1), for a fixed s, Eq. (8) is solved as

(9)
(10)

where y(tn)y(tn;s), f̂¯r=[ n=1NÂ2(tn)f̂r(tn) ]/[ n=1NÂ2(tn) ], and y¯=[ n=1NÂ2(tn)y(tn) ]/[ n=1NÂ2(tn) ]. Substituting Eqs. (9) and (10) into Eq. (8) reduces the three-dimensional problem to a one-dimensional minimization problem

(11)

or equivalently, the parameter s is solved by minimizing the following cost function:

(12)

Substitution of the estimated value ŝ into Eqs. (9) and (10) obtains the corresponding parameters β̂ and α̂. Then, the other motion parameters can be calculated with the formulas v̂=β̂c/α̂, R̂0=ŝv̂c/c2v̂2, f̂0=α̂(c2v̂2)/c2. This approach uses an accurately estimated CPA time t̂c in advance, so it not only simplifies the solution procedure of the nonlinear least squares method, but it also improves the robustness and precision of the motion parameter estimation.

The proposed scheme is applied to the analysis of a radiated acoustic noise from a truck with six cylinders. The truck travels along a straight road with a constant speed of 20 km/h. It has a four-stroke diesel engine and a tachometer reading of 2350 revolutions per minute (rpm). The output of the 35 m far microphone is sampled at the rate of 12 000 Hz, and 30 s of data are recorded during transit of the truck.18 The speed of sound propagation in the air is a constant 347 m/s during data recording.

Figure 3 shows the results of the acoustic data processing. The spectrogram of the recorded signal in Fig. 3(a) shows there are a harmonic series with a fundamental frequency of about 20 Hz, which corresponds to cylinder firing rate. Among the harmonic series, the strongest line corresponding to the engine firing rate is observed at a frequency of about 120 Hz. Therefore, the number of cylinders is given as six through dividing engine firing rate by cylinder firing rate. This number is confirmed by visual inspection because every sixth harmonic of the cylinder firing rate is a strong line. Because of Doppler effect, the IF of the engine firing rate spectral line is observed to change with time. The parameters of PCT-based method are set as the same as the previous numerical simulation. Before the criterion reaches the predetermined threshold, three iterations of PCT have been implemented. Figure 3(b) shows the generated TFDs for the first (the top subgraph) and third iterations (the bottom subgraph). The concentration of the time–frequency representation of the third iteration is much better than that of the first, which is essentially produced by the STFT. Correspondingly, Fig. 3(c) shows that the IF estimated by the PCT (denoted by black circle) is more accurate than the IF by the STFT (denoted by blue line). The normalized amplitude estimates are presented in Fig. 3(d). Figure 3(e) depicts the IFR curve according to the estimated coefficients of the polynomial kernel, and the estimated CPA time t̂c=15.22s. With the use of the estimated IF and amplitude, the source speed, closest distance, and rest frequency (i.e., engine firing rate) are estimated as v̂=21.4km/hr, R̂0=35.3m, and f̂0=118.64Hz, respectively. Then, the estimated cylinder firing rate is known as 19.8 Hz, and the estimated tachometer reading is equal to 2 × cylinder firing rate × 60 = 2373 rpm.

Fig. 3.

(Color online) Data analysis of a moving truck with the proposed scheme: (a) spectrogram of the recording signal, (b) TFDs of the first and last iterations of the PCT-based IF estimate method, (c) IF curves obtained by the STFT and PCT, (d) the estimated amplitude, and (e) the CPA time estimated by the IFR.

Fig. 3.

(Color online) Data analysis of a moving truck with the proposed scheme: (a) spectrogram of the recording signal, (b) TFDs of the first and last iterations of the PCT-based IF estimate method, (c) IF curves obtained by the STFT and PCT, (d) the estimated amplitude, and (e) the CPA time estimated by the IFR.

Close modal

Table 1 shows the compared results of the PCT and two existing schemes, which use the STFT and WVD, respectively, to estimate the IF of the Doppler signal and then solve the nonlinear least squares problem by Gauss–Newton algorithm. All the three methods provide acceptable estimates for the characteristic parameters of the moving truck. However, the estimation results of the PCT are more accurate than those of the STFT and WVD.

Table 1.

Compared results of the source parameter estimation obtained by the proposed and existing schemes.

Actual valuePCTSTFTWVD
v (km/h) 20 21.4 22.1 21.9 
R0 (m) 35 35.3 38.2 37.7 
TRa (rpm) 2350 2373 2375 2375 
Actual valuePCTSTFTWVD
v (km/h) 20 21.4 22.1 21.9 
R0 (m) 35 35.3 38.2 37.7 
TRa (rpm) 2350 2373 2375 2375 
a

TR denotes the tachometer reading.

The acoustical Doppler shift is used to estimate the parameters of a moving source. Central to the success of this frequency-based passive acoustic source motion parameter estimation scheme is the need for an accurate IF estimate. We consider the PCT method, which produces a high-quality TFD, with the energy concentrated on the IF law of the Doppler signal. Together with the accurate IF estimate, the CPA time is estimated through minimization of the IFR, which is also obtained by the PCT. With the known CPA time, estimating the other parameters with a simplified amplitude-weighted nonlinear least squares method is easy. The proposed scheme is validated through an analysis of the characteristics of a moving truck.

This work was supported by the National Natural Science Foundation of China under Grant No.11274253.

1.
B. G.
Ferguson
and
B. G.
Quinn
, “
Application of the short-time Fourier transform and the Wigner-Ville distribution to the acoustic localization of aircraft
,”
J. Acoust. Soc. Am.
96
(
2
),
821
827
(
1994
).
2.
B. G.
Quinn
, “
Doppler speed and range estimation using frequency and amplitude estimates
,”
J. Acoust. Soc. Am.
98
(
5
),
2560
2566
(
1995
).
3.
D. C.
Reid
,
A. M.
Zoubir
, and
B.
Boashash
, “
Aircraft flight parameter estimation based on passive acoustic techniques using the polynomial Wigner-Ville distribution
,”
J. Acoust. Soc. Am.
102
(
1
),
207
223
(
1997
).
4.
H. X.
Zou
,
Y. Q.
Chen
,
J. H.
Zhu
,
Q. H.
Dai
,
G. Q.
Wu
, and
Y. D.
Li
, “
Steady-motion-based Dopplerlet transform: Application to the estimation of range and speed of a moving sound source
,”
IEEE J. Ocean. Eng.
29
(
3
),
887
905
(
2004
).
5.
K. W.
Lo
and
B. G.
Ferguson
, “
Flight parameter estimation using instantaneous frequency measurements from a wide aperture hydrophone array
,”
IEEE J. Ocean. Eng.
39
(
4
),
607
619
(
2014
).
6.
B.
Boashash
, “
Estimating and interpreting the instantaneous frequency of a signal—Part II: Algorithms and applications
,”
Proc. IEEE
80
(
4
),
540
568
(
1992
).
7.
V.
Katkovnik
and
L. J.
Stankovic
, “
Periodogram with varying and data-driven window length
,”
Signal Process.
67
(
3
),
345
358
(
1998
).
8.
B.
Barkat
and
B.
Boashash
, “
Design of higher order polynomial Wigner-Ville distributions
,”
IEEE Trans. Signal Process.
47
(
9
),
2608
2611
(
1999
).
9.
B.
Boashash
,
Time Frequency Signal Analysis and Processing: A Comprehensive Reference
(
Elsevier
,
London
,
2003
).
10.
Y.
Yang
,
Z. K.
Peng
,
X. J.
Dong
,
W. M.
Zhang
, and
G.
Meng
, “
General parameterized time-frequency transform
,”
IEEE Trans. Signal Process.
62
(
11
),
2751
2764
(
2014
).
11.
K. C.
Tam
,
S. K.
Tang
, and
S. K.
Lau
, “
On the recovery of moving source characteristics using time-frequency approach
,”
Appl. Acoust.
73
(
4
),
366
378
(
2012
).
12.
S.
Peleg
and
B.
Friedlander
, “
The discrete polynomial-phase transform
,”
IEEE Trans. Signal Process.
43
(
8
),
1901
1914
(
1995
).
13.
S.
Barbarossa
,
A.
Scaglione
, and
G. B.
Giannakis
, “
Product high-order ambiguity function for multicomponent polynomial-phase signal modeling
,”
IEEE Trans. Signal Process.
46
(
3
),
691
708
(
1998
).
14.
P.
O'Shea
and
R. A.
Wiltshire
, “
A new class of multilinear functions for polynomial phase signal analysis
,”
IEEE Trans. Signal Process.
57
(
6
),
2096
2109
(
2009
).
15.
Z. K.
Peng
,
G.
Meng
,
F. L.
Chu
,
Z. Q.
Lang
,
W. M.
Zhang
, and
Y.
Yang
, “
Polynomial chirplet transform with application to instantaneous frequency estimation
,”
IEEE Trans. Instrumen. Meas.
60
(
9
),
3222
3229
(
2011
).
16.
Available at http://mathworld.wolfram.com/WeierstrassApproximationTheorem.html (Last viewed November 15, 2014).
17.
L. R.
Rabiner
,
R. W.
Schafer
, and
C. M.
Rader
, “
The Chirp z-transform algorithm
,”
IEEE Trans. Audio Electroacoust.
17
(
2
),
86
92
(
1969
).
18.
Data available at http://acousticstoday.org/wp-content/uploads/2014/04/truck30sec.wav (Last viewed December 3, 2014).