加载中...

定时同步学习及仿真

定时同步理论学习

基本介绍

什么是定时同步?

定时同步也叫符号同步、位同步、码元同步,只在数字通信中需要。

为什么需要定时同步?

数字通信中使用几个采样点来表示一个码元,在接收端只需检测一个采样点即可恢复该码元。在实际通信系统中,由于收发系统间存在传输时延、收发晶振间存在固有频差等问题,导致接收端无法在每个码元对应的最佳采样点处精确采样,从而引入定时误差,使得输出信噪比减小。

星座图比对

定时同步对数字信号在传输过程中产生的定时偏差进行纠正,使得信号在最佳采样点处采样。

定时同步的实现方法

3种方法,同源时钟法、数据辅助法、自同步法;

  • 同源时钟法

    收发时钟不同源,导致收发之间存在偏差,该方法将收发时钟同步到一个主时钟上,主时钟给出一个确定的信号,接收端通过计算与主时钟之间的时钟延时来确定定时误差,从而实现定时同步。
    该方法多用于低频通信中,目前的高速无线通信并中不常用。

  • 数据辅助法

    发送端在发送信号的同时发送一些特定的辅助数据,接收端利用这些辅助数据提取出定时信息从而完成定时同步。该方法虽然简单易行,但是会占用发送机的部分功率和信道的部分带宽,降低传输效率。

  • 自同步法

    接收端直接从接收到的数据中通过参数估计等方法提取定时信号以实现同步。该方法实现的计算复杂度相对较高,但是计算的定时误差相比前两种更加准确,应用更加广泛,成为目前无线通信的首要选择。
    自同步法的实现分为开环同步和闭环同步两种。
    闭环定时同步环路闭环同步法将本地时钟与接收到的信号进行比较,通过定时误差估计以及定时误差校正等环节使本地时钟与接收到的信号同步,从而得到正确的采样信号以完成定时同步。目前应用较广泛.

闭环定时同步环路

反馈型全数字定时同步环路

优点: · 同步精度高 · 适用于长时间连续数字通信系统
缺点: **· 需要一定的收敛时间 · **结构相对复杂

综合考虑,一般采用反馈型全数字定时同步环路

Gardner算法

主要原理:不改变采样时钟的情况下,通过插值算法调整接收信号,使信号达到最佳采样点,从而实现定时同步。
对接收到的信号按照本地时钟周期 Ts 进行采样,将采样后的信号送入插值滤波器中,内插出输出的插值点。定时误差检测器采用 Gardner 算法对输出信号进行定时误差估计,将得到的估计值送入环路滤波器中进行平滑处理,再将得到的定时误差送入数控振荡器(NCO),最后将计算出的整数采样时刻 mk 和分数间隔 µk 反馈回插值滤波器,从而完成定时同步。将该过程一直迭代直至定时误差为 0 或固定值,此时数据已经几乎不存在定时误差。

  1. 定时误差检测器

    插值滤波器求得内插输出点后,通过$Gardner$误差检测算法计算定时误差,以反馈进行纠正

    🏀Gardner定时误差检测算法不受频偏影响

    • 推导过程

      载波频偏并不影响Gardner算法估计出来的定时误差。因此,Gardner算法可以用在载波同步之前。

  2. 环路滤波器

    作用:滤除输出的高频分量和噪声,维持环路稳定性,提高系统的抗干扰能力,使输出更加平滑

  3. 定时控制器

    定时控制器由数控振荡器(NCO)和分数间隔计算器组成,通过对定时误差信号进行计算,用来控制插值运算。

  4. 差值滤波器

    滤波器设计

    基于Farrow结构的分段抛物线插值

    基于Farrow结构的立方插值

MATLAB仿真

仿真要求

仿真测试

BPSK系统代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
clc; clear; close all;
%%%%%%%%%在BPSK收发系统的基础上,接收端加上定时误差,%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% 用resample加 ,然后用gardner 算法进行同步,判决算误码率%%%%%%%%%%%%%%%%%
symbol_rate = 1000; %符号速率
sps = 4; %每个符号4个采样点
fs = sps * symbol_rate; % 采样率
span = 6;
fc = 4 * fs; % 载波频率
alfs = 0.35; % 滚降系数
data_len = 10; %数据长度
is_plot = 1; %是否打印输出

data = randi([0 1], 1, data_len); %生成0,1序列

s = (data - 1/2) * 2; %映射
data_up = upsample(s, sps); %上采样

fir_rcos = rcosdesign(alfs, span, sps, "sqrt"); %根升余弦
data_bx = conv(fir_rcos, data_up); %通过滤波器
data_bx_valid = data_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);

t = 0:length(data_bx_valid) - 1;
carrier = cos(2 * pi * fc / fs * t); %载波
s_up = data_bx_valid .* carrier; %上变频

r_n = awgn(s_up, 10, 'measured'); %加高斯白噪声

r_down = r_n .* carrier; %下变频

r_bx = conv(fir_rcos, r_down); %匹配滤波器接收
r_bx_valid = r_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);

r_data = r_bx_valid(1:sps:end) > 0;

plot(r_data - data);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
clc; clear; close all;
%%%%%%%%%在BPSK收发系统的基础上,接收端加上定时误差,%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% 用resample加 ,然后用gardner 算法进行同步,判决算误码率%%%%%%%%%%%%%%%%%
symbol_rate = 1000; %符号速率
sps = 4; %每个符号4个采样点
fs = sps * symbol_rate; % 采样率
span = 6;
fc = 4 * fs; % 载波频率
alfs = 0.35; % 滚降系数
N = 1e6; %数据长度
is_plot = 1; %是否打印输出
%data = repmat([0,1], 1,N/2);
data = randi([0 1], 1, N); %生成0,1序列

s = (data - 1/2) * 2; %映射
data_up = upsample(s, sps); %上采样

fir_rcos = rcosdesign(alfs, span, sps, "sqrt"); %根升余弦
data_bx = conv(fir_rcos, data_up); %通过滤波器
data_bx_valid = data_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);

t = 0:length(data_bx_valid) - 1;
carrier = cos(2 * pi * fc / fs * t); %载波
s_up = data_bx_valid .* carrier; %上变频

c1 = 5.41 * 10^(-3); c2 = 3.82 * 10^(-6); %环路滤波器系数

%ii = 1; %用来表示Ts的时间序号,指示n,n_temp,nco,

errbit = zeros(1, 10);
snr = zeros(1, 10);
i = 1;
for ebn0 = 1:10
snr(i) = ebn0 + 10 * log10(2) - 10 * log10(sps); %SNR和EbN0的转化
r_n = awgn(s_up, snr(i), 'measured'); %加高斯白噪声

r_down = r_n .* carrier; %下变频

r_bx = conv(fir_rcos, r_down); %匹配滤波器接收
r_bx_valid = r_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);
len_r_bx_valid = length(r_bx_valid);

yI_decision = (r_bx_valid(1:sps:end)) > 0;
errbit(i) = sum(abs(yI_decision-data))/N;
i = i +1;
end

if is_plot
% figure;
% subplot(2, 1, 1);
% stem(data, 'filled');
% subplot(2, 1, 2);
% stem(yI_decision, 'filled');
figure;
ebn0 = 1:10;
ber_theory = 0.5 * erfc(sqrt(power(10, ebn0 / 10))); %理论误码率
semilogy(ebn0, ber_theory, '-or');
grid on; hold on;

semilogy(ebn0, errbit, '-sb');
grid on; hold on;
xlabel('EbN0(dB)'); ylabel('BER')
title('BPSK误码率曲线图');
legend('理论误码率曲线', '仿真误码率曲线');
%plot(w); title('环路滤波器输出');

end

定时同步(Gardner算法)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
clc; clear; close all;
%%%%%%%%%在BPSK收发系统的基础上,接收端加上定时误差,%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% 用resample加 ,然后用gardner 算法进行同步,判决算误码率%%%%%%%%%%%%%%%%%
symbol_rate = 1000; %符号速率
sps = 4; %每个符号4个采样点
fs = sps * symbol_rate; % 采样率
span = 6;
fc = 4 * fs; % 载波频率
alfs = 0.35; % 滚降系数
N = 50000; %数据长度
is_plot = 1; %是否打印输出
%data = repmat([0,1], 1,N/2);
data = randi([0 1], 1, N); %生成0,1序列

s = (data - 1/2) * 2; %映射
data_up = upsample(s, sps); %上采样

fir_rcos = rcosdesign(alfs, span, sps, "sqrt"); %根升余弦
data_bx = conv(fir_rcos, data_up); %通过滤波器
data_bx_valid = data_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);

t = 0:length(data_bx_valid) - 1;
carrier = cos(2 * pi * fc / fs * t); %载波
s_up = data_bx_valid .* carrier; %上变频

errbit = zeros(1, 10);
snr = zeros(1, 10);

i = 1;
for ebn0 = 1:10
snr(i) = ebn0 + 10 * log10(2) - 10 * log10(sps); %SNR和EbN0的转化
r_n = awgn(s_up, snr(i), 'measured'); %加高斯白噪声

r_down = r_n .* carrier; %下变频

r_bx = conv(fir_rcos, r_down); %匹配滤波器接收
r_bx_valid = r_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);


%接收端加上定时误差
P = 40000;
Q = 40002;
r_data_er = resample(r_bx_valid, P, Q); %输入序列x在P/Q乘以原始采样速率下进行重新采样
len_r_data_er = length(r_data_er);
%c1 = 5.41 * 10^(-3); c2 = 3.82 * 10^(-6); %环路滤波器系数
c1 = 1/2^11; c2 = 1/2^23;
%gardner环参数设置
w = [0.5, zeros(1, N - 1)]; %环路滤波器输出寄存器,初值设为0.5
n = [0.7 zeros(1, N - 1)]; %NCO寄存器,初值设为0.7

u = [0.6, zeros(1, len_r_data_er - 1)]; %NCO输出的定时分数间隔寄存器,初值设为0.6
time_error = zeros(1, N); %Gardner提取的时钟误差寄存器
yI = zeros(1, 2 * N); %I路内插之后的输出数据
%yI数组长度是u和te的两倍,其原因在于yI(n-1/2),matlab中仿真不方便,所以取一半的数
%ii = 1; %用来表示Ts的时间序号,指示n,n_temp,nco,
k = 1; %用来表示Ti时间序号,指示u,yI,yQ
ms = 1; %用来指示T的时间序号,用来指示a,b以及w

for ii = 2:(len_r_data_er - 3)
n(ii + 1) = n(ii) - w(ms);

if (n(ii + 1) < 0)
u(k) = n(ii) / w(ms);
n(ii + 1) = mod(n(ii + 1), 1);

%内插滤波器
F1 = 0.5 * r_data_er(ii + 2) - 0.5 * r_data_er(ii + 1) - 0.5 * r_data_er(ii) + 0.5 * r_data_er(ii - 1);
F2 = 1.5 * r_data_er(ii + 1) - 0.5 * r_data_er(ii + 2) - 0.5 * r_data_er(ii) - 0.5 * r_data_er(ii - 1);
F3 = r_data_er(ii);
yI(k) = (F1 * u(k) + F2) * u(k) + F3;

% 时钟误差提取模块
if (mod(k, 2) == 1)
%每一个数据符号计算一次误差
if (k > 2)
time_error(ms) = yI(k - 1) * (yI(k) - yI(k - 2));
else
time_error(ms) = 0;
end

if (ms > 1)
w(ms + 1) = w(ms) + c1 * (time_error(ms) - time_error(ms - 1)) + c2 * time_error(ms);
else
w(ms + 1) = w(ms) + c1 * time_error(ms) + c2 * time_error(ms);
end

ms = ms +1;
end
k = k + 1;
end
end
%抽样判决
%yI_decision = (sign(yI(1:2:end)) + 1) / 2;
yI_decision = (yI(1:2:end)) > 0;
errbit(i) = sum(abs(yI_decision(301:N) - data(301:N))) / (N-300);
i = i +1;
end

if is_plot
% figure;
% subplot(2, 1, 1);
% stem(data, 'filled');
% subplot(2, 1, 2);
% stem(yI_decision, 'filled');
figure;
ebn0 = 1:10;
ber_theory = 0.5 * erfc(sqrt(power(10, ebn0 / 10))); %理论误码率
semilogy(ebn0, ber_theory, '-or');
grid on; hold on;
semilogy(ebn0, errbit, '-sb');
grid on; hold on;
xlabel('EbN0(dB)'); ylabel('BER')
title('BPSK误码率曲线图');
legend('理论误码率曲线', '仿真误码率曲线');
% figure;
% plot(w); title('环路滤波器输出');
figure;
plot(n);

end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
clc; clear; close all;
%%%%%%%%%在BPSK收发系统的基础上,接收端加上定时误差,%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% 用resample加 ,然后用gardner 算法进行同步,判决算误码率%%%%%%%%%%%%%%%%%
symbol_rate = 1000; %符号速率
sps = 4; %每个符号4个采样点
fs = sps * symbol_rate; % 采样率
span = 6;
fc = 4 * fs; % 载波频率
alfs = 0.35; % 滚降系数
N = 50000; %数据长度
is_plot = 1; %是否打印输出
%data = repmat([0,1], 1,N/2);
data = randi([0 1], 1, N); %生成0,1序列

s = (data - 1/2) * 2; %映射
data_up = upsample(s, sps); %上采样

fir_rcos = rcosdesign(alfs, span, sps, "sqrt"); %根升余弦
data_bx = conv(fir_rcos, data_up); %通过滤波器
data_bx_valid = data_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);

t = 0:length(data_bx_valid) - 1;
carrier = cos(2 * pi * fc / fs * t); %载波
s_up = data_bx_valid .* carrier; %上变频

errbit = zeros(1, 10);
snr = zeros(1, 10);
i = 1;

for ebn0 = 1:10
snr(i) = ebn0 + 10 * log10(2) - 10 * log10(sps); %SNR和EbN0的转化
r_n = awgn(s_up, snr(i), 'measured'); %加高斯白噪声

r_down = r_n .* carrier; %下变频

r_bx = conv(fir_rcos, r_down); %匹配滤波器接收
r_bx_valid = r_bx(1, (1 + length(fir_rcos)) / 2:end - (length(fir_rcos) - 1) / 2);
%r_bx_valid_down = r_bx_valid(1:sps:end);不需要下采样
%接收端加上定时误差
P = 40000;
Q = 40002;
r_data_er = resample(r_bx_valid, P, Q); %输入序列x在P/Q乘以原始采样速率下进行重新采样
len_r_data_er = length(r_data_er);
%c1 = 5.41 * 10^(-3); c2 = 3.82 * 10^(-6); %环路滤波器系数
c1 = 1/2^11; c2 = 1/2^23;
%gardner环参数设置
w = [0.25, zeros(1, N - 1)]; %环路滤波器输出寄存器,初值设为0.5
n = [0.7, zeros(1, N - 1)]; %NCO寄存器,初值设为0.7

u = [0.6, zeros(1, len_r_data_er - 1)]; %NCO输出的定时分数间隔寄存器,初值设为0.6
time_error = zeros(1, N); %Gardner提取的时钟误差寄存器
yI = zeros(1, N); %I路内插之后的输出数据
%ii = 1; %用来表示Ts的时间序号,指示n,n_temp,nco,(相当于mk)
ms = 1; %用来指示T的时间序号,用来指示a,b以及w
k = 1; %用来表示Ti时间序号,指示u,yI,yQ
m = zeros(1, N);

for ii = 2:(len_r_data_er - 3)
n(ii + 1) = n(ii) - w(ms);

if (n(ii + 1) < 0) %小于零相当于溢出,此时ii为溢出最后采样点的索引mk
n(ii + 1) = n(ii + 1) + 1;
u(k) = n(ii) / w(ms);
m(k) = ii;

%内插滤波器
F1 = 0.5 * r_data_er(m(k) + 2) - 0.5 * r_data_er(m(k) + 1) - 0.5 * r_data_er(m(k)) + 0.5 * r_data_er(m(k) - 1);
F2 = 1.5 * r_data_er(m(k) + 1) - 0.5 * r_data_er(m(k) + 2) - 0.5 * r_data_er(m(k)) - 0.5 * r_data_er(m(k) - 1);
F3 = r_data_er(m(k));
yI(k) = (F1 * u(k) + F2) * u(k) + F3;

if (k > 1)
mk_mid = floor((m(k) + m(k - 1) + u(k) + u(k - 1)) / 2);
F1 = 0.5 * r_data_er(mk_mid + 2) - 0.5 * r_data_er(mk_mid + 1) - 0.5 * r_data_er(mk_mid) + 0.5 * r_data_er(mk_mid-1);
F2 = 1.5 * r_data_er(mk_mid + 1) - 0.5 * r_data_er(mk_mid + 2) - 0.5 * r_data_er(mk_mid) - 0.5 * r_data_er(mk_mid-1);
F3 = r_data_er(mk_mid);
uk_mid = (m(k) + m(k - 1) + u(k) + u(k - 1)) / 2 - mk_mid;
yI_mid = (F1 * uk_mid + F2) * uk_mid + F3;
time_error(ms) = yI_mid * (yI(k) - yI(k - 1));
else
time_error(ms) = 0;
end

if (ms > 1)
w(ms + 1) = w(ms) + c1 * (time_error(ms) - time_error(ms - 1)) + c2 * time_error(ms);
else
w(ms + 1) = w(ms) + c1 * time_error(ms) + c2 * time_error(ms);
end
ms = ms +1;
k = k + 1;
end

end

%抽样判决
%yI_decision = (sign(yI(1:2:end)) + 1) / 2;
yI_decision = (yI(1:1:end)) > 0;
%errbit(i) = sum(abs(yI_decision(301:N) - data(301:N))) / (N - 300);
errbit(i) = sum(abs(yI_decision - data)) / N;
i = i +1;
end

if is_plot
% figure;
% subplot(2, 1, 1);
% stem(data, 'filled');
% subplot(2, 1, 2);
% stem(yI_decision, 'filled');
figure;
ebn0 = 1:10;
ber_theory = 0.5 * erfc(sqrt(power(10, ebn0 / 10))); %理论误码率
semilogy(ebn0, ber_theory, '-or');
grid on; hold on;
semilogy(ebn0, errbit, '-sb');
grid on; hold on;
xlabel('EbN0(dB)'); ylabel('BER')
title('BPSK误码率曲线图');
legend('理论误码率曲线', '仿真误码率曲线');
% figure;
% plot(w); title('环路滤波器输出');
figure;
plot(n);

end

借鉴学习

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
clc;
clear;
close all;

% --------------------------参数定义-------------------------- %
dataLength = 2^19;%数据长度
M = 2;%调制阶数

rolloffFactor = 0.35;%滚降系数
span = 6;%截断符号范围
numSPSForming = 4;%成型过采样倍数
symbolRate = 1e6;%符号速率
freqSampling = numSPSForming*symbolRate;%采样速率

freqCarrier = 2e6;%载波频率

EbN0 = 1:10;%dB
SNR = EbN0+10*log10(2)-10*log10(numSPSForming);%信噪比(dB)

% BL = 2e3;%环路噪声带宽
% Ts = 1/symbolRate;%符号持续时间
% wn = 8*0.707*BL/(4*0.707^2+1);
% C = 4+4*0.707*wn*Ts+(wn*Ts)^2;
% k1 = (4*(wn*Ts)^2+8*0.707*wn*Ts)/C;%环路滤波器系数
% k2 = 4*(wn*Ts)^2/C;%环路滤波器系数

BL = 2e3; % 环路带宽
k1 = 8*BL*1/symbolRate/3;
k2 = 32*(BL*1/symbolRate)^2/9;

% --------------------------发射序列-------------------------- %
dataOrig = randi([0,M-1],1,dataLength);%0,1等概率

% --------------------------映射-------------------------- %
dataPsk = 2*dataOrig-1;
% pskData = real(pskmod(origData,order));

% --------------------------成型-------------------------- %
dataUpSample = upsample(dataPsk,numSPSForming);%上采样
filterForm = rcosdesign(rolloffFactor,span,numSPSForming,"sqrt");%成型滤波器
signalBaseBand = conv(dataUpSample,filterForm,"same");

% --------------------------上变频-------------------------- %
t = (0:(length(signalBaseBand)-1))/freqSampling;
carrier = cos(2*pi*freqCarrier*t);
signalModulated = signalBaseBand.*carrier;

actualBitErrorRate = zeros(1,10);%实际误码率
filterMatched = conj(filterForm);%匹配滤波器
for cntSNR = 1:length(SNR)
% --------------------------加噪-------------------------- %
signalNoizeAdded = awgn(signalModulated,SNR(cntSNR),'measured');
% signalNoizeAdded = awgn(signalModulated,10,'measured');

% --------------------------下变频-------------------------- %
signalDownFreq = signalNoizeAdded.*carrier;
dataResmp = resample(signalDownFreq,1e4,1e4+1);

% --------------------------匹配滤波-------------------------- %
dataMatched = conv(dataResmp,filterMatched,"same");

% --------------------------定时同步-------------------------- %
dataRecv = dataMatched;

outSyn = zeros(1,dataLength);%定时同步输出
timeError = zeros(1,dataLength);%定时误差
outFilter = [0.25 zeros(1,dataLength-1)];%环路滤波器输出,初值0.25
outNCO = [0.35 zeros(1,length(dataRecv)-1)];%NCO输出,初值0.9
outNCOTemp = outNCO;
outInterval = zeros(1,dataLength);%分数间隔,初值0.6

mkArray = zeros(1,dataLength);%保存各mk值,供计算中间点
indexSyn = 1;%同步输出索引,指示插值输出、分数间隔、环路滤波器

for indexSmp = 1:length(dataRecv)-2%采样点索引,指示输入数据、数控振荡器
% --------------------------定时控制-------------------------- %
outNCOTemp(indexSmp+1) = outNCO(indexSmp)-outFilter(indexSyn);
if outNCOTemp(indexSmp+1)>0
outNCO(indexSmp+1) = outNCOTemp(indexSmp+1);
else%NCO溢出,插值
outNCO(indexSmp+1) = mod(outNCOTemp(indexSmp+1),1);
outInterval(indexSyn) = outNCO(indexSmp)/outFilter(indexSyn);%分数间隔计算

% --------------------------插值-------------------------- %
mkArray(indexSyn) = indexSmp;
f0 = dataRecv(indexSmp);
f1 = -1/3*dataRecv(indexSmp-1)-1/2*dataRecv(indexSmp)+dataRecv(indexSmp+1)-1/6*dataRecv(indexSmp+2);
f2 = 1/2*dataRecv(indexSmp-1)-dataRecv(indexSmp)+1/2*dataRecv(indexSmp+1);
f3 = -1/6*dataRecv(indexSmp-1)+1/2*dataRecv(indexSmp)-1/2*dataRecv(indexSmp+1)+1/6*dataRecv(indexSmp+2);
outSyn(indexSyn) = f3*outInterval(indexSyn)^3+f2*outInterval(indexSyn)^2+f1*outInterval(indexSyn)+f0;

% --------------------------定时误差检测-------------------------- %
if indexSyn > 1
mkMid = floor( (mkArray(indexSyn-1)+outInterval(indexSyn-1)+mkArray(indexSyn)+outInterval(indexSyn))/2 );%中间点mk
intervalMid = (mkArray(indexSyn-1)+outInterval(indexSyn-1)+mkArray(indexSyn)+outInterval(indexSyn))/2-mkMid;%中间点分数间隔
f0Mid = dataRecv(mkMid);
f1Mid = -1/3*dataRecv(mkMid-1)-1/2*dataRecv(mkMid)+dataRecv(mkMid+1)-1/6*dataRecv(mkMid+2);
f2Mid = 1/2*dataRecv(mkMid-1)-dataRecv(mkMid)+1/2*dataRecv(mkMid+1);
f3Mid = -1/6*dataRecv(mkMid-1)+1/2*dataRecv(mkMid)-1/2*dataRecv(mkMid+1)+1/6*dataRecv(mkMid+2);
midOutSyn = f3Mid*intervalMid^3+f2Mid*intervalMid^2+f1Mid*intervalMid+f0Mid;%中间点插值
timeError(indexSyn+1) = midOutSyn*(outSyn(indexSyn)-outSyn(indexSyn-1));%定时误差更新
else
timeError(indexSyn+1) = 0;
end

% --------------------------环路滤波-------------------------- %
outFilter(indexSyn+1) = outFilter(indexSyn)+k1*(timeError(indexSyn+1)-timeError(indexSyn))+k2*timeError(indexSyn+1);
indexSyn = indexSyn+1;
end
end

% plot(outNCO);
% figure(2);
% plot(outFilter);
% figure(3);
% plot(outInterval);
% figure(4);
% plot(timeError);
% --------------------------判决-------------------------- %
dataDecision = outSyn;
dataDecision(dataDecision>0) = 1;
dataDecision(dataDecision<0) = 0;

% --------------------------误码率计算-------------------------- %
cmp = dataDecision-dataOrig;
error = cmp(cmp~=0);
actualBitErrorRate(cntSNR) = length(error)/dataLength;
end

% --------------------------误码率比较-------------------------- %
theoBitErrorRate = berawgn(EbN0,'psk',M,'nondiff');
figure;
semilogy(EbN0,theoBitErrorRate);
hold on;
semilogy(EbN0,actualBitErrorRate);
xlabel('Eb/N0 (dB)');
ylabel('BER');
legend('Theoretical','Actual');
title('Bit Error Rate');
grid on;