How to use a FIR on a signal with a drifting central frequency?

Question

I have been working on this subject for a while and I couldn't get a proper answer.

I am trying to estimate the phase of the input signal :

$x(t)=Ae^{j(at^2+bt+c)}+n(t)+i(t)$

Where $n(t)$ is a AWGN, $i(t)$ is an interfiring signal and $a,b,c$ are constant values. As you can see, the central frequency of the useful signal changes over time ($f_c(t)=2at+b$).

First, I have to filter out the interfering signal so I implemented a low-pass FIR filter H and I compute the output of the filter as follow $y[n]=\sum{h[k].x[n-k]}$. However because of the moving frequency my signal will not stay in the filter passband very long.

One solution I had was to estimate the central frequency of the signal every 1 second and shift the input signal by this estimate (which means the I add to the FIFO of my FIR the samples $x[k]e^{-2\pi f_ck/F_e}$) but when I estimate the new central frequency I will have samples with different central frequencies in my FIFO and the filter doesn't like that (I can see phase jump at the time of frequency change). So can anybody tell me how I am supposed to tackle this issue ?

Thank you !

EDIT : Thank you very much for your answers, I realized while reading your answers that I simplified my problem too much. The actual evolution of the phase of my interest signal is much more complicated than a second order polynomial. My signal is more like $x(t)=Ae^{jd(t)}+n(t)+i(t)$ with $d(t)$ a very complex function taking into account the distance and the speed of the object emitting the signal (which is a pure carrier). Also, the $i(t)$ is a function representing other object emitting other pure carrier at different speed and position (so my signal is more like $x(t)=\sum{A_ke^{jd_k(t)}}+n(t)$). I want to get the phase of each signal separately so first of all I did a small accumulation of signal to compute the spectrum of the signal and get the approximate frequency of each signal. Then I apply several FIR filters, one for each signal hoping I would get each $A_ke^{jd_k(t)}$ separately. The problem I have is that the frequency of each signal moves a lot and I have to "adapt" my FIR filter (and I don't know how to do that...) I talked about a second order phase polynomial because it is a good model of the function $d(t)$ on a short term observation.

Thank you once again for the time and effort you put in your answers.

Answer 1

Without understanding more detailed results on what i(t) is, one solution is to estimate either the phase or frequency first, and then pass that result through a filter.

An approach to do the frequency estimation digitally directly on the I (real) and Q (imaginary) samples of x(t) is the cross product frequency discriminator. After hard limiting x(t) to remove the AM components (which would result in samples on the unit circle in the complex plane), the frequency can be estimated using:

$$F \propto I_1Q_2-I_2Q_1$$

Where $I_1$,$Q_1$ and $I_2$,$Q_2$ represent successive samples. This is because the cross product is proportional to the phase between two vectors, so in this case the phase is the phase versus time, and therefore frequency. This is shown in the math and diagram below.

Frequency is the change in phase versus time ($d\phi/dt$). To determine the phase between two hard limited successive samples in time $V_1$ and $V_2$ given by:

$V_1 = e^{j\phi_1} = I_1+jQ_1$

$V_2 = e^{j\phi_2} = I_2+jQ_2$

For small angles, the phase $\phi_2-\phi_1$ can be reasonably approximated by the imaginary component of the conjugate product $(V^*_1) (V_2)$:

$V^*_1V_2 = e^{-j\phi_1}e^{j\phi_2} =e^{j(\phi_2-\phi_1)}$

For small angles $cos(\phi) \approx \phi$,

So therefore for small angles

$\phi_2-\phi_1 \approx Im(e^{-j\phi_1}e^{j\phi_2})$

$= Im((I_1-jQ_1)(I_2+jQ_2))$

$=I_1Q_2-I_2Q_1$

The above is extremely simple with low processing requirements useful for hardware implementation but depends on the small angle criteria for accurate angle estimation. This will work well as long as the sampling rate is high enough that the phase between samples doesn't deviate significantly. Alternatively if further processing is feasible, an ATAN2 operation can be instead be done on the two successive I,Q samples for full 4 quadrant accurate phase estimation (and therefore frequency given the phase difference between two successive samples over a much wider frequency range).

A simple delay and multiply frequency discriminator could also be used, followed by a filtering of the resulting frequency versus time waveform. I have further detailed an example of using a simple delay and multiply frequency discriminator in this recent answer: Demodulation of FSK signal

Both approaches require hard limiting your waveform first to remove any amplitude variation of your signal. A differentiator as done through the differencing in these approaches is a high pass filter, the hard limiter removes all AM noise, and the post detection filtering will remove the noise induced "jumps" that you see in the result.

Answer 2

This is what I would try:

Ignore the $n(t)$ and $i(t)$ since you know nothing about them. I am assuming that $x$ is complex and the rest of the variables are all real. Take the complex natural log of both sides:

$$ \ln \left( | x(t) | \right) + j \arg( x(t) ) = \ln(A) + j ( a t^2 + b t + c ) $$

Separate the real and imaginary parts. Using the imaginary part:

$$ \arg( x(t) ) = a t^2 + b t + c $$

(Or you could have taken the $\arg$ of both sides.)

This is readily solvable with a set of sample points.

$$ \begin{bmatrix} \arg( x(t_0) \\ \arg( x(t_1) \\ \arg( x(t_2) \\ : \\ \arg( x(t_{N-1}) \\ \end{bmatrix} = \begin{bmatrix} t_0^2 & t_0 & 1 \\ t_1^2 & t_1 & 1 \\ t_2^2 & t_2 & 1 \\ : \\ t_{N-1}^2 & t_{N-1} & 1 \\ \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} $$

[Note: When taking the $arg$s, you may have to keep track of windings and add multiples of $2\pi$ depending on the size of your sample frame.]

As an matrix equation it looks like this:

$$ P = BC $$

$P$ are the point values, $B$ are your basis functions, and $C$ are your coefficients you want to solve for.

All the values are real so $B^*=B^T$.

$$ B^T P = B^T B C $$

$$ ( B^T B )^{-1} B^T P = ( B^T B )^{-1} ( B^T B ) C $$

$$ C = ( B^T B )^{-1} B^T P $$

You have now solved for a best fit $(a,b,c)$

Construct a signal $F=e^{j(at^2+bt+c)}$ from these values. Call your original signal $X$.

$$ X = A F $$

$$ F \cdot X = A F \cdot F $$

$$ A = \frac{F \cdot X }{ F \cdot F } $$

You have now solved for a best fit $A$. You can subtract $AF$ from $X$ to see $ n(t) + i(t) + e(t) $ where $e(t)$ is the fit error. This may give you a way to characterize $i(t)$, find its best fit and remove it from your original signal.

Answer 3

This is my second answer. Technically it is introducing a high pass filter, but a very simple one with interesting properties for the OP's signal type.

Starting with your signal:

$$ x(t)=Ae^{j(at^2+bt+c)}+n(t)+i(t) $$

For a chosen value of $d$ (a lag time):

$$ x(t-d)=Ae^{j(a(t-d)^2+b(t-d)+c)}+n(t-d)+i(t-d) $$

Define the filtered signal as the difference of the current minus the lag. This will mitigate the noise somewhat and more importantly remove the "DC" portion of $i(t)$ completely.

$$ \begin{aligned} y(t)&= x(t) - x(t-d) \\ &= Ae^{j(at^2+bt+c)}-Ae^{j(a(t-d)^2+b(t-d)+c)} \\ &+[n(t)-n(t-d)]+[i(t)-i(t-d)] \\ &\approx Ae^{j(at^2+bt+c)}-Ae^{j(at^2-2adt+a d^2+bt-bd+c)} \\ &= Ae^{j(at^2+bt+c)}-Ae^{j(at^2+(b-2ad)t+c-bd+a d^2)} \\ &= Ae^{j(at^2+(b-ad)t+c+\frac{d(a d-b)}{2} )}\left[e^{j(ad t-\frac{d(a d-b)}{2})} - e^{j(-adt+\frac{d(a d-b)}{2})} \right] \\ &= j 2 A e^{j(at^2+(b-ad)t+c+\frac{d(a d-b)}{2} )} \sin\left(ad t-\frac{d(a d-b)}{2}\right) \end{aligned} $$

From the magnitude of $y$:

$$ | y(t) | = \left| 2A\sin\left(ad t-\frac{d(ad-b)}{2}\right) \right| $$

If you are plotting real life data for visualization this one may be good. Parameter estimation is more difficult with this approach. The remarkable thing is that the frequency of the envelope is steady. The $t^2$ is isolated to the rotation. This allows you find $(A,a,b)$, but to find $c$ you still have to use the angle values.

From the angle of $y$:

$$ \begin{aligned} arg(y(t)) &= at^2+(b-ad)t+c-\frac{d(d-b)}{2} +\frac{\pi}{2} \\ &= at^2 + b't+c' \end{aligned} $$

Proceed as with my other answer with better results. This is my recommend approach.

---

To find the amplitude best do as follows:

Solve for (a, b', c') using the technique in the other answer and then solve for $(a,b,c)$. Use those in the definition of $y(t)$ with $A=1$ to generate the $F$ basis signal. Then you can solve for $A$ with

$$ A = \frac{F \cdot Y}{F \cdot F} $$

In response to the change in the OP's problem description.

Actually, your new information makes your situation a lot more solvable. You know a lot about $i(t)$.

The $a,b,c$ parameters can be determined on segments of your signal, overlapping if necessary. This forms a piecewise definition of your new $d_k(t)$ for the $k$ signal your are evaluating. For each segment evaluation for each signal you can select an optimal $d$ value to generate the $y$ signal.

The process goes like this:

On a patch of your signal:

1) Get a rough estimate of the dominant signal, found by selecting a very small $d$.

2) Substract the estimated dominant signal from the source signal

3) Find the second signal in the difference, estimate its parameters

4) Subtract the esimated second signal from the source signal, re-estimate the first.

5) Remove the re-estimated first from the signal, re-estimate the second

6) Remove both best estimates from source the signal, find the third signal.

7) Repeat, for each signal, remove the estimated other signals, and re-estimate, look for another signal.

This will converge rapidly. It is very similar to tonal decomposition with a DFT.

Here is the kicker, if you have a rough estimate of the adjacent signals frequency on the high side, you can set $d$ to its average cycle length within the segment and virtually eliminate any residue left from its misestimation and removal. Any higher frequncy signals than that will be significantly attenuated as well.

Your segment length selections have to follow DFT rules for the number of cycles per frame differences you need for proper resolution.

When you are processing the next segment, you should be able to start off with a signal list from the previous segment that is going to be close to the correct solution. You shouldn't start from scratch on each segment.

---

By doing a balanced difference, any phase shift is eliminated. This generalizes to any FIR constructed this way. The math is also a little easier to follow.

$$ \begin{aligned} y(t)&= x(t+d) - x(t-d) \\ &= Ae^{j(a(t+d)^2+b(t+d)+c)}-Ae^{j(a(t-d)^2+b(t-d)+c)} \\ &= Ae^{j(at^2+2adt+a d^2+bt+bd+c)}-Ae^{j(at^2-2adt+ad^2+bt-bd+c)} \\ &= Ae^{j(at^2+bt+c+a d^2)}\left[e^{j(2ad t+ bd)} - e^{j(-2ad t- bd)} \right] \\ &= j 2 A e^{j(at^2+bt+c)} e^{ja d^2} \sin\left( d(2a t+ b) \right) \end{aligned} $$

You can also divide the difference by the original signal to get the "exponential factor".

$$ \frac{y(t)}{x(t)} = e^{j\left( a d^2 + \frac{\pi}{2} \right) } 2 \sin\left(d(2a t+ b)\right) $$

$$ \arg \left( \frac{y(t)}{x(t)} \right) = \arg \left( \frac{x(t+d) - x(t-d)}{x(t)} \right) = a d^2 + \frac{\pi}{2} $$

For a simple single signal, the $a$ value can be read directly from the rotation factor.