Difference between revisions of "NOAA POES TIP Demodulation"

From NebarnixWiki
Jump to navigationJump to search
Line 29: Line 29:
 
==Matlab Script==
 
==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.  
+
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. Update: added some updates to the Manchester decoder
  
 
<pre style=" height: auto;
 
<pre style=" height: auto;
Line 43: Line 43:
 
clf;
 
clf;
 
%hfile = 'POES_56k250.raw';
 
%hfile = 'POES_56k250.raw';
hfile = 'pll_nocarrier_polyphased_equalized_loproto.raw';
+
hfile = 'pll_nocarrier_polyphased_equalized_loproto2.raw';
Nfft = 108473;
+
%hfile = 'pll_nocarrier_polyphased_equalized_loproto.raw';
 
fid = fopen(hfile,'rb');
 
fid = fopen(hfile,'rb');
 
y = fread(fid,'float32');
 
y = fread(fid,'float32');
 
%y = fread(fid);
 
%y = fread(fid);
  
%%
+
%fs = 56.25e3/3.0638856476079346557759626604434;
fs = 56.25e3/3.0638856476079346557759626604434;
+
%t=0 : 1/fs : (Nfft-1)/fs;
t=0 : 1/fs : (Nfft-1)/fs;
 
 
%realY = y(1:2:end);
 
%realY = y(1:2:end);
 
%imagY = y(2:2:end);
 
%imagY = y(2:2:end);
 
%y = y(1:2:end) + 1i*y(2:2:end);
 
%y = y(1:2:end) + 1i*y(2:2:end);
 
y = y(2:2:end) + 1i*y(1:2:end);
 
y = y(2:2:end) + 1i*y(1:2:end);
 
+
%%
%% Plot Constellation
 
 
clf;
 
clf;
 
colormap copper;
 
colormap copper;
 +
Nfft = 100000+Nfft;
 
x = y(Nfft:floor(1.2*Nfft));
 
x = y(Nfft:floor(1.2*Nfft));
 
c = linspace(1,10,length(x));
 
c = linspace(1,10,length(x));
 
scatter(real(x),imag(x),[],c);
 
scatter(real(x),imag(x),[],c);
 
+
%%
%% Convert constellation to bits
+
%convert constellation to bits
 
clear bitstream;
 
clear bitstream;
 
for idx = 1:numel(y)
 
for idx = 1:numel(y)
     if(y(idx) > 0)
+
     if(real(y(idx)) > 0)
 
         bitstream(idx) = '1';
 
         bitstream(idx) = '1';
 
     else
 
     else
Line 74: Line 73:
 
end
 
end
  
%% Convert to bits from raw manchester bits
+
%% convert to bits from raw manchester bits
 
idx2=1;
 
idx2=1;
 
clockmod = 0;
 
clockmod = 0;
 
idxerr=1;
 
idxerr=1;
 +
%idxerr2=1;
 
clear errx erry bitstream_manchester;
 
clear errx erry bitstream_manchester;
 
for idx = 2:numel(bitstream)
 
for idx = 2:numel(bitstream)
 
+
    if(mod(idx,2) ~= clockmod)   
 +
        if(bitstream(idx-1) == bitstream(idx))
 +
            errx(idxerr)=idx2;
 +
            erry(idxerr)=real(y(idx));
 +
            idxerr=idxerr+1; 
 +
            if(abs(real(y(idx-1))) > 0.8 && abs(real(y(idx))) > 0.8)               
 +
                clockmod = mod(idx,2); %only resync if we have confidence in the bit decisions
 +
            end
 +
           
 +
        end       
 +
    end
 +
   
 
     %check for bit boundary, if true, then grab the bit
 
     %check for bit boundary, if true, then grab the bit
 
     if(mod(idx,2) == clockmod)
 
     if(mod(idx,2) == clockmod)
         bitstream_manchester(idx2) = bitstream(idx);
+
         if(abs(real(y(idx))) > abs(real(y(idx+1)))) %use the strongest symbol to determine bit
 +
            if(bitstream(idx) == '1')               
 +
                bitstream_manchester(idx2) = '0';
 +
            else
 +
                bitstream_manchester(idx2) = '1';
 +
            end
 +
        else
 +
            if(bitstream(idx+1) == '1')               
 +
                bitstream_manchester(idx2) = '1';
 +
            else
 +
                bitstream_manchester(idx2) = '0';
 +
            end
 +
           
 +
        end
 +
        %bitstream_manchester(idx2) = bitstream(idx+1);
 
         idx2 = idx2+1;                       
 
         idx2 = idx2+1;                       
 +
       
 
         %if(bitstream(idx-1) == bitstream(idx))
 
         %if(bitstream(idx-1) == bitstream(idx))
 
         %    errx(idxerr)=idx2;
 
         %    errx(idxerr)=idx2;
 
         %    erry(idxerr)=real(y(idx));
 
         %    erry(idxerr)=real(y(idx));
 
         %    idxerr=idxerr+1;
 
         %    idxerr=idxerr+1;
         %end      
+
         %end          
    else
 
        if(bitstream(idx-1) == bitstream(idx))
 
            errx(idxerr)=idx2;
 
            erry(idxerr)=real(y(idx));
 
            idxerr=idxerr+1;
 
        end       
 
 
     end         
 
     end         
 
      
 
      
Line 108: Line 128:
 
end
 
end
  
%% plot autocorrelation to see if anything is here (broken, needs 1,0 not '1','0')
+
%% plot autocorrelation to see if anything is here
 
subplot(2,1,1);
 
subplot(2,1,1);
 
plot(xcorr(bitstream))
 
plot(xcorr(bitstream))
Line 131: Line 151:
 
k1 = strfind(bitstream_manchester, S1);
 
k1 = strfind(bitstream_manchester, S1);
 
k2 = strfind(bitstream_manchester, S2);
 
k2 = strfind(bitstream_manchester, S2);
fprintf([ '\n' num2str(numel(k2)+numel(k1)) ' detected\n' num2str(sum(mod(diff(k1),832)==0)) ' verified\n\n']);
+
fprintf([ '\n' num2str(numel(k2)+numel(k1)) ' detected\n' num2str(sum(mod(diff(k1),832)==0) + sum(mod(diff(k2),832)==0)) ' match length\n\n' num2str(idxerr) ' errors\n\n']);
  
 
subplot(2,1,1);
 
subplot(2,1,1);
 
plot(k1,1:length(k1),'o',k2,1:length(k2),'x',errx,erry*10,'.');
 
plot(k1,1:length(k1),'o',k2,1:length(k2),'x',errx,erry*10,'.');
 
subplot(2,1,2);
 
subplot(2,1,2);
plot(k1(2:end),diff(k1),'o',k2(2:end),diff(k2),'x',errx,erry*10,'.');</pre>
+
plot(k1(2:end),diff(k1),'o',k2(2:end),diff(k2),'x',errx,erry*10,'.');
 +
 
 +
 
 +
%% convert ascii to binary at matched syncword locations
 +
for frameIdx=1:numel(k1)
 +
    for frameByteIdx=0:103
 +
        byte=0;
 +
        for bit_idx=0:7     
 +
          if(bitstream_manchester(k1(frameIdx)+frameByteIdx*8+bit_idx)=='0')
 +
              byte = bitshift(byte,1);
 +
          else
 +
              byte = bitshift(byte,1);
 +
              byte = bitor(byte,1);
 +
          end
 +
        end
 +
        minorFrames(frameIdx,frameByteIdx+1)=byte;
 +
    end   
 +
end
 +
 
 +
plot(bitor(bitshift(bitand(minorFrames(:,5),1),8),minorFrames(:,6))) %plot minor frame counter
 +
 
 +
%%
 +
%Pull out HIRS3 data
 +
%16,17,
 +
%22,23,
 +
%26,27,
 +
%30,31,
 +
%34,35,
 +
%38,39,
 +
%42,43,
 +
%54,55,
 +
%58,59,
 +
%62,63,
 +
%66,67,
 +
%70,71,
 +
%74,75,
 +
%78,79,
 +
%82,83,
 +
%84,85,
 +
%88,89,
 +
%92,93
 +
 
 +
%offset by 1 because we start at 1 not zero
 +
%15,16,
 +
%21,22,
 +
%25,26,
 +
%29,30,
 +
%33,34,
 +
%37,38,
 +
%41,42,
 +
%53,54,
 +
%57,58,
 +
%61,62,
 +
%65,66,
 +
%69,70,
 +
%73,74,
 +
%77,78,
 +
%81,82,
 +
%83,84,
 +
%87,88,
 +
%91,92
 +
 
 +
HIRSdata(1,:) = (minorFrames(:,15));
 +
HIRSdata(2,:) = (minorFrames(:,16));
 +
 
 +
HIRSdata(3,:) = (minorFrames(:,21));
 +
HIRSdata(4,:) = (minorFrames(:,22));
 +
 
 +
HIRSdata(5,:) = (minorFrames(:,25));
 +
HIRSdata(6,:) = (minorFrames(:,26));
 +
 
 +
HIRSdata(7,:) = (minorFrames(:,29));
 +
HIRSdata(8,:) = (minorFrames(:,30));
 +
 
 +
HIRSdata(9,:) = (minorFrames(:,33));
 +
HIRSdata(10,:) = (minorFrames(:,34));
 +
 
 +
HIRSdata(11,:) = (minorFrames(:,37));
 +
HIRSdata(12,:) = (minorFrames(:,38));
 +
 
 +
HIRSdata(13,:) = (minorFrames(:,41));
 +
HIRSdata(14,:) = (minorFrames(:,42));
 +
 
 +
HIRSdata(15,:) = (minorFrames(:,53));
 +
HIRSdata(16,:) = (minorFrames(:,54));
 +
 
 +
HIRSdata(17,:) = (minorFrames(:,57));
 +
HIRSdata(18,:) = (minorFrames(:,58));
 +
 
 +
HIRSdata(19,:) = (minorFrames(:,61));
 +
HIRSdata(20,:) = (minorFrames(:,62));
 +
 
 +
HIRSdata(21,:) = (minorFrames(:,65));
 +
HIRSdata(22,:) = (minorFrames(:,66));
 +
 
 +
HIRSdata(23,:) = (minorFrames(:,69));
 +
HIRSdata(24,:) = (minorFrames(:,70));
 +
 
 +
HIRSdata(25,:) = (minorFrames(:,73));
 +
HIRSdata(26,:) = (minorFrames(:,74));
 +
 
 +
HIRSdata(27,:) = (minorFrames(:,77));
 +
HIRSdata(28,:) = (minorFrames(:,78));
 +
 
 +
HIRSdata(29,:) = (minorFrames(:,81));
 +
HIRSdata(30,:) = (minorFrames(:,82));
 +
 
 +
HIRSdata(31,:) = (minorFrames(:,83));
 +
HIRSdata(32,:) = (minorFrames(:,84));
 +
 
 +
HIRSdata(33,:) = (minorFrames(:,87));
 +
HIRSdata(34,:) = (minorFrames(:,88));
 +
 
 +
HIRSdata(35,:) = (minorFrames(:,91));
 +
HIRSdata(36,:) = (minorFrames(:,92));
 +
 
 +
</pre>
  
 
Here's the final constellation of bits. It is still crappy.  
 
Here's the final constellation of bits. It is still crappy.  
Line 151: Line 287:
  
 
==Todo==
 
==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.
+
* 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. Update: Reddit user _COD32_ was able to capture several minutes of very high quality data to play with! This has already helped a lot.  
  
* 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.  
+
* 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. Update: I have pulled the minor frames apart and plotted minor frame counter data. Look at that counter go! [[Image:Minor_frame_counter.png|250px]]<br clear=both />
  
 
* 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.

Revision as of 09:54, 31 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. http://nebarnix.com/sdr/nocarrier_lowpass_polyphase_equalizer.grc

POES GRC flow.JPG

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. Update: added some updates to the Manchester decoder

% reads  complex samples from the rtlsdr file
%
clear all;
clf;
%hfile = 'POES_56k250.raw';
hfile = 'pll_nocarrier_polyphased_equalized_loproto2.raw';
%hfile = 'pll_nocarrier_polyphased_equalized_loproto.raw';
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);
%%
clf;
colormap copper;
Nfft = 100000+Nfft;
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(real(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;
%idxerr2=1;
clear errx erry bitstream_manchester;
for idx = 2:numel(bitstream)
    if(mod(idx,2) ~= clockmod)    
        if(bitstream(idx-1) == bitstream(idx))
            errx(idxerr)=idx2;
            erry(idxerr)=real(y(idx));
            idxerr=idxerr+1;   
            if(abs(real(y(idx-1))) > 0.8 && abs(real(y(idx))) > 0.8)                
                clockmod = mod(idx,2); %only resync if we have confidence in the bit decisions
            end
            
        end        
    end
    
    %check for bit boundary, if true, then grab the bit
    if(mod(idx,2) == clockmod)
        if(abs(real(y(idx))) > abs(real(y(idx+1)))) %use the strongest symbol to determine bit
            if(bitstream(idx) == '1')                
                bitstream_manchester(idx2) = '0';
            else
                bitstream_manchester(idx2) = '1';
            end
        else
            if(bitstream(idx+1) == '1')                
                bitstream_manchester(idx2) = '1';
            else
                bitstream_manchester(idx2) = '0';
            end
            
        end
        %bitstream_manchester(idx2) = bitstream(idx+1);
        idx2 = idx2+1;                       
        
        %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
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) + sum(mod(diff(k2),832)==0)) ' match length\n\n' num2str(idxerr) ' errors\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,'.');


%% convert ascii to binary at matched syncword locations
for frameIdx=1:numel(k1)
    for frameByteIdx=0:103
        byte=0;
        for bit_idx=0:7       
           if(bitstream_manchester(k1(frameIdx)+frameByteIdx*8+bit_idx)=='0')
               byte = bitshift(byte,1);
           else
               byte = bitshift(byte,1);
               byte = bitor(byte,1);
           end
        end
        minorFrames(frameIdx,frameByteIdx+1)=byte;
    end    
end

plot(bitor(bitshift(bitand(minorFrames(:,5),1),8),minorFrames(:,6))) %plot minor frame counter

%% 
%Pull out HIRS3 data
%16,17,
%22,23,
%26,27,
%30,31,
%34,35,
%38,39,
%42,43,
%54,55,
%58,59,
%62,63,
%66,67,
%70,71,
%74,75,
%78,79,
%82,83,
%84,85,
%88,89,
%92,93

%offset by 1 because we start at 1 not zero
%15,16,
%21,22,
%25,26,
%29,30,
%33,34,
%37,38,
%41,42,
%53,54,
%57,58,
%61,62,
%65,66,
%69,70,
%73,74,
%77,78,
%81,82,
%83,84,
%87,88,
%91,92

HIRSdata(1,:) = (minorFrames(:,15));
HIRSdata(2,:) = (minorFrames(:,16));

HIRSdata(3,:) = (minorFrames(:,21));
HIRSdata(4,:) = (minorFrames(:,22));

HIRSdata(5,:) = (minorFrames(:,25));
HIRSdata(6,:) = (minorFrames(:,26));

HIRSdata(7,:) = (minorFrames(:,29));
HIRSdata(8,:) = (minorFrames(:,30));

HIRSdata(9,:) = (minorFrames(:,33));
HIRSdata(10,:) = (minorFrames(:,34));

HIRSdata(11,:) = (minorFrames(:,37));
HIRSdata(12,:) = (minorFrames(:,38));

HIRSdata(13,:) = (minorFrames(:,41));
HIRSdata(14,:) = (minorFrames(:,42));

HIRSdata(15,:) = (minorFrames(:,53));
HIRSdata(16,:) = (minorFrames(:,54));

HIRSdata(17,:) = (minorFrames(:,57));
HIRSdata(18,:) = (minorFrames(:,58));

HIRSdata(19,:) = (minorFrames(:,61));
HIRSdata(20,:) = (minorFrames(:,62));

HIRSdata(21,:) = (minorFrames(:,65));
HIRSdata(22,:) = (minorFrames(:,66));

HIRSdata(23,:) = (minorFrames(:,69));
HIRSdata(24,:) = (minorFrames(:,70));

HIRSdata(25,:) = (minorFrames(:,73));
HIRSdata(26,:) = (minorFrames(:,74));

HIRSdata(27,:) = (minorFrames(:,77));
HIRSdata(28,:) = (minorFrames(:,78));

HIRSdata(29,:) = (minorFrames(:,81));
HIRSdata(30,:) = (minorFrames(:,82));

HIRSdata(31,:) = (minorFrames(:,83));
HIRSdata(32,:) = (minorFrames(:,84));

HIRSdata(33,:) = (minorFrames(:,87));
HIRSdata(34,:) = (minorFrames(:,88));

HIRSdata(35,:) = (minorFrames(:,91));
HIRSdata(36,:) = (minorFrames(:,92));

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

Raw IQ Recordings

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. Update: Reddit user _COD32_ was able to capture several minutes of very high quality data to play with! This has already helped a lot.
  • 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. Update: I have pulled the minor frames apart and plotted minor frame counter data. Look at that counter go! Minor frame counter.png
  • 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.