costas环设计与使用（草稿）

costas环的结构

costas环的参数设计

costas环主要就是针对$c_1,c_2$的计算。
参数的计算公式，
$$\begin{gathered} c_1=\frac{2\xi w_{n}T}{K} \\ {c}_2=\frac{(w_{n}T{)}^2}{K} \end{gathered}$$
其中，
$$w_n=\frac{8\xi B_L}{4{\xi }^2+1}$$
式中，
$$\begin{gathered} \xi =0.707(或\xi =\frac{1}{\sqrt{2}}) \\ {B}_L\le 0.1R_b \\ T=NCO更新周期 \\ {R}_b=信息符号率 \end{gathered}$$

costas环设计例子

例1

用costas解调AM信号或BPSK信号。这里用了希尔伯特变换构成解析信号办法。

% 产生调幅信号
Fs = 16000; % 采样率
fm = 200;  % 调幅信号的频率
t_span = 0.5; % 仿真时间

t = 0:1/Fs:t_span;
Am = 1; % 信号幅度
m = Am*cos(2*pi*fm*t);
theta = 1; % 相位
Ac = 1; % 载波幅度
fc = 2e3; % 载波频率
frc = 1.9e3; % 接收端载波初始频率
input = Ac*sin(2*pi*fc*t+theta).*m; % 调幅调制信号

% 希尔伯特滤波器的设计
order = 32;
N = order + 1;
f = 0.1:0.001:0.901;
a = ones(length(f));
b = remez(order, f, a, 'hilbert');

[gd,f] = grpdelay(b,1,512,'whole',Fs);
groupdelay=gd(10);

% 解析信号生成
len = length(input);
hilbert_output = filter(b, 1, input);
analytic = input(1:len-groupdelay) + 1i*hilbert_output(groupdelay+1:len);

% 环路滤波器的参数计算
BL = 1000; % Arbitrarily set
zeta = 0.707;
k1 = mean(abs((Ac*m).^2));
alpha = 2*zeta*2*BL/(zeta+1/4/zeta)/Fs/k1;
beta = ((2*BL/(zeta+1/4/zeta)/Fs))^2/k1;

% 环路滤波器的频率响应
Hnum = [ k1*(beta+alpha) -k1*alpha];
Hden = [1 -(2-k1*(alpha+beta)) (1-k1*alpha)];
LF = filter(Hnum, Hden, [zeros(1,100) ones(1,2000)]);

% 初始化
out = [];
phi = [];
phi(1) = 0;
temp_out1=0;
temp_pre_out1=0;
temp_out2=0;
temp_out3=0;

%Simulation
for I=1:len-groupdelay
    phi(I)= temp_out3; % NCO的相位
    phase(I) = exp(-1i*phi(I)); % 求cos(theta)-j*sin(theta)
    c1(I) = real(analytic(I)*phase(I));
    c2(I) = imag(analytic(I)*phase(I));
    
    % 鉴相器,这里也可以选择乘法鉴相器
    out(I)= c1(I);% for bpsk out(I)=sign(c1(I));
    q(I) = sign(c1(I))*c2(I); 
    
    % 环路滤波器
    temp_out1=temp_pre_out1+q(I)*beta;
    temp_out2=alpha*q(I)+ temp_out1;
    
    % 模拟NCO的相位累加
    temp_out3=2*pi*frc/Fs+phi(I)+temp_out2;
    temp_pre_out1=temp_out1;
end

plot(t,m); % 绘制原始信号
hold on;
plot(t,input); % 绘制调制后信号
plot(t(1:len-groupdelay),out); % 绘制解调后信号
legend('AM','AM-modulated','AM-demodulated');
hold off;

figure;
% phase = 2*pi*fc*t
% delta_phase = 2*pi*fc*(t1-t0)
% fc = delta_phase/(2*pi)*Fs

delta_phase = phi(2:end)-phi(1:end-1);
real_freq = delta_phase/(2*pi)*Fs;
plot(t(1:length(real_freq)),real_freq);
xlabel('时间(s)');
ylabel('实时频率(Hz)');

例2