"Extracting the fundamental frequency from a time-domain signal" is basically a long-winded way of saying "demodulating the baseband signal from FM-encoded signal". Well, that's still long winded, but I digress.
In an FM radio, the contraption that accomplishes this task is referred to as a "discriminator". So the task here is to implement a software-based frequency discriminator. There are a few different ways to accomplish this.
The conceptually easiest is to count zero crossings. It is straightforward and it is relatively immune to amplitude noise. However, good results require a relatively high-sample rate, in order to get good accuracy on the exact location of each zero crossing (which may occur between samples when sampled at lower sample rates). It is also not continuous: you have no information about the signal between the zero crossings---so you have to guess. I ended up linearly interpolating in my implementation with acceptable results, but some sort of IIR or FIR filter would be more appropriate for certain domains.
As an aside, if you can guarantee that only one frequency will be present (with no noise) and you are sampling at a relatively high sample rate, then you can double the resolution of the zero-crossing detector by multiplying the signal with itself. This will effectively double the frequency of the signal. However, my own experimentation has shown this technique to be useless when there is noise present.
In hardware, some discriminators use phase-locked-loop (PLL) to discipline an oscillator to match the frequency detected. The adjustments used to adjust the oscillator is the demodulated baseband signal. It is possible to simulate this with software. I have not tried doing this.
A somewhat heavy-handed approach to implementing a discriminator is to do an FFT on the signal and then use a heuristic to determine the fundamental frequency from the FFT output. While somewhat straightforward, this was something I was trying to avoid. It just seemed like there should be a better solution.
My favorite method is somewhat different than those described above. The basic idea is this:
- Break the signal into quadrature (I and Q) signals
- Use I and Q to calculate the instantaneous angle via a simple arc tangent.
- Use the difference between the last two decoded angles to calculate angular velocity
- Convert the angular velocity into a frequency with some simple linear math.
#1 is a lot more simple than you might otherwise imagine. It is
easiest to describe it with code. For each sample that is input,
into the discriminator, you do the following:
// Constants
const float sample_rate;
// Persistent state
float prev_i = 0;
float prev_q = 0;
float last_angle = 0;
int phase = 0;
// Per-sample state
float v_i = 0;
float v_q = 0;
// Our input, a sample
float sample;
// Our output, a frequency
float ret;
switch(phase) {
case 0:
v_i = sample;
v_q = 0;
phase++;
break;
case 1:
v_i = 0;
v_q = sample;
phase++;
break;
case 2:
v_i = -sample;
v_q = 0;
phase++;
break;
case 3:
v_i = 0;
v_q = -sample;
phase = 0;
break;
}
We now have our I and Q signals, but there is a great deal of high-frequency noise on them, making them useless without additional filtering. While a FIR filter (via a convolution kernel) would be the preferable filter, I found a simple 2-pole IIR filter to be adequate for SSTV demodulation and is significantly faster than a larger (but more optimal) FIR filter. I pre-calculated the coefficients for a a two-pole low-pass IIR filter with a cutoff frequency of 0.25 times the sample rate. Applying the filter is then fairly easy:
v_i = iir_float_feed(filter_i, v_i);
v_q = iir_float_feed(filter_q, v_q);
Once you have your filtered I and Q signals, things get fairly straightforward:
ret = -self->last_angle;
self->last_angle = atan2f(v_q,v_i);
ret += self->last_angle;
"But hey!", you might exclaim, "How is this 'fast' if we are relying
on that slow atan2f() call?"
This is where the magic comes in. Have a look at this beauty, which completely replaces the code above:
ret = (v_q*prev_i - v_i*prev_q)/(v_i*v_i + v_q*v_q);
prev_i = v_i;
prev_q = v_q;
What we've done is collapse steps #2 and #3 into a single operation that involves only basic math. Now, if you were wondering, this code is indeed only an approximation---but it is a pretty damn good one based on my experience.
You can then convert the angular velocity to a frequency with some simple math:
ret = (1.0f-ret/M_PI)+sample_rate;
I've also found it helpful to also apply a filter to this output signal as well.
You could also extract the instantaneous amplitude with the help of
everyone's favorite religious zealot, Pythagorus: sort(v_i^2 + v_q^2). This is, however, not terribly useful for SSTV
demodulation.
The option 4 algorithm is great. It is fast: needs no trig functions and it is a perfect candidate for SIMD operations. It is super easy to use: one sample goes in, one frequency comes out. It is quite low-latency as well, with the latency limited primarily by the I and Q filters.