Parameter estimation of multi-component chirp signals based on discrete chirp Fourier transform and population Monte Carlo

Abstract

In this paper, a parameter estimation method for multi-component chirp signals with white Gaussian noise is proposed based on the modified discrete chirp Fourier transform (MDCFT) and the population Monte Carlo (PMC) methodology, in which the model order is unknown. By utilizing the integrability of linear parameters in the Bayesian model, this paper considers the posterior distribution of nonlinear parameters. MDCFT was adopted to calculate the chirpogram of the observed data, and clear peaks can be detected in the discrete chirp Fourier transform domain. The importance function (IF) was constructed according to the peaks, and the PMC algorithm was employed to evaluate the posterior distribution. The proposed method cannot only use the selected IF to generate the sample fitting target function in the parameter region of interest, but can also utilize samples and importance weights to update the IF adaptively. The simulation results indicated that the proposed method can realize joint Bayesian model selection and parameter estimation of multi-component chirp signals. Compared with the two existing methods based on Monte Carlo methodology, the proposed method exhibits improved performance.

Appendices

Appendix 1

The posterior distribution of $\left( {\varvec{\alpha }, {\varvec{f}},{\varvec{s}},\sigma _\varepsilon ^2 ,K} \right)$ can be expressed as follows according to the likelihood function of ${\varvec{y}}$ given in Eq. (3) and the a priori of each parameter given in Eqs. (6) and (7):

$$\begin{aligned}&p\left( {\varvec{\alpha } ,{\varvec{f}},{\varvec{s}},\sigma _\varepsilon ^2 ,K|{\varvec{y}}} \right) \!\propto \! p\left( {{\varvec{y}}|\varvec{\alpha } ,{\varvec{f}},{\varvec{s}},\sigma _\varepsilon ^2 ,K} \right) p\left( {\varvec{\alpha } |{\varvec{f}},{\varvec{s}},\sigma _\varepsilon ^2 ,K}\! \right) \nonumber \\&\quad \times p\left( {{\varvec{f}},{\varvec{s}}|K} \right) p\left( {\sigma _\varepsilon ^2 } \right) p\left( K \right) \nonumber \\&\quad \propto \frac{1}{\left( {\pi \sigma _\varepsilon ^2 } \right) ^{N}}\exp \left\{ {\!-\frac{1}{\sigma _\varepsilon ^2 }\left( {{\varvec{y}}\!-\!{\varvec{D}}\varvec{\alpha } } \right) ^{{\mathrm{H}}}\left( {{\varvec{y}}\!-\!{\varvec{D}}\varvec{\alpha } } \right) } \right\} \nonumber \\&\quad \times \frac{\left| {{\varvec{D}}^{{\mathrm{H}}}{\varvec{D}}} \right| }{\left( {\pi \delta ^{2}\sigma _\varepsilon ^2 } \right) }\exp \left\{ {-\frac{1}{\delta ^{2}\sigma _\varepsilon ^2 }\varvec{\alpha }^{{\mathrm{H}}}\left( {{\varvec{D}}^{{\mathrm{H}}}D} \right) \varvec{\alpha } } \right\} \nonumber \\&\quad \times \left( {\frac{1}{2}} \right) ^{K}\!\times \! \sigma _\varepsilon ^{2{\left( {-\upsilon _0 -1} \right) }}\exp \left[ {\frac{\!-\gamma _0 }{\sigma _\varepsilon ^2 }} \right] \!\times \! \frac{1}{K} \end{aligned}$$

(25)

Referring to the literature [12, 13, 18], the expression after integrating terms relevant to amplitudes $\varvec{\alpha }$ in Eq. (25) is given as

$$\begin{aligned}&p\left( {\varvec{\alpha } ,{\varvec{f}},{\varvec{s}},\sigma _\varepsilon ^2 ,K|{\varvec{y}}} \right) \nonumber \\&\quad \propto \frac{\left| {{\varvec{D}}^{{\mathrm{H}}}{\varvec{D}}} \right| }{\left( {\pi \delta ^{2}\sigma _\varepsilon ^2 } \right) }\exp \left\{ {-\frac{1}{\sigma _\varepsilon ^2 }\left( {\varvec{\alpha } -{\varvec{m}}_{\varvec{\alpha }} } \right) ^{{\mathrm{H}}}\varvec{{\varSigma }}_{\varvec{\alpha }}^{-1} \left( {\varvec{\alpha } -{\varvec{m}}_{\varvec{\alpha }} } \right) } \right\} \nonumber \\&\quad \times \frac{1}{\left( {\pi \sigma _\varepsilon ^2 } \right) ^{N}}\exp \left\{ {-\frac{1}{\sigma _\varepsilon ^2 }{\varvec{y}}^{{\mathrm{H}}}{\varvec{Py}}} \right\} \nonumber \\&\quad \times \sigma _\varepsilon ^{2\left( {-\upsilon _0 -1} \right) }\exp \left[ {\frac{-\gamma _0 }{\sigma _\varepsilon ^2 }} \right] \times \frac{1}{K\cdot 2^{K}} \end{aligned}$$

(26)

where

$$\begin{aligned}&\varvec{{\varSigma }}_{\varvec{\alpha }}^{-1} ={\varvec{D}}^{{\mathrm{H}}}{\varvec{D}}\left( {1+\delta ^{-2}} \right) \end{aligned}$$

(27)

$$\begin{aligned}&{\varvec{m}}_{\varvec{\alpha }} =\varvec{{\varSigma }}_{\varvec{\alpha }} {\varvec{D}}^{{\mathrm{H}}}{\varvec{y}} \end{aligned}$$

(28)

$$\begin{aligned}&{\varvec{P}}={\varvec{I}}_{N} -{{\varvec{D}}\left( {{\varvec{D}}^{{\mathrm{H}}}{\varvec{D}}} \right) ^{-1}{\varvec{D}}^{{\mathrm{H}}}}\big /{\left( {1+\delta ^{-2}} \right) } \end{aligned}$$

(29)

Equation (26) has a Gaussian function of $\varvec{\alpha }$ and inverse Gamma function of $\sigma _\varepsilon ^2$. After the integration of these two variables in (26), the posterior distribution of the observed data ${\varvec{y}}$ on $K, {\varvec{f}}$, and ${\varvec{s}}$ is

$$\begin{aligned} p\left( {{\varvec{f}},{\varvec{s}},K|{\varvec{y}}} \right) \propto \left( {\gamma _0 +{\varvec{y}}^{{\mathrm{H}}}{\varvec{Py}}} \right) ^{-\left( {N+\upsilon _0 } \right) }\times \frac{\left( {1+\delta ^{2}} \right) }{K\cdot 2^{K}} \end{aligned}$$

(30)

Appendix 2

We assume that $p\left( \varvec{\theta } \right)$ is the target function; $q^{\left( t \right) }\left( \varvec{\theta } \right)$ represents the IF in the $t$th iteration; and $\varvec{\theta }$ represents the unknown parameters. When applying the PMC methodology to the parameter estimation problem, whether the adopted IF can be approximated to the target function is important. Our objective is to minimize the Kullback divergence, which can be expressed as $K\left( {p|| q^{\left( t \right) }} \right) =\int {\log \left( {\frac{p\left( \varvec{\theta } \right) }{q^{\left( t \right) }\left( \varvec{\theta } \right) }} \right) p\left( \varvec{\theta } \right) \mathrm{d}\varvec{\theta }}$.

Assume $h_d \left( {\varvec{\theta };\alpha ,\varvec{\xi }} \right) =\frac{\alpha _d g\left( {\varvec{\theta };\varvec{\xi }_d } \right) }{\sum \nolimits _{d=1}^D {\alpha _d g\left( {\varvec{\theta };\varvec{\xi }_d } \right) } }$. In $t$th iteration, the intermediate quantity is constructed as

$$\begin{aligned}&L^{\left( t \right) }\left( {\alpha ,\varvec{\xi }} \right) \nonumber \\&\quad =\int {\sum \limits _{d=1}^D {h_d \left( {\varvec{\theta };\alpha ^{\left( t \right) },\varvec{\xi }_d^{\left( t \right) } } \right) \log \left( {\alpha _d^{\left( t \right) } g\left( {\varvec{\theta };\varvec{\xi }_d^{\left( t \right) } } \right) } \right) } p\left( \varvec{\theta } \right) \mathrm{d}\varvec{\theta }}\nonumber \\ \end{aligned}$$

(31)

The derivation in [26, 28, 29] indicates that for any $\alpha$ and $\varvec{\xi }$, when the intermediate quantity (31) increases, the target function also increases. The maximum $L^{\left( t \right) }\left( {\alpha ,\varvec{\xi }} \right)$ canobtain a closed solution. In the multivariate Gaussian distribution, the parameters include mean value and covariance matrix. The intermediate quantity can be expressed as

$$\begin{aligned}&L^{\left( t \right) }\left( {\alpha ,\varvec{\xi }} \right) \nonumber \\&\quad =\int \sum \limits _{d=1}^{D} h_d \left( {\varvec{\theta };\alpha ^{\left( t \right) },\varvec{\xi }_d^{\left( t \right) } } \right) \left\{ \log \left( {\alpha _d^{\left( t \right) } } \right) \right. \nonumber \\&\qquad \left. -\frac{1}{2}\left[ {\log \left| {\varvec{{\varSigma }}_d^{\left( t \right) } } \right| \!+\!\left( {\varvec{\theta }\!-\!\varvec{\mu }_d^{\left( t \right) } } \right) ^{\mathrm{T}}\left( {\varvec{{\varSigma }}_d^{\left( t \right) } } \right) ^{-1}\left( {\varvec{\theta }\!-\!\varvec{\mu }_d^{\left( t \right) } } \right) } \right] \right\} p\left( \varvec{\theta } \right) d\varvec{\theta }\nonumber \\ \end{aligned}$$

(32)

up to terms that do not depend on $\alpha ,\varvec{\xi }$.

When the above formula reaches the minimum, the following equations are satisfied:

$$\begin{aligned} \alpha _d^{\left( {t+1} \right) }&= \int {h_d \left( {\varvec{\theta };\alpha ^{\left( t \right) },\varvec{\xi }_d^{\left( t \right) } } \right) p\left( \varvec{\theta } \right) \mathrm{d}\varvec{\theta }} \end{aligned}$$

(33)

$$\begin{aligned} \mu _d^{\left( {t+1} \right) }&= \frac{\displaystyle \int {\varvec{\theta }^{\left( t \right) }h_d \left( {\varvec{\theta }^{\left( t \right) };\alpha ^{\left( t \right) },\varvec{\xi }_d^{\left( t \right) } } \right) p\left( \varvec{\theta } \right) \mathrm{d}\varvec{\theta }} }{\alpha _d^{t+1} } \end{aligned}$$

(34)

$$\begin{aligned}&\varvec{{\varSigma }}_d^{\left( {t+1} \right) }\nonumber \\&\quad =\frac{\displaystyle \int {\left( {\varvec{\theta }^{\left( t \right) }\!-\!\varvec{\mu }_d^{\left( {t\!+\!1} \right) } } \right) \left( {\varvec{\theta }^{\left( t \right) }\!-\!\varvec{\mu }_d^{\left( {t\!+\!1} \right) } } \right) ^{\mathrm{T}}h_d \left( {\varvec{\theta };\alpha ^{\left( t \right) },\varvec{\xi }_d^{\left( t \right) } } \right) p\left( \varvec{\theta } \right) \mathrm{d}\varvec{\theta }} }{\alpha _d^{t+1} }\nonumber \\ \end{aligned}$$

(35)

In practice, the numerator and the denominator are integral. By utilizing the samples and weights in each iteration for approximation, Eqs. (19)–(21) can be obtained.