Recovery of the starting times of delayed signals

TL;DR

This paper introduces a novel iterative method for accurately identifying the starting times of multiple delayed signals within finite data sequences, even under noisy conditions.

Contribution

The paper proposes a new iterative approach to locate the starting points of multiple delayed signals, including damped ones, in finite data sequences, with demonstrated noise robustness.

Findings

Effective in locating the longest delay first

Iterative method successfully identifies all starting points

Noise sensitivity tests show robustness under weak noise

Abstract

We present a new method to locate the starting points in time of an arbitrary number of (damped) delayed signals. For a finite data sequence, the method permits to first locate the starting point of the component with the longest delay, and then --by iteration-- all the preceding ones. Numerical examples are given and noise sensitivity is tested for weak noise.

Equations

$$G(z)=\sum _{n\ge 0}{s}_n{z}^n.$$

$$U(z)=\frac{G^{\prime }(z)}{G(z)}.$$

$$\displaystyle s_{n}=\left\{\begin{array}[]{ll}0&0\leq n<n_{k}\\ Ae^{i\omega\frac{T}{N}(n-n_{k})}\quad\omega=2\pi\nu+i\alpha,\quad n\geq n_{k};\end{array}\right.$$

$$G_k(z)=z^{n_k}\frac{A_k}{1-z{e}^{i{\omega }_k\frac{T}{N}}},$$

$$G_{signal}(z)={\Sigma }_{k=1}^K{G}_k(z)=\frac{{\Sigma }_{k=1}^K{z}^{n_k}{A}_k{\Pi }_{j=1,K}^{j\ne k}{Q}_j^{(1)}(z)}{{\Pi }_{j=1,K}{Q}_j^{(1)}(z)}=\frac{{\Sigma }_{k=1}^K{z}^{n_k}{Q}_k^{(K-1)}(z)}{Q_{signal}^{(K)}(z)}=\frac{P_{signal}^{(K-1+n_K)}(z)}{Q_{signal}^{(K)}(z)},$$

$$G_{total}(z)=\frac{P_{signal}^{(K-1+n_K)}(z)}{Q_{signal}^{(K)}(z)}+\frac{P_{noise}^{(H)}(z)}{Q_{noise}^{(H)}(z)}=\frac{P_{total}^{(K+H-1+n_K)}(z)}{Q_{total}^{(K+H)}(z)}=\frac{{\Pi }_{i=1}^{K+H-1+n_K}(z-z_i)}{{\Pi }_{j=1}^{K+H}(z-p_j)},$$

$$U(z)={\Sigma }_{j=1}^{K+H}\frac{-1}{z-p_j}+{\Sigma }_{1=1}^{K+H-1+n_K}\frac{+1}{z-z_i}.$$

$$s_n=A{e}^{i\omega \frac{T}{N}|n-n_k|}\omega =2\pi \nu +i\alpha ,\alpha >0;$$

$$G_k^{sym}(z)=A_k{e}^{i{\omega }_k\frac{T}{N}{n}_k}{\Sigma }_{n=0}^{n_k-1}{\left(z{e}^{-i{\omega }_k\frac{T}{N}}\right)}^n+A_k{z}^{n_k}{\Sigma }_{n=0}^{\infty }{\left(z{e}^{i{\omega }_k\frac{T}{N}}\right)}^n.$$

$$G_k^{sym}(z)=P^{(n_k-1)}(z)+z^{n_k}\frac{A_k}{1-z{e}^{i{\omega }_k\frac{T}{N}}}=\frac{Q_k^{(1)}(z)P^{(n_k-1)}(z)+A_k{z}^{n_k}}{Q_k^{(1)}(z)},$$

$$U(z)=\frac{{\Sigma }_{n=0}^{\infty }{b}_n{z}^n}{{\Sigma }_{n=0}^{\infty }{a}_n{z}^n}$$

$$c_n=b_n-\frac{1}{a_0}{\Sigma }_{k=1}^n{c}_{n-k}{a}_k$$

$${\gamma }_n=(n+1){\alpha }_{n+1}-{\Sigma }_{k=1}^n{\gamma }_{n-k}{\alpha }_k.$$

$$n_K<\frac{N-3}{2}$$

$$\displaystyle s_{n}=\left\{\begin{array}[]{ll}\cos(2*\pi*f_{1}*n)&0\leq n<n_{1}\\ \cos(2*\pi*f_{1}*n)+e^{-i\gamma(n-n_{1})}\cos(2*\pi*f_{2}*(n-n_{1}))&n\geq n_{1}\end{array}\right.$$

$$\displaystyle s_{n}=\left\{\begin{array}[]{ll}e^{i\omega_{1}n}+0.9e^{i\omega_{2}n}&0\leq n<n_{1}\\ e^{i\omega_{1}n}+0.9e^{i\omega_{2}n}+e^{-\gamma(n-n_{1})}e^{\omega_{3}*(n-n_{1})}&n\geq n_{1}.\end{array}\right.$$

$$\displaystyle s_{n}=\left\{\begin{array}[]{ll}e^{i\omega_{1}n}+0.9e^{i\omega_{2}n}&0\leq n<n_{1}\\ e^{i\omega_{1}n}+0.9e^{i\omega_{2}n}+e^{-\gamma(n-n_{1})}e^{\omega_{3}*(n-n_{1})}&n_{1}\leq n<n_{2}\\ e^{i\omega_{1}n}+0.9e^{i\omega_{2}n}+e^{-\gamma_{3}(n-n_{1})}e^{\omega_{3}*(n-n_{1})}+e^{-\gamma_{4}(n-n_{2})}e^{\omega_{4}*(n-n_{2})}&n\geq n_{2}\end{array}\right.$$

$$\displaystyle s_{n}=\left\{\begin{array}[]{ll}e^{i\omega_{1}n}+0.9e^{i\omega_{2}n}+e^{-\gamma(n_{1}-n)}e^{\omega_{3}*(n_{1}-n)}&0\leq n\leq n_{1}\\ e^{i\omega_{1}n}+0.9e^{i\omega_{2}n}+e^{-\gamma(n-n_{1})}e^{\omega_{3}*(n-n_{1})}&n>n_{1}.\end{array}\right.$$

Full text

I Introduction

In signal processing we often encounter situations where the onset of the signal –or of part of it– is delayed. Apart from cases when the starting point itself is a relevant parameter to be measured (for example in active radar and sonar detection systems), this delay often results in degradation of the other measured signal parameters, especially the amplitudes of the various frequency components. It is therefore convenient to locate the signal starting point so as to be able to optimize the analysis of the data sequence being studied.

With the exception of wavelets (or multiscale) techniques, which are known to be efficient at locating edges wav1 ; wav3 and some parametric models which take into account possible delays boy , most methods in the literature compare some measured quantity for a series of contiguous traveling windows, be it power timfre1 ; timfre2 , correlation to a library of known waveforms burst , or some other quantities mar .

The method we propose instead considers a single window; moreover –when different signal components start at different times– it locates first the onset time for the signal frequency component with the longest delay, i.e. the first data point where all the signal components are present. This is often the sought for information. By iteration it can then locate the starting points for all other components, should they be needed too.

Our method is based on the analytic properties of the generating function of the data series, which –as we shall see– depend on the signal delay. More in detail: given an infinite data series $s_{0},s_{1},s_{2},\cdots,s_{n},\cdots$, we can build its generating function, or “G-transform”:

[TABLE]

In practice, we always have data series that are truncated at order $N-1$, therefore consisting of $N<\infty$ data points. Assuming the total signal to be the discrete sum of a finite number of (damped) harmonic oscillations, we want to estimate the G-transform of the infinite time series by a suitable Padé Approximant: a rational function whose Taylor expansion around $z=0$ reproduces the G-transform up to order $N-1$. The problem is that any ratio $[I/L]$ of the order $I$ of the numerator to the order $L$ of the denominator such that $I+L=N-1$ is a candidate Padé Approximant; which is the best one?

We know that if all the signals are already present at $n=0$, then –even in the presence of noise– the best choice is a rational function of type $[M-1/M]$ noi , requiring $N=2M$ data points. From its parameters all the (complex) amplitudes and frequencies can be deducted, provided that $M\geq K$, where $K$ is the total number of frequencies in the signal, including complex conjugate pairs.

Things change if some of the frequency components of the signal start at some $n_{k}$’s larger than zero:

let’s assume we have $K$ signals starting at $n_{1},n_{2},\cdots,n_{K}$ with $0\leq n_{1}\leq n_{2}\leq\cdots\leq n_{K}$; as we shall see, the G-transform is then a rational function of the type $[M-1+n_{K}/M]$.

Our first aim is therefore to find the (discrete) time $n_{K}$ at which the last frequency component starts. To do this, we introduce the U-transform which is the logarithmic derivative of the G-transform:

[TABLE]

The poles of the U-transform are the poles and zeros of the G-transform and their residues are $-1$ for the G-transform poles and $+1$ for the zeros; the sum of the residues of the U-transform therefore is $n_{K}-1$.

If our aim is to find the best Padé Approximant for the given sample, we are then done, as $n_{K}$ is the sought for information.

If instead we want to find all the starting points $n_{k}$, $k=1,\cdots K$, we just have to repeat the procedure with a shorter sample of length $2n_{K}$ to find $n_{(K-1)}$, and so on, until the sample length is reduced to zero.

II The theoretical background of the method

Let us start considering a single discretized signal starting at an unknown time $n_{k}$:

[TABLE]

the signal can be transient (i.e. damped: $\alpha>0$), of constant amplitude ($\alpha=0$), or even have zero frequency (corresponding to a step function).

The G-transform of this signal reads:

[TABLE]

which is a $[n_{k}/1]$ rational function $z^{n_{k}}/Q_{k}^{(1)}(z)$, where $Q_{k}^{(j)}(z)$ is a polynomial in $z$ of order $j$ and the index indicates that it refers to the signal $k$. Summing $K$ such signals with different starting times, we get

[TABLE]

where $n_{K}=\max_{k=1,K}(n_{k})$ and we have indicated with $P$ the polynomial at the numerator and with $Q$ the one at the denominator.

Adding noise does not change the the substance of the above picture, as it means to add an equal number $H$ of poles and zeros (Froissart doublets dou1 ; dou2 ; dou3 ; dou4 ) so that the complete G-transform reads:

[TABLE]

where we have respectively called $z_{i}$ and $p_{j}$ the zeros and poles of $G_{total}(z)$, which is a $[(K+H-1+n_{K})/(K+H)]$ rational function, thus requiring $N=2K+2H+n_{k}$ data points to be reconstructed by Padé Approximant.

Now that we know the form of the G-transform, we can calculate the form of the corresponding U-transform, eq. 2:

[TABLE]

From this last formula, it is evident that $U(z)$ is a rational function of type $[M-1/M]$ with $M=2K+2H+n_{K}-1$; the $M$ residues of its poles can only assume the values $\pm 1$ and –most important– $n_{K}-1$ is equal to the number of poles with residue $+1$ minus the number of poles of residues $-1$ i.e. the sum of the residues itself, as we stated in the Introduction.

We conclude this section by noting that the above argument works not only for signals starting at $n_{1},n_{2},\cdots,n_{K}$, but also for signals peaking at those points. This is because the growing part of each signal ends at the starting point of the decaying part, and therefore its Z-transform is a polynomial. More in detail: let us consider a signal of the form

[TABLE]

it’s G-transform will be

[TABLE]

The first addendum is a polynomial in $z$ of order $n_{k}-1$, while the second addendum is again given by eq. 5; we can therefore write

[TABLE]

which again is a $[n_{k}/1]$ rational function.

III Some practical details for the numerics

The first step of our method is to compute the Taylor series of $U(z)$ from the one (eq. 1) defining $G(z)$, so as to be able to then build its Padé Approximant, from which to get the residues of its poles, and finally $n_{K}$.

From Ref. ref46 (section $0.31$), the coefficients $c_{n}/a_{0}$ of the Taylor expansion of the ratio of two expansions

[TABLE]

are given by the recurrence formula

[TABLE]

in our case the coefficients of the series at the denominator are $a_{n}=s_{n}$ from eq. 1; those of the series at the numerator are $b_{n}=(n+1)s_{n+1}$. It is therefore convenient to define $\alpha_{n}=s_{n}/s_{0}$ and $\gamma_{n}=c_{n}/s_{0}$ so as to have a recurrence formula for the coefficients themselves:

[TABLE]

From eq. 8, we have seen that the Padé for $U(z)$ has to be be of type $[M-1/M]$, with $M=(2K+2H+n_{K}-1)$ poles. We therefore need a $2M$ term expansion for $U(z)$. We’ll have to get to $\gamma_{2M-1}$, or $\alpha_{2M}$ in the recurrence formula eq. 14, for which we need $N=2M+1$ data points from the original series, as opposed to the $M+1$ points only, needed to reconstruct $G(z)$.

The maximum possible value for $n_{K}$ corresponds to the case when the data series contains only the signal starting at $n_{K}$, meaning $K=1$ and $H=0$; therefore, $M=n_{K}^{max}+1$ and $N=2M+1=2n_{K}^{max}+3$, which gives us the bound

[TABLE]

for the possible delays detectable from a truncated sequence of $N$ data points.

IV Numerical tests

We present here some numerical examples as a test of the method described above.

Our first example considers a signal of the form

[TABLE]

to which random noise uniform in $[-\mu,\mu]$ is added. In Figure 1 we show the difference between the number of residues around $+1$ (within a radius $r=100\cdot\mu$) and the number of residues around $-1$ (again within a radius $r=100\cdot\mu$), as a function of the delay $n_{1}$, for various levels of noise between $10^{-8}$ and $10^{-3}$. The plot is in agreement with our statement in section II that the delay equals said difference. As expected, the linear dependence of the difference on the delay fails at about $N/2$.

Our next example considers the case of three signals, one of which one is delayed:

[TABLE]

In this case, the noise we add is a complex Gaussian random noise with variance $\sigma=10^{-4}$ in both real and imaginary part. In Figure 2.a we show the signal for $n_{1}=30$: as the amplitudes of the three signals are comparable, the beginning of the delayed signal is barely visible. Figure 2.b shows the poles of the Padé approximant to $U(z)$ (crosses) and the poles and zeros of of the $[n-1/n]$ Padé approximant to $G(z)$ (circles and squares, respectively): clearly something is wrong with the assumption that $G(z)$ is a rational function of the $[n-1/n]$ type: instead of having the $G(z)$ poles distributed around the unit circle, $31$ poles form a circle inside the unit circle (two more are close but not on the circle itself). Finally Figure 2.c shows the cumulative number of residues as a function of the distance from $+1$ (black line) and from $-1$ (red line), together with their difference (blue line): the latter shows a maximum of $30$ at a distance of about $0.1$.

To check our claim that the difference between the number of poles with residue $+1$ and those with residues $-1$ equals the longest delay, we now pass to the case of two delayed signals:

[TABLE]

to which complex Gaussian random noise with variance $\sigma=10^{-4}$ in both real and imaginary part is again added. In Figure 3.a we show the signal for $n_{1}=22$ and $n_{2}=30$. Figure 3.b shows the poles of the Padé approximant to $U(z)$ (crosses) and the poles and zeros of of the $[n-1/n]$ Padé approximant to $G(z)$ (circles and squares, respectively): again the assumption that $G(z)$ is a rational function of the $[n-1/n]$ type produces a circle of poles inside the unit circle. Figure 3.c shows the cumulative number of residues as a function of the distance from $+1$ (black line) and from $-1$ (red line), together with their difference (blue line). This time the difference only reaches $26$ at a distance of about $0.1$ and only gets to $30$ for a distance close to $1$.

We finally consider the case when the delayed signal is present all the time, but peaks at $n=n_{1}-1$:

[TABLE]

Again complex Gaussian random noise with variance $\sigma=10^{-4}$ in both real and imaginary part is added and we consider $n_{1}=30$. The results are very similar to those of Figure 2; as expected the maximum difference is $29$.

V Conclusions

We have shown that the information about the starting points in time of a number of delayed signals can be directly retrieved from the residues of the poles of the U-transform of the signal: they can only be $\pm 1$ and their sum plus one equals the longest delay among those of the signal’s frequency components. We have then described how an iterative algorithm using this property can be implemented in practice to locate the starting points in time of all the frequency components of a signal. Numerical examples confirm the feasibility and stability of the proposed scheme.