Difference between revisions of "NOAA POES TIP Demodulation"

From NebarnixWiki
Jump to navigationJump to search
(Added data plots and updated matlab script)
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. Update: added some updates to the Manchester decoder
+
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. Here's a link to the crappy matlab script I have written to do most of the data processing and extraction.
  
<pre style=" height: auto;
+
http://wiki.nebarnix.com/w/images/6/69/Decode_syncd_IQ2.m
    max-height: 200px;
 
    overflow: auto;
 
    background-color: #eeeeee;
 
    word-break: normal !important;
 
    word-wrap: normal !important;
 
    white-space: pre !important;">
 
% 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));
 
 
 
</pre>
 
  
 
Here's the final constellation of bits. It is still crappy.  
 
Here's the final constellation of bits. It is still crappy.  
Line 285: Line 44:
 
* http://nebarnix.com/sdr/POES_56k250.raw 56.250khz IEEE 32-bit float raw IQ recording of the (I think) NOAA-18 beacon at 137.35Mhz  
 
* http://nebarnix.com/sdr/POES_56k250.raw 56.250khz IEEE 32-bit float raw IQ recording of the (I think) NOAA-18 beacon at 137.35Mhz  
 
* http://nebarnix.com/sdr/pll_nocarrier_polyphased_equalized_loproto.raw IEEE 32-bit float raw IQ output from the GNU radio script at 1 symbol per sample to pull into matlab.  
 
* http://nebarnix.com/sdr/pll_nocarrier_polyphased_equalized_loproto.raw IEEE 32-bit float raw IQ output from the GNU radio script at 1 symbol per sample to pull into matlab.  
 +
 +
==Data Extraction==
 +
I have updated the matlab script to extract a few of the instruments embedded in the stream!
 +
 +
The SEM-2 (Space Environment Monitor) contains two experiments, the MEPED (Medium Energy Proton and Electron Detector) and the TED (Total Energy Detector)
 +
The following MEPED data was taken as the NOAA15 spacecraft flew north from Hungary to Scandinavia and into the auroral oval.
 +
 +
[[Image:MEPED_fixed.png|450px]]
 +
 +
 +
The DCS/2 experiment is a re transmission from mobile platforms (sea buoys, arctic fox collars, sea ice monitors, weather balloons etc etc). The format is not provided, as the instrument is owned and managed by CNES of France. I plotted the raw stream to tease out patterns in the data. What juicy mysteries lay here :)
 +
 +
[[Image:DCSdata.png|450px]]
  
 
==Todo==
 
==Todo==

Revision as of 06:56, 3 September 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. Here's a link to the crappy matlab script I have written to do most of the data processing and extraction.

http://wiki.nebarnix.com/w/images/6/69/Decode_syncd_IQ2.m

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

Data Extraction

I have updated the matlab script to extract a few of the instruments embedded in the stream!

The SEM-2 (Space Environment Monitor) contains two experiments, the MEPED (Medium Energy Proton and Electron Detector) and the TED (Total Energy Detector) The following MEPED data was taken as the NOAA15 spacecraft flew north from Hungary to Scandinavia and into the auroral oval.

MEPED fixed.png


The DCS/2 experiment is a re transmission from mobile platforms (sea buoys, arctic fox collars, sea ice monitors, weather balloons etc etc). The format is not provided, as the instrument is owned and managed by CNES of France. I plotted the raw stream to tease out patterns in the data. What juicy mysteries lay here :)

DCSdata.png

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.