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:
- 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.
- The D/A converter AM modulates a 2400 hz tone (and therefore the modulation rate exceeds the carrier frequency) creating an analog audio stream.
- 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
- 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
- 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:
- Sample the real-valued audio signal at 5x the digital data encoding rate (20800 samples/sec)
- Work with one video line length of samples (10400 samples) at a time.
- 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. - (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.
- Then multiply by 2/5 and take the square root of each of the sum-of-squares to obtain the amplitudes for each byte.
- Digitize linearly to 8 bits.
- (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:
- Have the code work with one video line at a time (10400 samples)
- 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.
- apply python's real FFT. numpy.fft.rfft
- find the bin corresponding to the 2400 hz audio carrier.
- shift the FFT left so that the zero'th bin corresponds to 2400 hz.
- zero the zeroth bin, and pad with zeroes on the right according to the shift.
- Take the inverse real FFT numpy.fft.irfft
- Remove the demodulated signal from the middle as the inverse of step #2
- Downsample by 5 by taking every 5th sample and throwing the others away
- 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: