用matlab的2018b版本写一个gmsk相干解调

以下是一个简单的GMSK相干解调的Matlab代码示例：

%% 初始化参数
fc = 100; %载波频率
fs = 1000; %采样频率
Tb = 0.1; %比特周期
N = 1000; %采样点数
kf = 5; %调频斜率
Eb = 1; %比特能量
SNR = 10; %信噪比
phi_0 = 0; %初始相位
phi_hat = 0; %相位误差估计
theta_hat = 0; %频率误差估计
r = rand(1,N)>0.5; %生成随机比特序列

%% 生成调制信号
t = linspace(0,Tb,N);
f = kf/Tb; %调频频率
phi = pi*r; %生成不归零码
s = zeros(1,N);
for i = 2:N
    s(i) = s(i-1) + 2*pi*f*t(i-1) + phi(i); %相位连续积分
end
s = cos(s);

%% 加入高斯白噪声
sigma = sqrt(Eb/(2*10^(SNR/10))); %计算噪声标准差
noise = sigma*randn(1,N); %生成噪声信号
x = s + noise; %加噪声

%% 相干解调
I = x.*cos(2*pi*fc*t+phi_0); %乘以本地载波并取实部，得到I路信号
Q = -x.*sin(2*pi*fc*t+phi_0); %乘以本地载波并取虚部，得到Q路信号
I_filt = lowpass(I, fc/10, fs); %低通滤波
Q_filt = lowpass(Q, fc/10, fs); %低通滤波
d_phi = atan2(Q_filt, I_filt); %相位差
phi_hat = phi_hat + d_phi; %相位误差累积
theta_hat = mean(diff(phi_hat)); %频率误差估计

%% 绘制结果
figure;
subplot(2,1,1);
plot(t, s);
title('调制信号');
subplot(2,1,2);
plot(t, x);
title('接收信号');
figure;
plot(t, phi_hat);
title('相位误差');