Matlab实现科斯塔斯环，用于跟踪同步载波，请提供代码

科斯塔斯环是一种数字信号处理技术，用于同步解调数字通信信号。以下是Matlab实现科斯塔斯环的代码示例：

%% 参数设置
fc = 10e3;           % 载波频率
fs = 100e3;          % 采样率
Rs = 10e3;           % 符号速率
N = 10000;           % 仿真时长
EbN0 = 10;           % 信噪比

%% 生成调制信号
t = (0:N-1)'/fs;
data = sign(randn(N,1));
I = cos(2*pi*fc*t);
Q = sin(2*pi*fc*t);
s = I.*data + Q.*data;

%% 添加噪声
SNR = EbN0 + 10*log10(Rs/fs);
sigma = sqrt(1/(2*10^(SNR/10)));
noise = sigma*(randn(N,1)+1i*randn(N,1));
r = s + noise;

%% 科斯塔斯环解调
theta_hat = zeros(N,1);
I_acc = 0;
Q_acc = 0;
for n = 1:N
    I_acc = I_acc + r(n)*cos(theta_hat(n));
    Q_acc = Q_acc - r(n)*sin(theta_hat(n));
    theta_hat(n+1) = theta_hat(n) + 2*pi*(fc/Rs) + atan2(Q_acc, I_acc);
end

%% 绘制结果
figure;
subplot(2,1,1);
plot(t, real(s));
hold on;
plot(t, real(r));
xlabel('时间/s');
ylabel('幅值');
legend('调制信号', '接收信号');

subplot(2,1,2);
plot(t, mod(theta_hat,2*pi));
xlabel('时间/s');
ylabel('相位偏移');

代码中，首先设置了一些参数，包括载波频率fc、采样率fs、符号速率Rs、仿真时长N和信噪比EbN0。然后，生成调制信号，添加高斯噪声，并进行科斯塔斯环解调。最后，绘制了调制信号、接收信号和相位偏移随时间的变化。

注意，在实际应用中，需要根据具体的信号和系统参数进行调整。