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

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

%******************************通信声呐仿真BPSK方式**********************%几点假设：%（1）基于射线声学理论%（2）几何衰减按球面波传播衰减规律衰减，不考虑吸收衰减%（3）考虑水面和水底的反射（浅海声信道？）%（4）考虑在高斯白噪声背景下%（5）整个空间声速分布均匀close all;clear all;clc;c = 1500; %声速 Unit：m/sf_lfm = 28e3; %信号频率 Unit：Hzf1 = 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：dBv = 0;%水听器运动速度 正数表示运动方向靠近发射环能器，负数表示远离match_signal_points = 1024; %匹配信号点数 1024*1 3 Tao_Pilot = match_signal_points*Ts;%导频信号信号脉宽 Unit：sdelta_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+116*2+216*3+316*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：mH1 = 1; %发射换能器水深 Unit：m 1H2 = 1; %接收换能器水深 Unit：m 1D = 10; %接收与发射换能器水平距离 Unit：m 10Re_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阶 140KHzhn_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_signal_points,threshold);%%%绘图figure;subplot(2,1,1);plot(nTs,receive_signal);string = ['导频信号与调制的信号，信噪比SNR=',num2str(SNR),'dB'];title(string);subplot(2,1,2);plot(nTs,xcorre_xy,'-',nTs(arrive_set-32),xcorre_xy(arrive_set-32),'r*');title('相关运算');nF = (0:sample_num-1)*fs/sample_num;figure;subplot(2,1,1);plot(nTs,abs_xcorre_xy,'-',nTs(arrive_set-32),abs_xcorre_xy(arrive_set-32),'r*');title('取绝对值后');subplot(2,1,2);plot(nF,abs(fft(abs_xcorre_xy)));title('取绝对值后的频谱');figure;subplot(2,1,1);plot(nTs,LPF_abs_xcorre_xy,'-',nTs(arrive_set),LPF_abs_xcorre_xy(arrive_set),'r*');title('低通滤波后');subplot(2,1,2);plot(nF,abs(fft(LPF_abs_xcorre_xy)));title('低通滤波后频谱');%%%****************************************************************************************************************************************%*************************************************信道均衡（在低信噪比的时候效果更好）******************************************************%****************************************************************************************************************************************%根据搜索确定的到达时刻为起始时刻进行信道均衡M = 128;%FIR滤波器的系数个数 start_point = arrive_set-32-match_signal_points+1; %取第一帧的起始时刻 temp_signal = receive_signal(start_point:(start_point+match_signal_points-1)); %重合部分的导频信号desired_signal = match_signal; %期望信号 w = Channel_equalization(temp_signal,desired_signal,M);%求最优信道均衡系数figure;plot(w);title('均衡器系数');ww = w'; %nTs2 = (0:one_Frame_points-1)*Ts;start_point = arrive_set-32-match_signal_points; %取第一帧的起始时刻one_Frame_data = receive_signal(start_point:(start_point+one_Frame_points-1)); %取一帧数据one_Frame_data_equ = conv(one_Frame_data,ww); %利用得到的均衡系数，当前帧数据进行滤波去除多途影响one_Frame_data_equ = one_Frame_data_equ(1:(end-M+1)); %去掉尾部M个数据figure;subplot(2,1,1);plot(nTs2,one_Frame_data);title('均衡前数据帧');subplot(2,1,2);plot(nTs2,one_Frame_data_equ);title('均衡后数据帧');%%% 使用FFT技术进行解调[out_Demodulation,energy_f] = BFSK_Demodulation_FFT(one_Frame_data,match_signal_points,Code_w_points,non_data_points,Code_num);[out_Demodulation_equ,energy_f_equ] = BFSK_Demodulation_FFT(one_Frame_data_equ,match_signal_points,Code_w_points,non_data_points,Code_num);figure;plot(energy_f(1,:),'r-*');hold on;plot(energy_f(2,:),'b-*');legend('f0','f1');string = ['均衡前各信号能量,信噪比SNR=',num2str(SNR)];title(string);figure;plot(energy_f_equ(1,:),'r-*');hold on;plot(energy_f_equ(2,:),'b-*');legend('f0','f1');string = ['均衡后各信号能量,信噪比SNR=',num2str(SNR)];title(string);[D_Code,error_flag] = Linear_Block_Decode(out_Demodulation,G_gener,H_check); %译码[D_Code_equ,error_flag_equ] = Linear_Block_Decode(out_Demodulation_equ,G_gener,H_check); %译码str1 = [',无错';',有错可被纠正 ';',有错不可被纠正'];for i = 1:Bytes_numstring = ['编码值',num2str(Code(i)),',均衡前译码值',num2str(D_Code(i)),str1(error_flag(i)+1,:),',均衡后译码值',num2str(D_Code_equ(i)),str1(error_flag_equ(i)+1,:),',信噪比SNR=',num2str(SNR),'dB\n'];fprintf(string);end

其中Com_Sonar_data_generation函数的内容如下（具体的模型介绍参见浅水声信道模型的建立）：

function receive_signal = Com_Sonar_data_generation(c,fs,f0,phase,BT,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num)%**********************************通信声呐水听器数据生成函数************************************%几点假设：%（1）基于射线声学理论%（2）几何衰减按球面波传播衰减规律衰减，不考虑吸收衰减%（3）考虑水面和水底的反射（浅海声信道？）%（4）考虑在高斯白噪声背景下%（5）整个空间声速分布均匀%***************************************输入参数说明******************************************% c 声速 Unit：m/s% fs 采样率 Unit：Hz% f0 信号频率（LFM时表示其中心频率） Unit：Hz% phase 调制相位Unit: 弧度% BT 时间增益（等于1时表示CW信号，大于1时表示LFM信号）% Tao信号脉宽Unit：s% Sample_time 采样时长Unit：s假设信号发射的时刻为零时刻% SNR信噪比 Unit：dB% H 水深Unit：m% H1 发射换能器水深 Unit：m% H2 接收换能器水深 Unit：m% D 接收与发射换能器水平距离 Unit：m% Re_coef_surf 水面反射系数 % Re_coef_bottom水底反射系数 存在半波损失% Reflex_num 考虑最大的反射次数%***************************************输出参数说明******************************************% receive_signal以信号发射为起始采样点的水听器接收信号%**************最后更新 1019%***更新内容%*首次编写%**************最后更新 1104%***更新内容%* (1)加入BT输入参数%* (2)可生成LFM信号%**************最后更新 0221%***更新内容%* (1)加入相位调制输入参数 phase%%%**********************************************************************************************************************H1b = H - H1; %发射换能器到水底的距离 Unit：mH2b = H - H2; %接收换能器到水底的距离 Unit：mTs = 1/fs;%采样时间间隔sample_num = fix(Sample_time*fs); %采样总点数nTs = (0:sample_num-1)/fs; %离散的采样时刻sample_num1 = fix(Tao*fs); %信号的采样点数nTs1 = (0:sample_num1-1)/fs; %信号的离散的采样时刻receive_signal = zeros(size(nTs));%用于存储水听器接收的信号d0 = sqrt((H1-H2)^2+D^2);%两个换能器的直线距离S0 = 1/d0;%直达波的声压幅值Noise_var = 10^(-SNR/20)*S0; %白噪声的方差%% 计算到达回波信号real_time = d0/c; %信号的实际到达时刻signal_start_time = fix(real_time*fs)+1; %信号第一个采样点的时刻delta_t = signal_start_time*Ts-real_time; %得到信号的第一个采样点和实际信号开始时刻点的差值(总是落后于信号开始点)delta_f = BT/Tao; %LFM信号的频率变化量if BT==1 %CW信号temp_signal = S0*sin(2*pi*f0*(nTs1+delta_t) + phase); %模拟采样elseif abs(BT)>1%LFM信号 正数表示频率由小变大，负数表示频率由大变小t = nTs1+delta_t - Tao/2; %标准表达式的形式k = delta_f/Tao; %频率变化率temp_signal = S0*sin(2*pi*f0*t+pi*k*t.^2); %模拟采样endreceive_signal(signal_start_time:(signal_start_time+sample_num1-1)) = temp_signal;receive_signal = receive_signal + Noise_var*randn(size(nTs));%加入噪声if Reflex_num>0 %考虑反射时for i = 1:Reflex_num %考虑i次反射的情况if mod(i,2) == 0%偶数次时%---------------------------------------------------第一次从水面反射时--------------------------------------------------- M = i*H+H1-H2; %中间系数k1 = sqrt(1+(D/M)^2); %中间系数d1 = k1*H1; %发射换能器距离水面第一个反射点的距离d2 = k1*H;%界面相邻两点间的距离d3 = k1*(H-H2);%水听器距离边界最后一个反射点最近的距离di = d1+d2*(i-1)+d3; %反射波总路程Si = 1/di*Re_coef_surf^(i/2)*Re_coef_bottom^(i/2); %反射波的声压幅值real_time = di/c; %信号的实际到达时刻signal_start_time = fix(real_time*fs)+1; %信号第一个采样点的时刻delta_t = signal_start_time*Ts-real_time; %得到信号的第一个采样点和实际信号开始时刻点的差值(总是落后于信号开始点)if BT==1 %CW信号temp_signal = Si*sin(2*pi*f0*(nTs1+delta_t) + phase); %模拟采样elseif abs(BT)>1%LFM信号 t = nTs1+delta_t - Tao/2; %标准表达式的形式temp_signal = Si*sin(2*pi*f0*t+pi*k*t.^2); %模拟采样endreceive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =...receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + temp_signal; %模拟采样%---------------------------------------------------第一次从水底反射时--------------------------------------------------- Mb = i*H+H1b-H2b; %中间系数k1b = sqrt(1+(D/Mb)^2); %中间系数d1b = k1b*H1b; %发射换能器距离水底第一个反射点的距离d2b = k1b*H; %界面相邻两点间的距离d3b = k1b*(H-H2b); %水听器距离边界最后一个反射点最近的距离dib = d1b+d2b*(i-1)+d3b; %反射波总路程Sib = 1/dib*Re_coef_surf^(i/2)*Re_coef_bottom^(i/2);%反射波的声压幅值real_time = dib/c;%信号的实际到达时刻signal_start_time = fix(real_time*fs)+1; %信号第一个采样点的时刻delta_t = signal_start_time*Ts-real_time; %得到信号的第一个采样点和实际信号开始时刻点的差值(总是落后于信号开始点)if BT==1 %CW信号temp_signal = Sib*sin(2*pi*f0*(nTs1+delta_t) + phase); %模拟采样elseif abs(BT)>1%LFM信号 t = nTs1+delta_t - Tao/2; %标准表达式的形式temp_signal = Sib*sin(2*pi*f0*t+pi*k*t.^2); %模拟采样endreceive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =...receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + temp_signal; %模拟采样else %奇数次时%---------------------------------------------------第一次从水面反射时--------------------------------------------------- M = (i-1)*H+H1+H2; %中间系数k1 = sqrt(1+(D/M)^2);%中间系数d1 = k1*H1;%发射换能器距离水面第一个反射点的距离d2 = k1*H;%界面相邻两点间的距离d3 = k1*H2;%水听器距离边界最后一个反射点最近的距离di = d1+d2*(i-1)+d3; %反射波总路程Si = 1/di*Re_coef_surf^((i-1)/2+1)*Re_coef_bottom^((i-1)/2); %反射波的声压幅值real_time = di/c; %信号的实际到达时刻signal_start_time = fix(real_time*fs)+1; %信号第一个采样点的时刻delta_t = signal_start_time*Ts-real_time; %得到信号的第一个采样点和实际信号开始时刻点的差值(总是落后于信号开始点)if BT==1 %CW信号temp_signal = Si*sin(2*pi*f0*(nTs1+delta_t) + phase); %模拟采样elseif abs(BT)>1%LFM信号 t = nTs1+delta_t - Tao/2; %标准表达式的形式temp_signal = Si*sin(2*pi*f0*t+pi*k*t.^2); %模拟采样endreceive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =...receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + temp_signal; %模拟采样%---------------------------------------------------第一次从水底反射时--------------------------------------------------- Mb = (i-1)*H+H1b+H2b;%中间系数k1b = sqrt(1+(D/Mb)^2); %中间系数d1b = k1b*H1b; %发射换能器距离水底第一个反射点的距离d2b = k1b*H; %界面相邻两点间的距离d3b = k1b*H2b; %水听器距离边界最后一个反射点最近的距离dib = d1b+d2b*(i-1)+d3b;%反射波总路程Sib = 1/di*Re_coef_surf^((i-1)/2)*Re_coef_bottom^((i-1)/2+1); %反射波的声压幅值real_time = dib/c; %信号的实际到达时刻signal_start_time = fix(real_time*fs)+1; %信号第一个采样点的时刻delta_t = signal_start_time*Ts-real_time; %得到信号的第一个采样点和实际信号开始时刻点的差值(总是落后于信号开始点)if BT==1 %CW信号temp_signal = Sib*sin(2*pi*f0*(nTs1+delta_t) + phase); %模拟采样elseif abs(BT)>1%LFM信号 t = nTs1+delta_t - Tao/2; %标准表达式的形式temp_signal = Sib*sin(2*pi*f0*t+pi*k*t.^2); %模拟采样endreceive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =...receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + temp_signal; %模拟采样 endendendreceive_signal = receive_signal/S0; %将直达波的幅值归一化

其中Linear_Block_Encode为线性分组码编码函数，具体内容如下：

function code = Linear_Block_Encode(Byte_code,G)%*************************(7,4,3)线性分组码编码子函数***********************%* 码字长度n=7 信息组长度k=4 分组码的最小距离d=3%***输入参数% Byte_code： 字节数据% G： 分组码生成矩阵%***输出参数% code：编码输出码字 0 1码流%***最后更新 1114%**更新内容 首次编写%说明：为实现一个字节数据传输，采用了两个(7,4,3)线性分组码%******************************************%%Bytes_num = length(Byte_code); %求字节个数code = zeros(1,14*Bytes_num); %生成编码结果存储矩阵offset = 0;for i = 1:Bytes_num%将一字节数据转换成 0 1 的码流形式msg = [ rem(Byte_code(i),2),...rem(fix(Byte_code(i)/2),2),...rem(fix(Byte_code(i)/4),2),...rem(fix(Byte_code(i)/8),2),...rem(fix(Byte_code(i)/16),2),...rem(fix(Byte_code(i)/32),2),...rem(fix(Byte_code(i)/64),2),...rem(fix(Byte_code(i)/128),2)];%计算编码输出code((offset+1):(offset+7)) = rem(msg(1:4)*G,2);%相乘后去2的余数 低4位编码code((offset+8):(offset+14)) = rem(msg(5:8)*G,2);%相乘后去2的余数 高4位编码offset = offset + 14;end

其中BFSK_Modulation为BFSK调制函数，具体内容如下：

function out = BFSK_Modulation(Code,signal_Pilot,sig0,non_data_points,Code_w_points)%*****************************BFSK调制子函数*****************************%***输入参数% Code 二进制编码后的码流，低位在前高位在后 LSB ~~~ MSB% signal_Pilot导频信号% sig0 代表码元信号 第一行代表码值‘0’，第二行代表码值‘1’ % non_data_points 导频信号与第一位数据间隔的点数% Code_w_points一个码元占用的点数%***输出参数% out 编码信号%***最后更新 0524%**更新内容 首次编写%******************************************%%code_num = size(Code,2); %计算码元个数out = signal_Pilot; %加入导频信号offset = non_data_points;%for i=1:code_numif Code(i) == 1%码元值为1时temp = sig0(2,:);out = out + [temp((end-offset+1):end),temp(1:(end-offset))];else %码元值为0时temp = sig0(1,:);out = out + [temp((end-offset+1):end),temp(1:(end-offset))];endoffset = offset + Code_w_points;%下一位码元的数据点偏移end

其中Get_Frame_Start_Time为信号到达时刻搜索函数，具体内容如下：

function arrive_set = Get_Frame_Start_Time(LPF_abs_xcorre_xy,sub_frame_points,threshold)%***************从匹配滤波包络信号中获取一帧数据的起始点**********************%****输入参数% LPF_abs_xcorre_xy：匹配滤波后的包络信号% sub_frame_points： 子帧的点数，也即导频信号点数% threshold： 相对阈值%****输出参数% arrive_set：信号的开始点，返回-1表示没有找到信号起始点即可能没有信号%***最后更新 1118%**更新内容 首次编写%******************************************sample_num = size(LPF_abs_xcorre_xy,2);%搜寻最大值i=1;while i<=(sample_num-sub_frame_points+1)temp_signal = LPF_abs_xcorre_xy(i:(i+sub_frame_points-1));%取出当前数据帧total_energy = sqrt(sum(temp_signal.^2));temp_signal = temp_signal/total_energy; %能量归一化[max_value,max_idx] = max(temp_signal); %找出最大值和最大值的坐标if max_value > threshold %最大值大于一定的值才能算有效的最大值arrive_set = max_idx + i - 1; %记录其在原始帧中的位置return; %找到一个后返回endi = i+fix(sub_frame_points/2); %下次的起点 endarrive_set = -1;%表示没有找到大于相对阈值的最大值

其中Channel_equalization为信号均衡函数，详细的理论介绍参考博客浅谈自适应滤波器，具体内容如下：

function w = Channel_equalization(signal,desired_signal,M)%*******************************信道均衡子函数******************************%****输入参数% signal： 待均衡信号% desired_signal： 期望信号% M： FIR滤波系数%****输出参数% w： FIR最优系数%***最后更新 1118%**更新内容 首次编写%**************************************************************************% signal = signal((M+1):end); % 去掉多途不重合的部分% desired_signal = desired_signal((M+1):end); % 去掉多途不重合的部分% % signal = signal(end:-1:1); %倒序 使序号小的表示时刻大的% desired_signal_points = size(desired_signal,2); %期望信号点数% X = zeros(M,desired_signal_points);%构造信号矢量矩阵% h = desired_signal(end:-1:1)'; %期望信号% for i=1:M%X(i,1:(end-i+1)) = signal(i:end);%构造第i行% end% % w = (X*X')\(X*h); %求最优系数%%%****************************常规的RLS算法*********************************% lamda = 1;%定义遗忘因子% I = eye(M); %生成对应的单位矩阵% c = 1; %小正数 保证矩阵P非奇异% % Signal_Len = size(signal,2); % % Signal_Len = Signal_Len - M; %一开始的M点不要% % Signal_Len = 5*M; %一开始的M点不要% % for i=1:Signal_Len%%输入数据%if i == 1 %如果是第一次进入% P_last = I/c;% w_last = zeros(M,1); %end%d = desired_signal(i+M); %输入新的期望信号 列向量 从第M个点开始 即去掉了多途重合部分%x = signal((M + i):-1:(i+1))';%输入新的信号矢量 列向量%%算法正体%K = (P_last * x)/(lamda + x'* P_last * x); %计算增益矢量%y = x'* w_last;%计算FIR滤波器输出%Eta = d - y; %计算估计的误差%w = w_last + K * Eta;%计算滤波器系数矢量%P = (I - K * x')* P_last/lamda;%计算误差相关矩阵%%变量更替%P_last = P;%w_last = w;%%滤波结果存储% %y_out(i) = y;% %Eta_out(i) = Eta;% %w_out(i,:) = w';% end%%%****************************横向快速RLS算法*********************************Signal_Len = size(signal,2); %获取信号长度% Signal_Len = fix(Signal_Len/2);lambda = 1;%定义遗忘因子delta = 0.01; for i = 1:Signal_Lenif i==1w_f_last = zeros(M,1);w_b_last = zeros(M,1);w_last = zeros(M,1);Phi_last = zeros(M,1);gamma_last = 1;xi_f_last = delta;xi_b_last = delta;x_N_1 = zeros(M+1,1);end%输入数据if i<= Mx_N_1(1:i) = signal(i:-1:1)'; %输入新的信号矢量else x_N_1 = signal(i:-1:(i-M))'; %输入新的信号矢量end d = desired_signal(i); %输入新的期望信号 列向量%算法主体e_f = x_N_1' * [1;-w_f_last]; %(1)epsilon_f = e_f * gamma_last; %(2)xi_f = lambda * xi_f_last + e_f * epsilon_f; %(3)w_f = w_f_last + Phi_last * epsilon_f; %(4)Phi_N_1 = [0;Phi_last] + e_f/(lambda * xi_f_last)*[1;-w_f_last]; %(5)Phi_N_1_N_1 = Phi_N_1(end);gamma_N_1 = (lambda * xi_f_last)/xi_f * gamma_last;%(6)e_b = lambda * xi_b_last * Phi_N_1_N_1; %(7)gamma = 1/(1/gamma_N_1 - Phi_N_1_N_1*e_b); %(8)epsilon_b = e_b * gamma; %(9)xi_b = lambda * xi_b_last + epsilon_b * e_b; %(10)Phi = Phi_N_1 - Phi_N_1_N_1 * [-w_b_last;1]; %(11)Phi = Phi(1:end-1); w_b = w_b_last + Phi * epsilon_b; %(12)x = x_N_1(1:M);y = w_last'* x;e = d - y; %(13)epsilon = e * gamma; %(14)w = w_last + Phi * epsilon; %(15)%变量更替xi_f_last = xi_f;w_f_last = w_f;gamma_last = gamma;xi_b_last = xi_b;Phi_last = Phi;w_b_last = w_b;w_last = w;end

其中BFSK_Demodulation_FFT为基于FFT的解调函数，具体内容图下：

function [out_Demodulation,energy_f] = BFSK_Demodulation_FFT(one_Frame_data,match_signal_points,Code_w_points,non_data_points,Code_num)%**************************基于FFT技术的BFSK解调函数*************************%***输入参数% one_Frame_data 提取后的数据帧的原始信号% match_signal_points 引导信号的点数% Code_w_points 一个码元占用的点数% non_data_points导引信号与数据信号间隔的点数% Code_num 码元个数%***输出参数% out_Demodulation解调 输出0 1 码流 二进制% energy_f 码元的个频率分量能量%***最后更新 0524%**更新内容 首次编写%******************************************%%offset = match_signal_points + non_data_points;out_Demodulation = zeros(1,Code_num);energy_f = zeros(2,Code_num); %当前码元各频率分量的能量for i=1:Code_numcode_fft = fft(one_Frame_data((offset+1):(offset+Code_w_points))); %求FFT energy_f(1,i) = abs(code_fft(27)); % '0'码的能量energy_f(2,i) = abs(code_fft(29)); % '1'码的能量if energy_f(1,i) > energy_f(2,i) %满足码值‘0’的条件out_Demodulation(i) = 0;elseout_Demodulation(i) = 1;endoffset = offset + Code_w_points;end

其中Linear_Block_Decode为线性分组码译码函数，具体内容如下：

function [D_code,error_flag] = Linear_Block_Decode(rec,G,H)%************************(7,4,3)线性分组码译码子函数**********************%* 码字长度n=7 信息组长度k=4 分组码的最小距离d=3%****输入参数% rec： 译码器接收的硬判决码字 14的整数倍个元素% G： 分组码生成矩阵4行7列矩阵% H： 分组码的校验矩阵 3行7列矩阵%****输出参数% D_code： 输出译码字节结果 正常范围是0~255 若有错则输出-1% error_flag： 错误标志 0：表示无错 1：有错但已经纠正 2：有错但无法纠正%***最后更新 1114%**更新内容 首次编写%***最后更新 0220%**更新内容 函数可一次自动译码处多个字节%******************************************%%num_rec = size(rec,2); %求总的 01 码流个数Bytes_num = fix(num_rec/14);%求需要解码的字节数D_code = zeros(1,Bytes_num);%存储输出的字节error_flag = zeros(1,Bytes_num); %初始化错误标志D_code_0_1 = zeros(1,8); %生成0 1译码码流输出offset = 0;for i = 1:Bytes_num[temp,error_flag(i)] = Block_decoder7_4_3(rec((offset+1):(offset+7)),G,H);%解码出低4位数据及错误标志if error_flag(i) == 2 %解码出错返回D_code(i) = -1; continue; else %否则正常输出D_code_0_1(1:4) = temp;end[temp1,error_flag(i)] = Block_decoder7_4_3(rec((offset+8):(offset+14)),G,H);%解码出高4位数据及错误标志if error_flag(i) == 2 %解码出错返回D_code(i) = -1; continue; else %否则正常输出D_code_0_1(5:8) = temp1;end%组合生成译码字节结果D_code(i) = D_code_0_1(1)*2^0 +... D_code_0_1(2)*2^1 +...D_code_0_1(3)*2^2 +...D_code_0_1(4)*2^3 +...D_code_0_1(5)*2^4 +...D_code_0_1(6)*2^5 +...D_code_0_1(7)*2^6 +...D_code_0_1(8)*2^7;offset = offset + 14; %改变码流偏移量end

运行后结果如下：

窗口打印输出如下：

如果觉得《大话水声通信技术---（BFSK仿真）》对你有帮助，请点赞、收藏，并留下你的观点哦！