大话水声通信技术---（BFSK仿真）

       在之前的理论篇中，笔者梳理了水声通信相关的理论知识体系。本次笔者给出了一套基于BFSK的水声通信系统，该系统已经在实际的硬件中得到了验证。

%******************************通信声呐仿真BPSK方式**********************
%几点假设：
%（1）基于射线声学理论
%（2）几何衰减按球面波传播衰减规律衰减，不考虑吸收衰减
%（3）考虑水面和水底的反射（浅海声信道？）
%（4）考虑在高斯白噪声背景下
%（5）整个空间声速分布均匀

close all;
clear all;
clc;

c = 1500;            %声速             Unit：m/s
f_lfm = 28e3;        %信号频率         Unit：Hz
f1 = 28e3;           %代表1  第29个点
f0 = 26e3;           %代表0  第27个点
fs = 128e3;        %采样率 最高频率的5倍   Unit：Hz 
Ts = 1/fs;

Code_w_points = 128;           %码元宽度所占的点数  
Code_w = Code_w_points*Ts;     %码元宽度  总是不小于码元信号宽度，因为有多途
Tao = 0.5e-3;                  %码元信号宽度  保证信号是整周期发射的
non_data_points = 5*Code_w_points;           %导频信号与数据信号隔的点数

Sample_time = 0.2;   %采样时长         Unit：s      假设信号发射的时刻为零时刻
SNR = 40;             %信噪比           Unit：dB
v = 0;               %水听器运动速度 正数表示运动方向靠近发射环能器，负数表示远离

match_signal_points = 1024;           %匹配信号点数 1024*1  3 
Tao_Pilot = match_signal_points*Ts;     %导频信号信号脉宽         Unit：s
delta_f = 6e3;             %带宽

BT = delta_f*Tao_Pilot;    %时间增益
k = delta_f/Tao_Pilot;     %频率变化率  Hz/s

%一帧数据的点数 = 导频信号点数+无信号点数+7码元信号点数+128点的尾部噪声数据
Code_num = 14*5;                                         %码元个数
Code_start = match_signal_points + non_data_points;      %一帧数据中编码数据开始的点索引值
one_Frame_points = Code_start + Code_num*Code_w_points + 128;  
Bytes_num = fix(Code_num/14);                            %一帧包含的字节数
%***********************线性分组码相关的参数**********************
Code = [16*10+10 
        16*1+1
        16*2+2
        16*3+3
        16*4+4];    %待编码的数据

odd_even = 1;     %奇校验
%基本参数
n_c = 7;           %输出码字长度
k_c = 4;           %输入信息组长度
d_c = 3;           %分组码最小距离
rate_c = k_c/n_c;  %码率

%定义校验矩阵（也可以用Matlab函数直接生成）
H_check=[ 1 0 0 1 0 1 1;
          0 1 0 1 1 1 0;
          0 0 1 0 1 1 1];

%定义生成矩阵
G_gener=[ 1 1 0 1 0 0 0;
          0 1 1 0 1 0 0;
          1 1 1 0 0 1 0;
          1 0 1 0 0 0 1];

%**********************模拟水池的环境************************
H = 5;       %水深                    Unit：m
H1 = 1;    %发射换能器水深           Unit：m       1
H2 = 1;    %接收换能器水深           Unit：m       1
D = 10;    %接收与发射换能器水平距离  Unit：m   10

Re_coef_surf = -0.5;     %水面反射系数  
Re_coef_bottom = 0.01;   %水底反射系数   
Reflex_num = 4;         %考虑最大的反射次数
%**********************************************************
%%
%****************************************************************************************************************************************
%*********************************************************编码与调制**********************************************************************
%****************************************************************************************************************************************
sample_num = fix(Sample_time*fs);       %采样总点数
nTs = (0:sample_num-1)/fs;              %离散的采样时刻

nTs1 = (0:Ts:(Tao_Pilot-Ts)) - Tao_Pilot/2;
match_signal = sin(2*pi*f_lfm*nTs1 + pi*k*nTs1.^2);
% match_signal_points = size(match_signal,2);             %匹配信号点数

%生成水听器接收数据
receive_signal_Pilot = Com_Sonar_data_generation(c,fs,f_lfm*(1+2*v/c),0,BT,Tao_Pilot,Sample_time,140,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);   %生成导频信号 
receive_signal0_0 = Com_Sonar_data_generation(c,fs,f0*(1+2*v/c),0/180*pi,1,Tao,Sample_time,140,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);  %生成码元信号 '0'
receive_signal1_0 = Com_Sonar_data_generation(c,fs,f1*(1+2*v/c),0/180*pi,1,Tao,Sample_time,140,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);  %生成码元信号 '1'

Code_E = Linear_Block_Encode(Code,G_gener);        %编码数据
sig0 = [receive_signal0_0;receive_signal1_0];
receive_signal = BFSK_Modulation(Code_E,receive_signal_Pilot,sig0,match_signal_points + non_data_points,Code_w_points);  %二元频移键控调制

receive_signal = [receive_signal((end-499):end),receive_signal(1:(end-500))];     %在前面补零
receive_signal = receive_signal + (10^(-SNR/20))*randn(1,sample_num);             %加入一定信噪比的噪声

xcorre_xy = conv(receive_signal,match_signal(end:-1:1));                  %匹配滤波
xcorre_xy = xcorre_xy(1:(end-match_signal_points+1));                     %裁剪掉后match_signal_points个数据


%%
%****************************************************************************************************************************************
%********************************************************提取包络************************************************************************
%****************************************************************************************************************************************
%***************************低通滤波 Fp=2KHz ,Fstop=10KHz 64阶  140KHz
hn_m = [0.0004    0.0005    0.0008    0.0011    0.0013    0.0015    0.0015    0.0012    0.0007   -0.0001   -0.0013   -0.0028...
        -0.0045   -0.0063   -0.0080   -0.0093   -0.0101   -0.0100   -0.0088   -0.0062   -0.0022    0.0034    0.0105    0.0189...
        0.0284    0.0385    0.0488    0.0588    0.0678    0.0756    0.0814    0.0851    0.0863    0.0851    0.0814    0.0756...
        0.0678    0.0588    0.0488    0.0385    0.0284    0.0189    0.0105    0.0034   -0.0022   -0.0062   -0.0088   -0.0100...
        -0.0101   -0.0093   -0.0080   -0.0063   -0.0045   -0.0028   -0.0013   -0.0001    0.0007    0.0012    0.0015    0.0015...
        0.0013    0.0011    0.0008    0.0005    0.0004];

abs_xcorre_xy = abs(xcorre_xy);                      %取绝对值
LPF_abs_xcorre_xy =  filter(hn_m,1,abs_xcorre_xy);     %低通滤波取包络

%%
%****************************************************************************************************************************************
%***********************************************搜索信号到达的时刻************************************************************************
%****************************************************************************************************************************************
threshold = 0.12;
arrive_set = Get_Frame_Start_Time(LPF_abs_xcorre_xy,match_si