直接来干货！同学们一看就知道这个代码的分量了！！！

%%%%%%%%%%%%%%%         MSK相干解调符号同步仿真程序             %%%%%%%%%%%%%

%%%%%%%%%%%%      File: mskwaveform_search1andtrack_rel1.m        %%%%%%%%%

%%%%%%%%     date: 2009-07-08  author: haige nanjing research         %%%%%

%%%                           程序说明

% 本程序为mskwaveform_search1andtrack_sim1的升级版，改变了跟踪时的误码率计算方式，

% 使得跟踪不正确时可以不统计误码率指标。

% 程序采用时域捕获算法完成初始频偏估计和位同步估计，采用DDPLL和ad hoc feedback scheme完成载波恢复和位同步。

% 程序在位同步完成后进行载波恢复，载波恢复针对每个符号进行。

% 位跟踪环路完全采用gardner的环路结构，定时误差鉴别器采用内插点进行计算，环路滤波器的参数采用了理想二阶环的参数。

% 本地8倍的符号率，内抽出4个点。本算法要求在一个符号内至少抽取出4个采样点。

% 在程序在msk_timingtrk1_DDPLL_coh_rel3和msk_search1_rel2合并的基础上进行波形仿真。

% 本程序采用规范帧头和BTR码，输出软信息结果。threshold已经赋值。

%%%                           仿真环境

% 信号源设置：含有帧标志，帧标志和信息都采用MSK调制。

% 采样速率：符号率的整数倍，奈奎斯特采样

% 采样偏差：0    载波偏差：可变    符号率: 2.4KHz

% BTR码采用1100序列，总共128比特，帧头42比特，000001010001110010101001100000101110010011。

% mil-std-188-182标准中600bps时，引导序列长284比特，1200bps时，引导序列长358比特，

% 由于采用soqpsk调制，1个符号2比特。

clear all;

close all;

%%*****************            程序主体            ***********************%%

format long;

%****************           commom  parameter            ******************

%%%************** 系统参数

sr = 2400;          % Symbol rate

ml = 1;             % Number of modulation levels

br = sr.*ml;        % Bit rate

IPOINT = 8;         % Number of oversamples

%%%************** 锁相环参数

damp =0.707;

Bp = 100;

tp = 1/sr;

Kp = 1;

Wp = 2*Bp/(damp+1/(4*damp));

alpha1 = 8*damp*Wp*tp/(Kp*(4+4*damp*Wp*tp+(Wp*tp)^2));

beta1  = 4*(Wp*tp)^2/(Kp*(4+4*damp*Wp*tp+(Wp*tp)^2));

%%%************** 位跟踪环参数

damp = 0.707;

Bd = 0.5;

td = 1/sr;

Kd = 1;

Wd = 2*Bd/(damp+1/(4*damp));

alpha2 = 2*damp*Wd*td;

beta2  = Wd^2*td^2;

%%%**************  loop  parameter

nloop = 1e2;             % Number of simulation loops

%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%              发射机部分            %%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%

%******************      BTR code and frame  start    *********************

frame_head_sample = [  1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 ];  % 32 bits

% frame_head_sample = [ 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 ];

% frame_head_sample = randint(1,32);  % 32 bits

frame_head = [frame_head_sample frame_head_sample frame_head_sample frame_head_sample];   % BTR code  all in all for 128 bits

detect_flag1 = [ 0 0 0 0 0 1 0 1 0 0 0 1 1 1 0 0 1 0 1 0 1 0 0 1 1 0 0 0 0 0 1 0 1 1 1 0 0 1 0 0 1 1];  % frame of start

%**********************         data  source        ***********************

data = randint(1,5000);

%***********************         creat frame        ***********************

frame = [randint(1,20) frame_head detect_flag1 data]; % 帧头前加装无用数据

len_frame = length(frame)-20;  % Number of symbols that simulates in each loop

%**************************     MSK modulation    *************************

data_source = 2*frame - 1;

for k=1:IPOINT

data_sample(k:IPOINT:k+IPOINT*(length(data_source)-1)) = data_source;

end

phase = zeros(1,length(data_sample)+1);

ich = zeros(1,length(data_sample)+1);

qch = zeros(1,length(data_sample)+1);

for ii=2:length(data_sample)+1

phase(ii) = phase(ii-1) + pi/2*data_sample(ii-1)./IPOINT;

end

ich = cos(phase);

qch = sin(phase);

%%%************** Add  timing offset

toffset = 5;

% 采样偏差造成的移位影响 (16*2.4e3+10)*4000/(16*2.4e3) = 4001

% (8*2.4e3+10)*4000/(8*2.4e3) = 4002  (8*2.4e3+5)*4000/(8*2.4e3) = 4001

ideal_samplefre = IPOINT*sr;

samplefre = ideal_samplefre + toffset;

offset_start = 1.0;  %%% 采样偏差的起始时刻

time1 = ideal_samplefre/samplefre.*(offset_start:length(ich)+offset_start-1);

%%%  生成以IPOINT*sr采样频率采样的数据

ich_sample = interp1((1:length(ich)),ich,time1,'spline');

qch_sample = interp1((1:length(qch)),qch,time1,'spline');

%%%************** Add carrier frequency offset  (assuming perfect tracking)

freoffset = 120;      % 频偏

pha_offset = 0;      % 相偏

%%%************** Generate complex signal

time = (0:length(ich_sample)-1)/IPOINT/sr;

% 基带信号用复信号的形式来表达

tra_sig = (ich_sample + j*qch_sample).*exp(j*(2*pi*freoffset.*time+pha_offset));

fprintf('%%设置起始位置偏差%5.2f 时钟偏差%dHz TLL噪声带宽%5.2f 频偏%5.2fHz PLL噪声带宽%5.2f 采样频率%d倍符号率 %d个符号循环%d次\n',(offset_start-1),toffset,Bd,freoffset,Bp,IPOINT,len_frame,nloop);

% fprintf('%%******************* 跟踪仿真结果 ******************* \n');

fprintf('%%\t信噪比\t误比特率\n');

当年的程序注释写的非常棒，因此就没有写算法文档。下一堂课给出接收机部分的内容！！！

20190217 完成初稿！