Petite analyse d'une captation d'une minute avec un Orange Pi Zero et une clé RTL-SDR V3 autour de 162 MHz de 16000000 échantillons à 262144 échantillons par seconde (environ une minute donne un fichier de 32 MB).
/usr/bin/rtl_sdr -f 162000000 -s 262144 -n 16000000 core.iq
Je prends juste le quart du fichier pour accélérer le traitement dans GNU Octave (dd).
Pour lire les données IQ 8 bits fournies par rtl_sdr, j'utilise une fonction dans loadFile.m :
function y = loadFile(filename)
% y = loadFile(filename)
%
% reads complex samples from the rtlsdr file
%
fid = fopen(filename,'rb');
y = fread(fid,'uint8=>double');
y = y-127;
y = y(1:2:end) + i*y(2:2:end);
Avec plotspec.m,
% plotspec(x,Ts) plots the spectrum of the signal x
% Ts = time (in seconds) between adjacent samples in x
function plotspec(x,Ts)
N=length(x); % length of the signal x
t=Ts*(1:N); % define a time vector
ssf=(ceil(-N/2):ceil(N/2)-1)/(Ts*N); % frequency vector
fx=fft(x(1:N)); % do DFT/FFT
fxs=fftshift(fx); % shift it for plotting
subplot(2,1,1), plot(t,x) % plot the waveform
xlabel('seconds'); ylabel('amplitude') % label the axes
subplot(2,1,2), plot(ssf,abs(fxs)) % plot magnitude spectrum
xlabel('frequency'); ylabel('magnitude') % label the axes
je peux observer que j'ai bien enregistré des messages et qu'ils se situent bien à -25 KHz et + 25 KHz autour de 162 MHz.
Je vais maintenant tacher d'isoler un message à +25 KHz (up.iq8) et un autre à -25 KHz (down.iq8). Pour ce faire, je vais visualiser la magnitude du signal.
Zoomons...
Zoomons encore...
On peut maintenant extraire juste les messages du fichier core.iq8. Sachant qu'il y a deux bytes par échantillon et que le graphe nous donne les offsets absolus.
$ # 134000..141000 & 274000..281000
$ dd if=core.iq8 of=up.iq8 bs=2000 skip=134 count=7
$ dd if=core.iq8 of=down.iq8 bs=2000 skip=274 count=7
En cherchant un peu, on trouve un message à +25 Khz dont on s'intéresse à l'évolution de la phase :
En regardant la phase et en zommant :
Et la même chose pour un message à -25 KHz :
qui donne une évolution de la phase descendante :
On voit que, dans les deux cas, on a des sauts quand la phase passe au dessus de 360° ou en dessous de 0°. On va rendre tout cela continu pour, ensuite, soustraire/ajouter le changement de phase dû au +/- 25 KHz.
% y = contzup(phi)
%
% rend la phase continue
%
function y = contzup(phi)
N=length(phi); % length of the signal x
phi=phi+(pi/2);
y=phi;
t=0;
for i = 1:N-1
if ( phi(i) > phi(i+1) )
t = t + 2*pi;
endif
y(i+1) = (phi(i+1) + t)-(((2*pi)*25000*i)/262144);
% y(i+1) = (phi(i+1) + t);
endfor
Détail, le début avec le préambule et les premiers octets.
avec le grid à 9600 samples/sec :
Comme je ne parviens pas vraiment à écrire un programme me permettant d'obtenir les trames NMEA à partir des IQ sortant de la clé RTL-SDR en bouquinant sur la théorie. J'envisage de partir d'un programme qui fonctionne et de l'analyser... Par exemple,
rtl_ais donne d'assez bons résultats. D'un terril au dessus de Liège avec une antenne dipôle, je capte des bateaux sur la Meuse jusque Tihange (Huy) dans un sens et au delà de Maastricht dans l'autre. Et même quelques message provenant du côté de Genk, sur le canal Albert. Le code n'est pas très long et est assez lisible mais, c'est vraiment de la sorcellerie. Il va falloir du temps pour s'en imprégner et comprendre son fonctionnement.
À suivre...
