Daily Archives: 9/28/2011

An impractical classic SSTV decoder…

A few days ago, I posted a .WAV file for a classic 8s SSTV image and asked if anyone could decode it. Nobody replied (I wasn’t surprised) so I set about writing my own demodulator.

Since I’m inherently lazy, here was my idea: generate the complex signal using my previously debugged Hilbert transform code. Then, for every adjacent pair of complex numbers, determine the angle between the adjacent samples by normalizing the numbers (so they have unit length) and then computing the dot product between the two. If you’ve done as much raytracing as I have, you know that the dot product of two vectors is product of their lengths and the cosine of the angle between them. So, you take the dot product, pass that through the acos function, and you’ll get an angle between 0 and $latex \frac{\pi}{2}$. You can then multiply that by the sample rate (and divide through by $latex 2\pi), and you can recover the frequency.

So, that’s what I did. If the frequency was in the band between 1500 and 2300 hz, I mapped that into a grayscale value. For values which were below 1500, I set the pixels to black. For values > 2300, I set the pixels to white. I generated one output pixel for every input sample. That means the image is too wide, but that’s no big deal: I resized it to 120 columns wide using some of the netpbm utilities (command line graphics programs).

And, here’s the image:

That’s a bit teeny, so let’s blow it up to 512×512:

Not too bad for a first try. Some thoughts: the sync pulses seem really wide to me. The spec says the sync pulses should be 5ms wide, which means they are about 7.5% as the line length. That’s a lot of the image to leave behind. I’ll have to work on that. I’m not currently doing anything to detect the vertical or horizontal sync pulses, but I’ve done that before with my NOAA weather satellite code, so I don’t anticipate any real problems. All in all, not a bad first try.

Addendum: The large image shows some JPEG compression artifacts since I expanded it from the small JPEG image without trying to save it at high quality. I can do better.

Addendum2: I decreased the overall filter size to length 11. It didn’t work very well, because the original file was at 22050, which pushes makes the portion of the spectrum we are interested in down into the region of the filter reponse which isn’t very uniform. But if we resample the sound file down to 8Khz, we filter works just fine with just 11 taps. Check it out:

The image is slightly skewed because the line length is no longer an integer number of samples, and the accumulated error causes the sync to drift slightly. A real implementation would track this accurately.

Difficulties with the Hilbert Transform…

Well, it wasn’t so much a difficulty with the Hilbert transform as a difficulty with my understanding. But with the help of my good friend Tom, my understanding was soon put right, and I thought it might make an interesting (in other words, horribly boring to anyone but myself) post, and at the very least, it would be good for me to document this little bit of mathematical fun.

Warning: complex numbers ahead!

This came up as a subtask of my slow scan television decoder. I’ll probably have a more detailed posting regarding the underlying theory, but the basic idea is as follows. The “pure” (without noise) SSTV signal is a frequency modulated signal with constant amplitude. But what does that really mean? After all, the instantaneous value of the signal varies: it’s a sine wave after all. What’s useful is to think of the analytic signal. Instead of being a single value at each moment in time, we think of a signal as a complex number, with a real and imaginary part. Let’s say the real part is $latex a$, and the imaginary part is $latex b$. The analytic signal can be represented as the complex number $latex a+b\imath$. The amplitude of the signal is then the $latex \sqrt{a^2+b^2}$.

If doesn’t make any sense to you, don’t fret. It’s mostly because we do a horrible job of teaching complex numbers, and I haven’t really done any better. You might be thinking, “what does this have to do with actual signals? The signals that I have are described by a single value, which varies over time. Why are you hurting my brain with complex numbers?”

In an attempt to stimulate you to think about it more, I’ll respond with the Socratic method, and as “what does the amplitude of a signal mean if they are described by only a single value?” In particular, how can you tell if a signal is “constant amplitude”? What does that even mean?

Okay, I’ll drop that topic for another time. Let’s just imagine that you accept my definition of constant amplitude, and we want to construct a complex signal from the real valued signal we receive from the radio. Let’s just pretend that the input signal gets copied into the real portion of the complex signal. In other words, the put values become the values for $latex a$ in our $latex a+b\imath$ above. For the purposes of our discussion, let’s say that the signal is some sine wave, $latex sin(t)$, where $latex t$ is a time variable (we could scale it to make it higher and lower frequency if we like). Let’s say the we’d like the amplitude of that signal is 1. Can we find values for $latex b$ that makes the amplitude $latex \sqrt{a^2+b^2} = 1$? If we square both sides, we see that $latex a^2 + b^2 = 1$, and since $latex a = sin(t)$, we see that $latex b^2 = 1-sin^2(t)$. If you remember any of your high school trig, you might remember that make $latex b = cos(t)$ ($latex sin^2(t) + cos^2(t) = 1$).

What this means is that for any given sine wave at a given frequency, we can get the complex version of the signal by copying a cosine signal of the same frequency into the imaginary part. If you think of cosine as just a phase with a 90 degree phase shift, you’ll see that if we had a magic filter that could make a 90 degree phase shift, you could easily construct the complex signal from the real signal.

But what if our signal isn’t made up of sine waves? Well, you can decompose the signal into the sum of a bunch of different sines (and cosines) of varying amplitude (proof left for the motivated reader). Applying the magic 90 degree phase shifter to all the sines, and then adding them back up, we’ll get the imaginary part of the signal we need.

None of this has anything to do with the problem really, it’s just background. It’s remarkably hard to explain this stuff to anyone who doesn’t understand it already, and I fear I’ve not done anything good. Oh well, onto my problem.

The way that you can implement the 90 degree phase shifter is by using some magic called the Hilbert transform, which can be implemented as an FIR filter. I wanted to create a black box that I feed values of my real signal in, and get (after some delay) the complex signal out. The code isn’t all that hard really, I’ll go ahead and put it here:

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>

#define NZEROS  (50)

static float xv[NZEROS+1];
static float yv[NZEROS+1];

static float xcoeffs[NZEROS+1] ;

double
window(int n, int N)
{
    double a0 = 0.355768 ;
    double a1 = 0.487396 ;
    double a2 = 0.144232 ;
    double a3 = 0.012604 ;

    return a0 - a1 * cos(2.0 * M_PI * n / (N - 1))
              + a2 * cos(4.0 * M_PI * n / (N - 1)) 
              - a3 * cos(6.0 * M_PI * n / (N - 1)) ;
}

void 
oddhilb (float *hilb, int n)
{
    float f;
    int i, j, k, m;

    m = j = k = (n - 1) >> 1;

    for (i = 1; i <= m; i += 2) {
        f = 2.0 / (i * M_PI) ;

        hilb[j++] = 0.0;
        hilb[j++] = f;
        hilb[k--] = 0.0;
        hilb[k--] = -f;
    }
    /* and now... window it...
     */
    FILE *fp = fopen("unwindowed.dat", "w") ;
    for (i=0; i<=n; i++)
        fprintf(fp, "%lf\n", hilb[i]) ;
    fclose(fp) ;

    for (i=0; i<=n; i++) {
        hilb[i] *= window(i, n+1) ;
    }

    fp = fopen("windowed.dat", "w") ;
    for (i=0; i<=n; i++)
        fprintf(fp, "%lf\n", hilb[i]) ;
    fclose(fp) ;
}

complex 
filterit(float val)
{
    float sum; 
    int i;

    for (i = 0; i < NZEROS; i++) {
        xv[i] = xv[i+1];
        yv[i] = yv[i+1];
    }
    xv[NZEROS] = val ;

    for (i = 0, sum = 0.; i <= NZEROS; i++) 
        sum += (xcoeffs[i] * xv[i]);
    yv[NZEROS] = sum ;
    return xv[NZEROS/2] + I * yv[NZEROS] ;
}

main(int argc, char *argv[])
{
    float f ;
    complex c ;
    int i ;
    double t = 0. ;

    oddhilb(xcoeffs, NZEROS+1) ;

    for (i=0; i<4000; i++) {
        f = sin(t) ;
        t += 0.5 ;
        c = filterit(f) ;
        printf("%lf %lf\n", creal(c), cimag(c)) ;
        // printf("%lf\n", cabs(c)) ;
    }

}

The program implements the Hilbert transform as an FIR filter. At each step, we compute a new value $latex f$ for our input signal. We then pass it to our filter stage, which returns our complex valued signal. If we then plot the real and imaginary parts of our signal, we should get a circle. And, in fact we do, now:

But when I started, I didn’t. The radius of the circle varied over quite a bit (maybe 10% or so). The reason was I forgot to apply the windowing function to the Hilbert transform filter coefficients. Plotting the two sets of coefficients, we see that the filtered versions fall to zero faster away from the middle.

We can use the FFT on these coefficients to see what the frequency response is.

You can see that the unwindowed filter has very uneven response over the range of frequencies, while the windowed signal is nice and flat over all but the region very near zero and the Nyquist frequency. This means that when we feed it nice sine waves, we’ll get a nice analytic signal with very constant amplitude.

Addendum(s): There is a small “off by one” error that adds an extra zero coefficient to the filters above. I’ll tidy it up a bit more when I get a chance. As Tom pointed out, for use in my SSTV decoder, it’s almost certain that just 11 coefficients would work just fine (and be faster). We could also optimize out the multiply by zeros in the convolution.