How to estimate the frequency characteristics of a nonlinear chirp?

Given a partial time-domain recording of a noise-corrupted chirp that asymptotically approaches a frequency, and the function that determines the non-linearity, how can I accurately estimate the final frequency it approaches?

Here's an example of what I'm talking about. A radio receiver hears noise until around sample x=5100, then the chirp starts at f=0, and quickly sweeps up to an unknown frequency, but then fades out around x=5500 to some annoying filter ringing caused by a hardware front-end band-pass filter. X: sample count, Y: ADC count

Ideal solution will have the following characteristics:

1. Accurately predict the final frequency
2. Work over a range of final frequencies
3. Be invariant to different starting frequencies
4. Be invariant to different stopping frequencies
5. Relatively quick to compute

Approaches Tried So Far:

1. Estimating instantaneous frequency using zero-crossing, then curve-fitting. It seems that noise and filter-ringing causes errors in the frequency calculation that make it very sensitive to slight changes in waveform shape. It works better but not well if the signal is resampled to a higher sample-rate first.

2. Estimating instantaneous frequency using autocorrelation, then curve-fitting using the strongest lags from autocorrelation Function. This works well given knowledge of the incoming frequency so that the autocorrelation window size can be picked ahead of time, but fails when the ACF window isn't wide enough to see the whole signal. If the ACF is too wide, it'll sometimes see a harmonic as the strongest lag.

3. Curve-fitting the raw signal. This tends not to converge if there's even a little error in the initial conditions, or noise in the signal.

4. FFT/STFT, then curve-fitting the bins. What I ran into here was that it seemed the bins were destroying the shape of the signal. I couldn't find a good size of bin that gave enough frequency-resolution and time-resolution.

An example recording is available here.

• Have you tried a PLL-based approach? – Marcus Müller Dec 3 '18 at 11:06
• I haven't tried a PLL-based approach yet. – rsaxvc Dec 3 '18 at 20:09

How about just using $$\dfrac{d\phi}{dt}$$ to compute the instantaneous frequency?

GNU Octave code (I just made up a sample rate of 2 Msps):

y = hilbert(x);
z = y - mean(real(y));
dphi_dt = diff(unwrap(arg(z)));
Fs = 2e6;
f_inst = (Fs/2)/pi * dphi_dt;
k = [0:length(f_inst)-1];
t = k/Fs;
plot(k, f_inst);
axis([5000 5800 -2e5 2e5]);
title('Instantaneous Frequency');
xlabel('Sample index');
ylabel('Frequency (Hz)');
grid on; • I think this is exactly what I'm looking for and will likely accept the answer in a few days once I've had time to test it. – rsaxvc Dec 3 '18 at 21:28
• Ok. I played a little fast and loose with removing the DC component of the analytic signal and also in just winging in tbe hilbert function to make the analytic signal. Normally I build a single sided complex BPF that goes from DC to $+Fs/2$. – Andy Walls Dec 3 '18 at 22:25
• There is also a slightly more clever way to compute $d\phi/dt$ without having to use diff() and unwrap(). One takes the $\mathrm{Arg}(z[n]z^*[n-1])$, but I didn't feel like writing a loop in Octave/Matlab. – Andy Walls Dec 3 '18 at 22:31
• I should also mention that the first difference is an approximate derivative filter, that is best close to DC but isn't very good the closer you get to $\pm F_s/2$. A good derivative filter isn't hard, but how to choose/design is another separate question that requires your performamce requirememts s to be better specified. – Andy Walls Dec 3 '18 at 22:53
• Ok. But because your input signal was real, you can not make a distinction between positive and negative frequencies. So with only real input, you can only estimate frequencies in $[0, +F_s/2]$. If you have both I & Q on input, then you can talk about negative frequencies too. – Andy Walls Dec 4 '18 at 3:25