cotas_carrier_recovery.m

美国弗州大学的讲义,用cotas环实现载波恢复

% AM信号生成部分
Fs = 16000; fm = 200; t_span = 0.5;
t = 0:1/Fs:t_span;
Am = 1; m = Am*sin(2*pi*fm*t);
theta = 1;Ac = 1;fc = 2000;
input = Ac*sin(2*pi*fc*t+theta).*m; %Input is a modulated signal
% Hilbert滤波器设计部分
order = 32;
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)
% 生成解析信号部分
hilbert_output = filter(b, 1, input);
%groupdelay should be accounted for generating a analytic signal
len = length(input);
analytic = input(1:len-groupdelay) + j*hilbert_output(groupdelay+1:len);
% 环路滤波器设计部分
BL = 10; % Arbitrarily set
zeta = 0.707; % Set in Specs
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
ratio=alpha/beta
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)]);
% 同步部分
%Initialization
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;
    phase(I) = exp(-i*phi(I));
    c1(I) = real(analytic(I)*phase(I));
    c2(I) = imag(analytic(I)*phase(I));
    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;
    temp_out3=2*pi*(fc-5)/Fs+phi(I)+temp_out2;
    temp_pre_out1=temp_out1;
end
x = 1:len-groupdelay;
plot(x,q);