How to improve an FFT AM demodulator for NOAA APT satellite imagery?

Question

My question involves proper FFT Technique for AM demodulation, in the special case where the modulation rate exceeds the carrier frequency.

This question has a specific hobbyist application, the decoding of weather satellite pictures using minimal computing hardware and inexpensive receivers.

Warning: this is necessarily long, and you might want to skim it to see if you have time to read it. Or, perhaps just focus on what seems problematic.

APT: How the NOAA satellites send low-resolution pictures

A simplified version of the APT encoding process is as follows:

1. An 8 bit grayscale image that includes sync data, Visual and IR imagery and test/telemetry bars is fed to an 8 bit D/A converter at a rate of 2 lines per second, where each line is 2080 bytes, for an encoding rate of 4160 bytes/sec.
2. The D/A converter AM modulates a 2400 hz tone (and therefore the modulation rate exceeds the carrier frequency) creating an analog audio stream.
3. The analog audio stream FM modulates an RF carrier around 137 Mhz.

The above 1-3 summarizes more detailed descriptions at Wikipedia: Automated Picture Transmission, and NOAA KLM Users Guide 4.2 - APT System

The satellites also send high resolution data on an S-band microwave downlink, but that is beyond the scope of this question.

Hobbyist decoding process

1. The rtl_fm program of the RTL-SDR package is used to command an inexpensive dongle to receive the RF signal and record audio samples of the decoded FM to a .wav file at a rate of 20800 samples/sec (5 x data rate). The resolution of the RTL-SDR dongle is 8 bits/sample. The raw digital audio is converted to a .wav file with sox
2. A python program takes the real-valued digital audio sample from the .wav, demodulates the AM data from the audio, and arrange the decoded bytes to produce the received grayscale image.

Most of the programs for step #2 are inexpensive but proprietary, thus the desire to develop an open source decoder. Also open source is atpdec by Thierry Leconte F4DWV but it would not gracefully deal with missing samples and I was curious to see how far I could get writing my own.

For completeness, here is a stable github repository for the version of my code under discussion and the APT audio signal of the pictures posted.

I hope to describe the procedures adequately and include snippets below.

AM RMS Demodulation

The most sensitive method I have found so far is as follows:

1. Sample the real-valued audio signal at 5x the digital data encoding rate (20800 samples/sec)
2. Work with one video line length of samples (10400 samples) at a time.
3. In each line, group the samples into groups of 5, to correspond to the same transmitted digital value. Square and sum the samples within the groups of 5. The idea is that if S(t) = A cos(wt) then (S^2)(t) = A^2 (cos^2)(wt) = 0.5 A^2 + 0.5 A^2 cos(2wt) from the cos^2 trig identity and omega with the NOAA satellites is such that 5 omega is a little bigger than 180 degrees, so maybe the last term summed up over the 5 samples is small enough to get a noisy A out if we are lucky.
4. (Noise reduction) For each line of video, deduct the minimum sum of squares from each of the sums. If you consider a noisy discrete signal R(t) = S(t) + e(t) where e(t) ~ i.i.d. N (0, sigma^2) is a gaussian noise then it is clear that an RMS estimator of A is biased upwards by gaussian noise. Since we know some of the As in a line of video are low (near black), especially in the sync pulses, then we might want to push the estimates back down towards 0 using the minimum of all the estimates as a poor, but practical, estimate of noise.
5. Then multiply by 2/5 and take the square root of each of the sum-of-squares to obtain the amplitudes for each byte.
6. Digitize linearly to 8 bits.
7. (Advanced) Jiggle to try to get the 5x binning of samples to line up with what was sent by examining sync signals etc. Estimate signal phase and use that to apply a correction which greatly reduces the near-vertical-line visual artifacts.

Python snippet:

# _G5() reshapes an array to 2D - groups of 5 - i.e. (len/5 x 5)
ssq = np.sum(_G5(np.square(self.signal[start:end])), axis=1)
if denoise is True:
    denoiseAmt = np.min(ssq)
    ssq -= denoiseAmt
if phase is None:  # None means do not apply phase correction
    demod = np.sqrt(0.4*ssq)

Here is a picture decoded with this AM RMS method without step #7.

Although the picture is not perfectly aligned because of a doppler-like effect, and quite noisy, it is obviously the southeastern USA and a portion of the Atlantic Ocean.

FFT-based AM demodulation

Since the modulation rate exceeds the tone frequency, I decided to try to demodulate the upper sideband only.

The FFT-based demodulator is coded to operate as follows:

1. Have the code work with one video line at a time (10400 samples)
2. Create a sample of size 2^15 = 32768 consisting of zeros, the 10400 samples, and more zeroes so that the audio samples are in the middle of the zeroes.
3. apply python's real FFT. numpy.fft.rfft
4. find the bin corresponding to the 2400 hz audio carrier.
5. shift the FFT left so that the zero'th bin corresponds to 2400 hz.
6. zero the zeroth bin, and pad with zeroes on the right according to the shift.
7. Take the inverse real FFT numpy.fft.irfft
8. Remove the demodulated signal from the middle as the inverse of step #2
9. Downsample by 5 by taking every 5th sample and throwing the others away
10. Digitize

Python code:

def _demodAM_by_FFT(self, start=None, end=None, warn=True, minpower2=15):
    if start is None:
        start = 0
    if end is None:
        end = len(self.signal)
    # find a power of 2 that encloses the sample, then go up another power
    # minimum window of 2^15 = 32768
    power2 = max(minpower2, 1+math.ceil(math.log(end-start)/math.log(2.0)))
    fftsize = 2**power2
    binNOAACarrier = math.floor(fNOAACarrierHz*fftsize/self.rate)
    padding = fftsize-(end-start)
    padl = padding/2
    padr = padding-padl
    sample = np.concatenate( (np.zeros(padl, np.float),\
                              self.signal[start:end],\
                              np.zeros(padr, np.float)) )
    sample_fft = np.fft.rfft(sample)
    high_bin = np.argmax(np.abs(sample_fft))
    if warn:
        if high_bin!=binNOAACarrier:
            print "warning: in _demodAM_by_FFT unexpected fft peak"
            print " peak occurs in bin ",high_bin
            print " expected carrier in bin ",binNOAACarrier
        else:
            print "fft peak as expected"
    cutoff = binNOAACarrier
    demod_fft = np.concatenate( (np.copy(sample_fft[cutoff:]),\
                                 np.zeros(cutoff, np.complex) ) )
    demod_fft[0] = 0.0
    demod = np.fft.irfft(demod_fft)[padl:(padl+end-start)]
    demod5 =_G5(demod)
    demod5x2 = demod[:,2]
    return demod5x2

If one looks very carefully, the test/telemetry bars can be identified but the image is either not there or badly distorted.

The FFT based image does not change as the centered, zero-padded sample line is increased from 2^15 to 2^17. Demodulating the entire signal instead of line by line also does not noticeably improve the image.

Although data analysis plays a role in some of my work, I do not usually use FFT for anything. So I must admit some unfamiliarity with proper set up of FFT.

I am looking for advice on whether an FFT demodulation method is salvageable for noisy signals such as this, and what techniques the code above is missing (such as windowing techniques or pre-filtering the samples in a certain way) to have a better chance at success.

FOLLOW-UP

For completeness and comparison, here are some images obtained from the Hilbert transform technique suggested by the answer below.

Image from abs value of hilbert transform, single video line signal inputs

And if instead, all of the signal data is processed in one operation, (about 19M samples for a 15 minute LEO satellite pass), and then a median filter is applied before downsampling, we get something a bit clearer:

Answer 1

As @Jason R says, your _demodAM_by_FFT method pretty much just applies a (bad) filter to the raw data. It isn't an AM demodulator.

To demodulate the signal, you need to apply a hilbert transform. You can either implement this yourself, or you could use the scipy.signal.hilbert function, since you are already using scipy. Implementing the hilbert transform is usually done by way of the FFT.

You need to do (more or less) the following:

1. Hilbert transform of a set of samples

2. For each sample, take the square root of the sum of the squares of the real and imaginary parts of the hilbert result for that sample.

3. (2.) is your demodulated signal, still with 5 samples per pixel.

Those three steps recover the modulated amplitude from the audio signal. I may have explained that poorly. Lets try some pseudocode:

hilb=scipy.signal.hilbert(signal) demod=numpy.abs(hilb)

I don't think you will want to do the whole signal at once (way too much data for one FFT,) but breaking it into chunks of a 1024 (or 4096) samples should be OK. Reassemble before filtering and downsampling.

You could at this point apply a low pass filter (cutoff over 1040Hz) and resample to get your image data.

You could alternatively apply some filtering at this point to clean up the signal before doing the lowpass and downsampling. A median filter might be a good place to start (scipy.signal.medfilt)

I would first implement a proper demodulator with the hilbert function and get that working before I started filtering.

There's probably more to the decoding once you get the demodulated data, but demodulation is the first step.

In case you are interested in implementing the hilbert transform with the FFT, here is a site that I found to be fairly easy to understand.