GMSK调制 MATLAB代码

%span 符号数
%sps 每个符号样本数
%BT 时间积
%Fs 采样频率
%Fb 码元速率
Nbit = 1000 ;
BT = 0.5 ;
span = 8 ;
sps = 130 ;
Fb = 50e3 ;
Fc = 4 * Fb ;
Fs = sps * Fb ;
K = ceil(Fs / 2 / Fc - 1) ;% 不重要，这里是采用了多载波发送，K表示可用载波数量
fcmat = (1 : K)' .* repmat(Fc, K, 1) ;
% CRC校验码
crc_polynominal = [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1] ;%生成多项式16次
crcgenerator = comm.CRCGenerator(crc_polynominal);  
msgmat = zeros(K, Nbit + length(crc_polynominal) - 1) ; % 信息长度
for i = 1:size(msgmat, 1)
    msgmat(i, :) = crcgenerator(randi([0, 1], Nbit, 1))' ; % 一个长度为1000+16
end
modwave = zeros(1, sps * size(msgmat, 2)) ;
msg = msgmat(1,:);
Len = span+1 ; % symbol time span after Gaussian pulse shaping
%% 高斯滤波
csps = 200 * sps ; % high sample rate for discrete integral
% 升采样
gaussfilter = gaussdesign(BT, span, csps) ; 
rectbasewav = rectpulse([1], csps) ; % 阶跃函数
gt = conv(gaussfilter, rectbasewav) ;% 得到gt
deltaphi = zeros(Len, sps) ;
for i = 1:Len  % 一个信号可以影响到四个码元，长度为4*2+1 = 9
    for j = 1:sps       
        deltaphi(Len+1-i, j) = sum(gt(1+(i-1)*csps : 1+(i-1)*csps+j*csps/sps-1)) / sum(gt) * pi/2 ;
    end 
end  % 求出 gt的积分，长度为 9 T

%% MSK调制
bimsg = 2*msg - 1 ; % 转成双极性
padbimsg = [zeros(1, (Len-1)/2), bimsg, zeros(1, (Len-1)/2)] ;  % 前后增加4个长度
phitrace = zeros(1,length(msg)) ;
lastphi = 0 ; 
for k = 1:length(msg)
    kthphi = lastphi + padbimsg(k : k+Len-1) * deltaphi ;
    lastphi = kthphi(end) ;
    phitrace = [phitrace kthphi] ; 
end
t = (0:length(msg)*sps-1)/Fs ;
modwave = (cos(phitrace) + 1i*sin(phitrace)) .* exp(1i*2*pi*Fc*t) ;