Difference between revisions of "NOAA POES TIP Demodulation"

From NebarnixWiki
Jump to navigationJump to search
(added matlab code)
Line 29: Line 29:
 
I decided to process the data in Matlab because I wanted access to all of the analytical tools until I could figure out what I was doing.  
 
I decided to process the data in Matlab because I wanted access to all of the analytical tools until I could figure out what I was doing.  
  
(post matlab script here)
+
(code posted at bottom of page)
  
 
Here's the final constellation of bits. It is still crappy.  
 
Here's the final constellation of bits. It is still crappy.  
Line 45: Line 45:
  
 
* Make the GNU radio script more tolerant of noise. Airplanes flying overhead transmit ACARS packets that blast the sat signal out of the water from time to time. The clock sync seems to lose lock easily and doesn't like to re-aquire it after a few hard knocks.
 
* Make the GNU radio script more tolerant of noise. Airplanes flying overhead transmit ACARS packets that blast the sat signal out of the water from time to time. The clock sync seems to lose lock easily and doesn't like to re-aquire it after a few hard knocks.
 +
 +
==Matlab Source==
 +
 +
<pre>
 +
% reads  complex samples from the rtlsdr file
 +
%
 +
clear all;
 +
clf;
 +
%hfile = 'POES_56k250.raw';
 +
hfile = 'pll_nocarrier_polyphased_equalized_loproto.raw';
 +
Nfft = 108473;
 +
fid = fopen(hfile,'rb');
 +
y = fread(fid,'float32');
 +
%y = fread(fid);
 +
 +
%%
 +
fs = 56.25e3/3.0638856476079346557759626604434;
 +
t=0 : 1/fs : (Nfft-1)/fs;
 +
%realY = y(1:2:end);
 +
%imagY = y(2:2:end);
 +
%y = y(1:2:end) + 1i*y(2:2:end);
 +
y = y(2:2:end) + 1i*y(1:2:end);
 +
 +
%% Plot Constellation
 +
clf;
 +
colormap copper;
 +
x = y(Nfft:floor(1.2*Nfft));
 +
c = linspace(1,10,length(x));
 +
scatter(real(x),imag(x),[],c);
 +
 +
%% Convert constellation to bits
 +
clear bitstream;
 +
for idx = 1:numel(y)
 +
    if(y(idx) > 0)
 +
        bitstream(idx) = '1';
 +
    else
 +
        bitstream(idx) = '0';
 +
    end
 +
end
 +
 +
%% Convert to bits from raw manchester bits
 +
idx2=1;
 +
clockmod = 0;
 +
idxerr=1;
 +
clear errx erry bitstream_manchester;
 +
for idx = 2:numel(bitstream)
 +
 +
    %check for bit boundary, if true, then grab the bit
 +
    if(mod(idx,2) == clockmod)
 +
        bitstream_manchester(idx2) = bitstream(idx);
 +
        idx2 = idx2+1;                     
 +
        %if(bitstream(idx-1) == bitstream(idx))
 +
        %    errx(idxerr)=idx2;
 +
        %    erry(idxerr)=real(y(idx));
 +
        %    idxerr=idxerr+1;
 +
        %end       
 +
    else
 +
        if(bitstream(idx-1) == bitstream(idx))
 +
            errx(idxerr)=idx2;
 +
            erry(idxerr)=real(y(idx));
 +
            idxerr=idxerr+1;
 +
        end       
 +
    end       
 +
   
 +
    %look for two same bits in a row
 +
    %this can only happen on bit boundary
 +
    %if we find this, resync to this boundary
 +
    %if(bitstream(idx-1) == bitstream(idx))
 +
    %    if(bitstream(idx) == 0)
 +
    %        clockmod = mod(idx,2);
 +
    %    end
 +
    %end
 +
end
 +
 +
%% plot autocorrelation to see if anything is here (broken, needs 1,0 not '1','0')
 +
subplot(2,1,1);
 +
plot(xcorr(bitstream))
 +
subplot(2,1,2);
 +
plot(xcorr(bitstream_manchester))
 +
 +
%% test autocorrelation of random number generator as reference
 +
testarray = rand(1,100000);
 +
for idx = 1:numel(testarray)
 +
    if(idx>0.5)
 +
        testarray2(idx) = 1;
 +
    else
 +
        testarray2(idx) = 0;
 +
    end
 +
end
 +
plot(xcorr(testarray2))
 +
 +
%% Look for syncword and its inverse (in case of phase reversal)
 +
S1 = '11101101111000100000100';
 +
S2 = '00010010000111011111011';
 +
 +
k1 = strfind(bitstream_manchester, S1);
 +
k2 = strfind(bitstream_manchester, S2);
 +
fprintf([ '\n' num2str(numel(k2)+numel(k1)) ' detected\n' num2str(sum(mod(diff(k1),832)==0)) ' verified\n\n']);
 +
 +
subplot(2,1,1);
 +
plot(k1,1:length(k1),'o',k2,1:length(k2),'x',errx,erry*10,'.');
 +
subplot(2,1,2);
 +
plot(k1(2:end),diff(k1),'o',k2(2:end),diff(k2),'x',errx,erry*10,'.');</pre>

Revision as of 08:03, 28 August 2015

Problem Statement

While listening to NOAA APT signals I noticed a strange telemetry signal at 137.35 and 137.77 that tracked NOAA POES sats. At first I thought it was a second sat, but then I realized that the doppler shifts matched too closely. I found the handbook for the satellite constellation, which makes mention of a DBD or TIM signal

POES waterfall.jpg

Research

I can't actually find any instance of anyone using SDR to decode this signal. I decided that this needed to change :)

I found these documents which describe the modulation and format [1] [2]

Modulation Format

I discovered that this is a very old modulation scheme, dating back to at least the late 60's. The scheme is a BPSK that uses a modulation index of 67 degrees to preserves the carrier for tracking purposes. It is also Manchester encoded, which means that even though the stream is 8320 bps, it is encoded using double the number of symbols, or 16,640 sps.

GNU Radio Script

I wrote a GNU radio workflow to prepare the data for demodulation. Coming soon.

File Source->Throttle->AGC2->Carrier Tracking PLL->High Pass Filter (to remove carrier)->Polyphase Clock Sync w/ lowpass proto->CMA Equalizer->File Sink (post GNU radio file and block image here)

Matlab Script

I decided to process the data in Matlab because I wanted access to all of the analytical tools until I could figure out what I was doing.

(code posted at bottom of page)

Here's the final constellation of bits. It is still crappy.

POES Constellation.jpg

The circles are the locations where the sync word, "11101101 1110001 0000" appear. The bottom plot is a diff of the first, showing intervals of 832bits (minor frame length) which proves this is good data. The green dots are where there are Manchester encoding errors (two 1's or two 0's in a row not on a bit boundary). Note that there is no error correct -- unless I find a better way to clean up the phase constellation (or get a better signal) -- or weed out 'barely' bad bits, I probably don't have more than a couple perfect minor frames... it takes 32 seconds or 3200 minor frames to make up a whole major frame. I have some work left to do smile emoticon

Syncword manchester errors.jpg

Todo

  • Get better data. The signal is weak, only about 1 watt, which means that getting 32 seconds (major frame) of GOOD clean data is somewhat difficult. Alternatively I could make the demodulator more robust, though I'm not quite sure how to do that just yet.
  • Write a routine to pull the bytes apart using the sync word and then shuffle them into ident bins. This shouldn't be hard, just tedious.
  • Make the GNU radio script more tolerant of noise. Airplanes flying overhead transmit ACARS packets that blast the sat signal out of the water from time to time. The clock sync seems to lose lock easily and doesn't like to re-aquire it after a few hard knocks.

Matlab Source

% reads  complex samples from the rtlsdr file
%
clear all;
clf;
%hfile = 'POES_56k250.raw';
hfile = 'pll_nocarrier_polyphased_equalized_loproto.raw';
Nfft = 108473;
fid = fopen(hfile,'rb');
y = fread(fid,'float32');
%y = fread(fid);

%%
fs = 56.25e3/3.0638856476079346557759626604434;
t=0 : 1/fs : (Nfft-1)/fs;
%realY = y(1:2:end);
%imagY = y(2:2:end);
%y = y(1:2:end) + 1i*y(2:2:end);
y = y(2:2:end) + 1i*y(1:2:end);

%% Plot Constellation
clf;
colormap copper;
x = y(Nfft:floor(1.2*Nfft));
c = linspace(1,10,length(x));
scatter(real(x),imag(x),[],c);

%% Convert constellation to bits
clear bitstream;
for idx = 1:numel(y)
    if(y(idx) > 0)
        bitstream(idx) = '1';
    else
        bitstream(idx) = '0';
    end
end

%% Convert to bits from raw manchester bits
idx2=1;
clockmod = 0;
idxerr=1;
clear errx erry bitstream_manchester;
for idx = 2:numel(bitstream)

    %check for bit boundary, if true, then grab the bit
    if(mod(idx,2) == clockmod)
        bitstream_manchester(idx2) = bitstream(idx);
        idx2 = idx2+1;                       
        %if(bitstream(idx-1) == bitstream(idx))
        %    errx(idxerr)=idx2;
        %    erry(idxerr)=real(y(idx));
        %    idxerr=idxerr+1;
        %end        
    else
        if(bitstream(idx-1) == bitstream(idx))
            errx(idxerr)=idx2;
            erry(idxerr)=real(y(idx));
            idxerr=idxerr+1;
        end        
    end        
    
    %look for two same bits in a row
    %this can only happen on bit boundary
    %if we find this, resync to this boundary
    %if(bitstream(idx-1) == bitstream(idx))
    %    if(bitstream(idx) == 0)
    %        clockmod = mod(idx,2);
    %    end
    %end
end

%% plot autocorrelation to see if anything is here (broken, needs 1,0 not '1','0')
subplot(2,1,1);
plot(xcorr(bitstream))
subplot(2,1,2);
plot(xcorr(bitstream_manchester))

%% test autocorrelation of random number generator as reference
testarray = rand(1,100000);
for idx = 1:numel(testarray)
    if(idx>0.5)
        testarray2(idx) = 1;
    else
        testarray2(idx) = 0;
    end
end
plot(xcorr(testarray2))

%% Look for syncword and its inverse (in case of phase reversal)
S1 = '11101101111000100000100';
S2 = '00010010000111011111011';

k1 = strfind(bitstream_manchester, S1);
k2 = strfind(bitstream_manchester, S2);
fprintf([ '\n' num2str(numel(k2)+numel(k1)) ' detected\n' num2str(sum(mod(diff(k1),832)==0)) ' verified\n\n']);

subplot(2,1,1);
plot(k1,1:length(k1),'o',k2,1:length(k2),'x',errx,erry*10,'.');
subplot(2,1,2);
plot(k1(2:end),diff(k1),'o',k2(2:end),diff(k2),'x',errx,erry*10,'.');