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

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.

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.

• Hello, thanks for the complete answer. However I feel like you focus on how to estimate the frequency of the signal (which is indeed an important point but I solved that issue with a simple FFT computation of my signal). The complexity of the question relies in the filtering operation for a frequency drifting signal (I edited my question by the way). – Ryan Nov 28 '19 at 12:54
• Hi @Ryan- Yes I think Cedron's answer is best in that it is specifically answering your question on how to estimate PHASE not frequency as I focused. A simple FFT will give you the best estimate in the presence of white noise of only the average frequency over the FFT's time sample block. Cendron is giving you the best solution for the presence of a frequency ramp (linear fit of a freq slope since your question suggested quadratic phase). Regardless in either case my recommendation with regards to the noise holds in that you can filter after you translate to phase. His LMS solution does that. – Dan Boschen Nov 28 '19 at 13:09
• @Ryan So I reread your edits--- it sounds like this is an FMCW radar and you would like to resolve distance versus time? In that case frequency would be your unit of interest (frequency = range) and then I think my answer would be very applicable to you since it is giving you frequency versus unit time. The noise on this answer is reduced with additional filtering but the challenge in any case is you have unknown dynamics so your measurement (filter) needs to pass anything within the possible dynamics you could have. Anything else you should filter. – Dan Boschen Nov 28 '19 at 13:14
• The more you know about the dynamics the better you can filter. – Dan Boschen Nov 28 '19 at 13:14

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.

• I like this answer best! It IS the least squares estimator to what the phase is versus time – Dan Boschen Nov 27 '19 at 14:44
• @DanBoschen Thanks. Hang on though, I'm working on a second "slick trick" option which should be better. – Cedron Dawg Nov 27 '19 at 14:59
• Hi ! This answer is indeed very good, however it works only if you can neglect i(t) because if you cannot isolate your signal, your discriminator ($arg(x(t0)$) will not work properly. Which is why I wanted to add a FIR filter around my signal of interest and this is where my troubles start, I don't know how to filter a signal if its frequency changes over time. I edited my question to clarify it. – Ryan Nov 28 '19 at 12:59
• @Ryan My second answer should largely address your concerns, check it out. As Dan mentions, you may be better off modelling $d_k(t)$ as piecewise linear with standard DFT techniques. – Cedron Dawg Nov 28 '19 at 16:20

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.

$$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.

• Hi Cedron, great material (still digesting!). Just a detail (correct if I'm wrong), in the $y(t)$ development shouldn't it be $\frac{d(b-d)}{2}$ (instead of $\frac{d(d-b)}{2})$? – David Nov 30 '19 at 16:58
• @David Good catch, I negated the term rather than switch the order. There was also a missing $a$ on the $d^2$ term. I have corrected it above. See also the followup about balanced differences. Yes, lots to digest. Thanks. – Cedron Dawg Nov 30 '19 at 18:08