【无标题】ADC-B报文识别&CRC检码&采样率转换

zz不倦

已于 2024-04-09 15:15:51 修改

阅读量219

点赞数 2

分类专栏：通信 ADS-B检测&定位&教研&译码文章标签： matlab 信息与通信

于 2024-04-09 15:12:20 首次发布

本文链接：https://blog.csdn.net/lzh_made_wall/article/details/137550025

版权

通信 ADS-B检测&定位&教研&译码专栏收录该内容

1 篇文章 0 订阅

订阅专栏

加载一段1090ES信号，对标准112位报文(详情请自查)通过阈值进行识别，并且提取标准报文段。

clc
clear all
file = "D:\ADS_B_data.dat";
fid = fopen(file, 'r');
N = 327680;
s_array = fread(fid,[327680, 375],'int16');

% 数据处理采样率转换

% 12.8MHz*5/4=16MHz
si = 18;
thr = 80;
s = s_array(:, si);
s = s(2:2:end) + sqrt(-1) * s(1:2:end);
Fs = 1.28e7;
T = 1/Fs;
s_n = 0: T: T*(length(s)-1);
s_amp = abs(s);
%%　插值滤波　滤波插值
% interp
inter_Fs = 5*Fs;
inter_T = 1/inter_Fs;
int_s = zeros(5*length(s), 1);
int_s(1:5:end) = s;
int_s_amp = abs(int_s);
int_n = 0:inter_T:(5*length(s)-1)*inter_T;
%filt1
v = lowpass(int_s, Fs*5/64 *10, inter_Fs);  % 10MHz
v_amp = abs(v);
%filt2
v1 = lowpass(v, Fs*5/64 *5, inter_Fs);  % 5MHz
v1_amp = abs(v1);
% deci
dec_Fs = 5*Fs/4;
dec_T = 1/dec_Fs;
dec_s = v1(1:4:end);
dec_n = 0: dec_T: (5*length(s)/4-1)*dec_T;
dec_amp = abs(dec_s);
figure(1); stem(dec_amp, "b", "filled"); title("Decimated Sequence");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 根据门限值检测LEP（上升沿）并记录点索引
LEP_index = [];
for i = 1:length(dec_amp)
    if dec_amp(i) > thr
        LEP_index = [LEP_index, i];
    end
end
% disp(length(LEP_index));
% 记录连续continue_n个有效样本值的点为有效高脉冲VPP（高脉冲带宽）
VPP_index = [];
continue_n = 5;
for i = 1: length(LEP_index)-continue_n  % 防止超出索引
    if sum(LEP_index(i:i+4))-continue_n*LEP_index(i)==10
        VPP_index = [VPP_index, LEP_index(i)];
    end
end

% 选择只记录VPP首值索引
repet_index = []; % 记录VPP中其他索引值
for i = 2:length(VPP_index)
    if VPP_index(i-1)+1 == VPP_index(i) 
       repet_index = [repet_index, i];
    end
end
VPP_index(repet_index) = []; % 将VPP中其他索引值置0
%disp(VPP_index);
% VPP间隔点数
Delta_ns = [];
for i = 2: length(VPP_index)
    Delta_n = VPP_index(i) - VPP_index(i-1);
    Delta_ns = [Delta_ns, Delta_n];
end
%disp(Delta_ns);
% VPP 间隔时间
VPP_n = dec_n(VPP_index);
Delta_ts = [];
for i = 2: length(VPP_n)
    Delta_t = VPP_n(i) - VPP_n(i-1);
    Delta_ts = [Delta_ts, Delta_t];
end
% disp(sum(Delta_ts));
% 根据报头间距记录报头位置 % 根据间距查找尾部
%disp(Delta_ts);
Err_toler = 3*T;  % 误差三个点
head_index = [];
end_index = [];
for i = 1: length(Delta_ts)-2
    if abs(Delta_ts(i)-1e-6) < Err_toler & abs(Delta_ts(i+1)-2.5*1e-6) < Err_toler & abs(Delta_ts(i+2)-1e-6) < Err_toler
        head_index = [head_index, VPP_index(i)];
        for j = i+4:length(Delta_ts)
            if Delta_ts(j) > 2.2*1e-6
               end_index = [end_index, VPP_index(j)];
               break        
            end
        end
    end
end
% 补齐end_index最后一值
if length(head_index) ~= length(end_index) 
    end_index = [end_index, VPP_index(end)];
end
% 去除非112位长度的报文编码
Code_Len = end_index-head_index;
head_index(Code_Len < 112*16) = []; % 将不符合长度的去掉
end_index(Code_Len < 112*16) = [];
disp([si, head_index, end_index, end_index - head_index]);
if isempty(head_index)
    error("未检出合格ADS_B编码！");
end
% 求取自适应阈值　＝　一个VPP幅值均值与底噪差的一半并加上底噪均值。
signal_mean = sum(dec_amp(head_index(1):head_index(1)+7))/8;
noise_mean = sum(dec_amp(head_index(1)-8:head_index(1)-1))/8;
thr_ada = (signal_mean-noise_mean)/2+noise_mean; 
disp(thr_ada);
% 转换成二值信号
binary_sequence = [];
for i = head_index(1): head_index(1)+16*120-1
    if dec_amp(i) > thr_ada
        binary_sequence(i-head_index(1)+1) = 1;
    else 
        binary_sequence(i-head_index(1)+1) = 0;
    end
end    
%disp(binary_sequence)
figure(2);subplot(211); plot(binary_sequence); %axis([-20, 1050, -0.5, 3]); title("Binary sequence");
% 计算报文编码（前高后低为1，前低后高为0）
ADS_B_code = [];
for i = 16*8+1: 16: length(binary_sequence)
    if mean(binary_sequence(i:i+7)) > 0.8 & mean(binary_sequence(i+8:i+15)) < 0.2
        ADS_B_code = [ADS_B_code, 1];
    else 
        ADS_B_code = [ADS_B_code, 0];
    end
end
%disp(ADS_B_code);
figure(2);subplot(212); stem(ADS_B_code, "b", "filled"); 
%Code_ME_Array = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
Code_DF = bin2dec(num2str(ADS_B_code(1:5)));  % 数值到char
if Code_DF == 17
    disp("DF=17 用于S模式应答机发出ADS-B报文");
    Code_ME = ADS_B_code(33: 88);
    ME_Type = bin2dec(num2str(Code_ME(1:5)));



elseif Code_DF == 18
    disp("DF=18 用于非S模式应答机发出ADS-B报文或TIS-B报文！")
elseif Code_DF == 19
    disp("DF=19用于军事用途！")
else
    disp([Code_DF, "报文格式不符合！"])
end

%disp(ADS_B_code(end-23: end))

获取到标准112位报文后需进行 CRC校验（原理待更新…）

K = ADS_B_code;
G = [1 1 1 1  1 1 1 1  1 1 1 1  1 0 1 0  0 0 0 0  0 1 0 0  1];% 25位的生成多项式
R = mod2div(K, G); % 模二除法函数
if ~R
    disp("报文校验通过！");
    %disp(R);
else
    disp("报文校验不通过！");
    disp(R);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 模二除法函数
function r = mod2div(K, G)  % 余数 被除数 除数
    while length(K) ~= length(G)-1
        if K(1) == 1
            r = xor(K(1: length(G)), G);
            K = [r(2:length(G)), K(length(G)+1: end)];
        else
            r = xor(K(1: length(G)), zeros(1, length(G)));
            K = [r(2:length(G)), K(length(G)+1: end)];
    
        end
    end
    r = K; % 被除数位数不足除数除数时即为余数
end