Difference between revisions of "NOAA POES TIP Demodulation"
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 = ' | + | hfile = 'pll_nocarrier_polyphased_equalized_loproto2.raw'; |
− | + | %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); | ||
− | + | %% | |
− | %% | ||
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 |
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 |
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 |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
end | end | ||
Line 108: | Line 128: | ||
end | end | ||
− | %% plot autocorrelation to see if anything is here | + | %% 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)) ' | + | 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
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
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.
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
Raw IQ Recordings
- 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.
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!
- 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.