载波频偏粗补偿

发表于|更新于|工程实践
|浏览量:

载波粗捕获算法研究

引言

地面信关站基带高速数传设备在运行过程中通常会有一些大的频偏，导致载波同步锁相环无法收敛。为解决这一技术难题，系统需预先实施频偏粗补偿处理。粗补偿一般辅助序列和非辅助序列两种，辅助序列通常依靠导频的自相关进行频偏估计，而非辅助序列则一般通过去调制的方法，其不依赖于导频序列，但其复杂度会更高。

为提升基带数传设备的场景适应性与系统兼容性，本设计采用非辅助序列实现频偏粗补偿。该方案通过信号高次方非线性变换剥离调制相位，结合FFT频谱分析完成频偏捕获。具体通过构建M次方倍频器消除QPSK等数字调制影响，使残留载波频偏转化为单频信号，最终经FFT峰值检测获取粗估计值。

符号速率600Msps，中频1.2G时，主要达成以下指标：

  • 载波捕获范围：≥ ± 2MHz；
  • 多普勒变化率适应范围：±75kHz/s
  • 调制方式主要包括：BPSK、QPSK、8PSK、16QAM、16APSK、32APSK

理论分析

由于时间关系，相关公式无法完整推导并附上，以下仅根据去调制原理对其进行阐述：
所谓去调制，即通过倍频的方式将调制信号变为单音，单音所在的频率即为频偏，将其补偿即可

  • 对于MPSK调制，其通过M次方倍频即可，消除相位调制，得到$M\Delta f$的单音，进一步可得到频偏；
  • 对于QAM调制，其通过4次方即可得到谱线，即$4\Delta f$的单音；
  • 对于APSK调制，其本质上是由多个半径不同的PSK组成，所以理论上去调制的方法是：倍频次数为根据各环相位点数的最小公倍数
    • 例如16APSK（4+12点环）：LCM(4,12)=12次方；32APSK（需根据具体环结构计算，如4+12+16点环：LCM=48）；
    • 特殊现象：32APSK在16次方、32次方是也有谱线，其原因在于此时去掉了外圈和内圈，且外全都幅值较大，高次方倍频对幅值大的更明显，所以谱线也比较明显
    • 此外，根据论文李东波《一种高效的 APSK 调制信号载波频偏估计算法》，发现其内环均为QPSK调制，根据幅度判决，找出内环数据，外环数据置零，置零处理后的信号数据由于只保留了内环上的数据，可将其看作 QPSK 信号（数据点间随机插入了零数据），然后使用4倍频即可。❗注意：0不能扣除，0的位置保留了原信号相位信息的连续性。

代码附录

以下给出不同调制方式载波粗捕获及细补偿的MATLAB程序。

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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
% @ DATA: 2025-03-06-09:57:39
% @ Author: Poster
% @ Description: 验证各种调制方式的大频偏粗捕获
% @ Notes:
clc; clear; close all;
%----------------------parameters--------------------------%
IsPlot = 1; % ? 是否绘图
fs = 5e9; %:采样频率
fc = 1.2e9; %:载波频率
Rs = 600e6; %:符号速率
N = 1e5; %:符号数
%==========================================================%
%% Transmitter Processing
%---------------Modulation Signal Generation--------------%
ModulationType = 6; %:调制方式
switch ModulationType
    case 1
        M = 2; %:调制阶数
        Symbols = randi([0 M - 1], 1, N); %:随机产生符号
        ModulationSignal = pskmod(Symbols, M, 0);
        fprintf('ModulationType is BPSK\n');
    case 2
        M = 4; %:调制阶数
        Symbols = randi([0 M - 1], 1, N); %:随机产生符号
        ModulationSignal = pskmod(Symbols, M, pi / 4);
        fprintf('ModulationType is QPSK\n');
    case 3
        M = 8; %:调制阶数
        Symbols = randi([0 M - 1], 1, N); %:随机产生符号
        ModulationSignal = pskmod(Symbols, M, pi / 8);
        fprintf('ModulationType is 8PSK\n');
    case 4
        M = 16; %:调制阶数
        Symbols = randi([0 M - 1], 1, N); %:随机产生符号
        ModulationSignal = qammod(Symbols, M, 'gray');
        fprintf('ModulationType is 16QAM\n');
    case 5
        M = 16; %:调制阶数
        Symbols = randi([0 M - 1], 1, N); %:随机产生符号
        R1 = 1; R2 = 2;
        ModulationSignal = apskmod(Symbols, [4 12], [R1 R2]);
        fprintf('ModulationType is 16APSK\n');
    case 6
        M = 32; %:调制阶数
        Symbols = randi([0 M - 1], 1, N); %:随机产生符号
        R1 = 1; R2 = 2; R3 = 3;
        ModulationSignal = apskmod(Symbols, [4 12 16], [R1 R2 R3]);
        fprintf('ModulationType is 64QAM\n');
    otherwise
        error('Invalid modulation type.');
end
%==========================================================%
%---------------------RRC Shaping Filter------------------%
alpha = 0.35; sps = 4;
rcosFIR = rcosdesign(alpha, 8, sps, 'sqrt');
ModulationSignalUp = upsample(ModulationSignal, sps);
SignalShaping = conv(ModulationSignalUp, rcosFIR);
%==========================================================%
%----------------------DUC Processing------------------%
%:正交上变频
t = (0:length(SignalShaping) - 1) / fs;
SignalUpconverter = SignalShaping .* exp(1i * 2 * pi * fc * t);
SignalSend = real(SignalUpconverter);
%% Receiver Processing
%==========================================================%
%----------------------DDC Processing------------------%
Offset = 2e6; %->频偏
RecvSignal_I = SignalSend .* cos(2 * pi * (fc - Offset) * t);
RecvSignal_Q = SignalSend .* sin(2 * pi * (fc - Offset) * t);
RecvSignal = RecvSignal_I - 1i * RecvSignal_Q;
% RecvSignal = SignalUpconverter .* exp(-1i * 2 * pi * (fc-Offset) * t);
%==========================================================%
%--------------------------匹配滤波--------------------------%
RecvSignalMatched = conv(RecvSignal, rcosFIR);
%==========================================================%
%--------------------------降采样--------------------------%
RecvSignalDown = RecvSignalMatched(length(rcosFIR):sps:end - length(rcosFIR) - 1);
%==========================================================%
%--------------------------去调制--------------------------%
switch ModulationType
    case 1
        order = 2;
        RecvSignalRomove = RecvSignalDown .^ order;
    case 2
        order = 4;
        RecvSignalRomove = RecvSignalDown .^ order;
    case 3
        order = 8;
        RecvSignalRomove = RecvSignalDown .^ order;
    case 4
        order = 4;
        RecvSignalRomove = RecvSignalDown .^ order;
    case 5
        order = 12;
        RecvSignalRomove = RecvSignalDown .^ order;
    case 6
        order = 16;
        RecvSignalRomove = RecvSignalDown .^ order;
    otherwise
        error('Invalid modulation type.');
end
%--------------基于FFT的信号信号载波频偏估计-------------------%
% NFFT = 16384; %->FFT点数
% FreResolution = fs / NFFT; %:频率分辨率
% fprintf('FreResolution is %f kHz\n', FreResolution / 1e3);
FreResolution = 100e3; %:频率分辨率
fs_actual = fs / sps; %:抽完了, 采样率对应降低
NFFT = pow2(nextpow2((fs_actual / FreResolution))); % 确保FFT点数为2的幂次方
fprintf('Actual FreResolution is %f kHz\n', fs_actual / NFFT / 1e3);
RecvSignalRomoveFFT = fft(RecvSignalRomove, NFFT);
%* 计算频偏
[~, FreIndex] = max(abs(RecvSignalRomoveFFT));
if FreIndex > NFFT / 2
    FreIndex = FreIndex - NFFT; %:处理负频偏
end
FreOffset = (FreIndex - 1) * fs_actual / NFFT;
switch ModulationType
    case 1
        FreOffsetValid = FreOffset / order;
        fprintf('The estimation of frequency offset is %f kHz\n', FreOffsetValid / 1e3);
    case 2
        FreOffsetValid = FreOffset / order;
        fprintf('The estimation of frequency offset is %f kHz\n', FreOffsetValid / 1e3);
    case 3
        FreOffsetValid = FreOffset / order;
        fprintf('The estimation of frequency offset is %f kHz\n', FreOffsetValid / 1e3);
    case 4
        FreOffsetValid = FreOffset / order;
        fprintf('The estimation of frequency offset is %f kHz\n', FreOffsetValid / 1e3);
    case 5
        FreOffsetValid = FreOffset / order;
        fprintf('The estimation of frequency offset is %f kHz\n', FreOffsetValid / 1e3);
    case 6
        FreOffsetValid = FreOffset / order;
        fprintf('The estimation of frequency offset is %f kHz\n', FreOffsetValid / 1e3);
    otherwise
        error('Invalid modulation type.');
end
FreOffsetError = FreOffsetValid - Offset;
fprintf('The estimation of frequency offset error is %f kHz\n', FreOffsetError / 1e3);
%==========================================================%
%--------------------------频偏粗补偿--------------------------%
t_new = (0:length(RecvSignalDown) - 1) / fs_actual;
RecvSignalDownComp = RecvSignalDown .* exp(-1i * 2 * pi * FreOffsetValid * t_new);
%==========================================================%
%--------------------------锁相环--------------------------%
RecvCarrierSync = pllSync(RecvSignalDownComp, ModulationType, fs_actual);

%% Figure Plot
if IsPlot
    scatterplot(ModulationSignal); title('发端信号星座图');
    % scatterplot(RecvSignalDown); title('接收端经匹配滤波下采样后的信号星座图');
    figure(3);
    % plot(abs(RecvSignalRomoveFFT));
    plot(fs_actual * (-NFFT / 2:NFFT / 2 - 1) / NFFT, abs(fftshift(RecvSignalRomoveFFT)));
    title('频偏估计'); xlabel('频率(Hz)'); ylabel('幅度'); grid on;
    scatterplot(RecvSignalDownComp); title('接收端经频偏粗补偿后的信号星座图');
    scatterplot(RecvCarrierSync(3e4:end)); title('接收端经频偏细补偿后的信号星座图');

end
%% Function Definitions
function PLL_OUT = pllSync(RecvSignalDownComp, ModulationType, fs)
    %-> 锁相环频偏细补偿
    % RecvSignalDown: 接收端经匹配滤波下采样且经过粗补偿
    % ModulationType: 调制方式
    % fs: 采样率
    % PllOut: 锁相环输出
    switch ModulationType
        case 1
            M = 2;
            PLL_OUT = MPSKPLL(M, RecvSignalDownComp, fs) * exp(-1i * pi / 2); % MPSK鉴相初相为π/2
        case 2
            M = 4;
            PLL_OUT = MPSKPLL(M, RecvSignalDownComp, fs);
        case 3
            M = 8;
            PLL_OUT = MPSKPLL(M, RecvSignalDownComp, fs);
        case 4
            PLL_OUT = QAMPLL(RecvSignalDownComp, fs);
        case 5
            M = [4 12];
            PLL_OUT = MAPSKPLL(M, RecvSignalDownComp);
        case 6
            M = [4 12 16];
            PLL_OUT = MAPSKPLL(M, RecvSignalDownComp);
    end
    function MPSK_Output = MPSKPLL(M, RecvSignalDownComp, fs)
        N_Valid = length(RecvSignalDownComp);
        theta = zeros(1, N_Valid);
        PLL_Detector = zeros(1, N_Valid); % 鉴相器输出
        PLL_NCO = zeros(1, N_Valid); % NCO输出
        PLL_LOOP_OUT = zeros(1, N_Valid); % 环路滤波器输出
        A = zeros(1, N_Valid);
        B = zeros(1, N_Valid);
        I_PLL = zeros(1, N_Valid);
        Q_PLL = zeros(1, N_Valid);
        %: 计算环路滤波器的参数
        BL = 27e6;
        Wn = BL / 0.53; % 固有振荡频率
        T = 1 / fs;
        PLL_C1 = (4 * (Wn * T) ^ 2 + 4 * sqrt(2) * Wn * T) / (4 + 2 * sqrt(2) * Wn * T + (Wn * T) ^ 2);
        PLL_C2 = (4 * (Wn * T) ^ 2) / (4 + 2 * sqrt(2) * Wn * T + (Wn * T) ^ 2);
        % 初始化第一个输出值
        PLL_OUT = [RecvSignalDownComp(1) zeros(1, N_Valid - 1)];
        % 遍历每个样本点进行处理
        for i = 2:N_Valid
            PLL_OUT(i) = RecvSignalDownComp(i) * exp(-1i * PLL_NCO(i));
            %* 数字鉴相器
            I_PLL(i) = real(PLL_OUT(i));
            Q_PLL(i) = imag(PLL_OUT(i));
            theta(i) = atan2(Q_PLL(i), I_PLL(i));
            phi = mod(M * theta(i) - pi, 2 * pi);
            if (phi <= pi)
                PLL_Detector(i) = phi / M;
            else
                PLL_Detector(i) = (phi - 2 * pi) / M;
            end
            %* 环路滤波器
            A(i) = PLL_Detector(i) * PLL_C2 + A(i - 1);
            B(i) = PLL_C1 * PLL_Detector(i);
            PLL_LOOP_OUT(i) = A(i - 1) + B(i);
            %* 更新NCO输出
            PLL_NCO(i + 1) = PLL_NCO(i) + PLL_LOOP_OUT(i);
            PLL_NCO(i + 1) = mod(PLL_NCO(i + 1), 2 * pi);
        end
        MPSK_Output = PLL_OUT;
    end
    function QAM_output = QAMPLL(RecvSignalDownComp, ~)
        N_Valid = length(RecvSignalDownComp);
        theta = zeros(1, N_Valid);
        PLL_Detector = zeros(1, N_Valid); %鉴相器输出
        PLL_NCO = zeros(1, N_Valid); %NCO输出
        PLL_LOOP_OUT = zeros(1, N_Valid); %环路滤波器输出
        R = zeros(1, N_Valid); %幅度平方
        A = zeros(1, N_Valid);
        B = zeros(1, N_Valid);
        I_PLL = zeros(1, N_Valid);
        Q_PLL = zeros(1, N_Valid);
        Standard_Amplitude = sqrt((2 * 4 + 10 * 8 + 18 * 4) / 16); %标准幅度
        Timing_OUT_Nomal = Standard_Amplitude * (RecvSignalDownComp / (sqrt(mean(abs(RecvSignalDownComp) .^ 2)))); %归一化到标准星座点
        PLL_OUT = [Timing_OUT_Nomal(1) zeros(1, N_Valid - 1)];
        %准确求法
        % BL = 15e6; %环路噪声带宽
        % Wn = BL / 0.53; %固有振荡频率
        % T = 1 / fs;
        % PLL_C1 = (4 * (Wn * T) ^ 2 + 4 * sqrt(2) * Wn * T) / (4 + 2 * sqrt(2) * Wn * T + (Wn * T) ^ 2);
        % PLL_C2 = (4 * (Wn * T) ^ 2) / (4 + 2 * sqrt(2) * Wn * T + (Wn * T) ^ 2);
        PLL_C1 = 1/2 ^ 12; PLL_C2 = PLL_C1 ^ 2/2;

        for i = 2:N_Valid
            PLL_OUT(i) = Timing_OUT_Nomal(i) * exp(-1i * PLL_NCO(i));
            R(i) = imag(PLL_OUT(i)) .^ 2 + real(PLL_OUT(i)) .^ 2;
            % * 数字鉴相器
            I_PLL(i) = real(PLL_OUT(i));
            Q_PLL(i) = imag(PLL_OUT(i));
            theta(i) = atan2(Q_PLL(i), I_PLL(i));
            phi = mod(4 * theta(i) - pi, 2 * pi);
            if (phi <= pi)
                theta(i) = phi / 4;
            else
                theta(i) = (phi - 2 * pi) / 4;
            end
            if (R(i) <= (3 + sqrt(5)) || R(i) >= (7 + 3 * sqrt(5)))
                PLL_Detector(i) = theta(i);
            else
                PLL_Detector(i) = theta(i - 1);
            end
            % * 环路滤波器
            A(i) = PLL_Detector(i) * PLL_C2 + A(i - 1); %频率响应函数
            % A(i) = (PLL_Detector(i) - PLL_Detector(i-1)) * PLL_C2;
            B(i) = PLL_C1 * PLL_Detector(i); %相位响应函数
            PLL_LOOP_OUT(i) = A(i - 1) + B(i);
            % * NCO输出
            PLL_NCO(i + 1) = PLL_NCO(i) + PLL_LOOP_OUT(i);
            PLL_NCO(i + 1) = mod(PLL_NCO(i + 1), 2 * pi);
        end
        QAM_output = PLL_OUT;
    end
    function MAPSK_Output = MAPSKPLL(M, RecvSignalDownComp)
        PLL_C1 = 1/2 ^ 4; PLL_C2 = PLL_C1 ^ 2/2;
        Timing_OUT_Nomal = (RecvSignalDownComp / (sqrt(mean(abs(RecvSignalDownComp) .^ 2)))); %归一化
        N_Valid = length(Timing_OUT_Nomal);
        PLL_Detector = zeros(1, N_Valid); %鉴相器输出
        PLL_NCO = zeros(1, N_Valid); %NCO输出
        PLL_LOOP_OUT = zeros(1, N_Valid); %环路滤波器输出
        A = zeros(1, N_Valid);
        B = zeros(1, N_Valid);
        PLL_OUT = [Timing_OUT_Nomal(1) zeros(1, N_Valid - 1)];
        for i = 2:N_Valid
            PLL_OUT(i) = Timing_OUT_Nomal(i) * exp(-1i * PLL_NCO(i));
            % *数字鉴相器
            if (sum(M) == 16)
                PLL_OUT_Pro = PLL_OUT(i) ^ 3;
                PLL_Detector_Pro = PLL_OUT_Pro * (sign(real(PLL_OUT_Pro)) - 1i * sign(imag(PLL_OUT_Pro)));
                PLL_Detector(i) = imag(PLL_Detector_Pro);
            else
                PLL_OUT_Pro_Before = PLL_OUT(i) ^ 4;
                PLL_OUT_Pro = PLL_OUT_Pro_Before; %!外圈初相不为0
                % PLL_OUT_Pro = PLL_OUT_Pro_Before * exp(1i * pi / 4); %!初相为0, 相移pi/4
                PLL_Detector_Pro = PLL_OUT_Pro * (sign(real(PLL_OUT_Pro)) - 1i * sign(imag(PLL_OUT_Pro)));
                PLL_Detector(i) = imag(PLL_Detector_Pro);
            end
            %环路滤波器
            A(i) = PLL_Detector(i) * PLL_C2 + A(i - 1);
            B(i) = PLL_C1 * PLL_Detector(i);
            PLL_LOOP_OUT(i) = A(i) + B(i);
            %NCO输出
            PLL_NCO(i + 1) = PLL_NCO(i) + PLL_LOOP_OUT(i);
            PLL_NCO(i + 1) = mod(PLL_NCO(i + 1), 2 * pi);
        end
        MAPSK_Output = PLL_OUT;
    end
    function MAPSK_Output_2 = MAPSKPLL_2(M, RecvSignalDownComp)
        if (sum(M) == 16)
            R1 = 1; R2 = 2; %#ok<*NASGU>
            Standard_Amplitude = sqrt((4 * R1 ^ 2 + 12 * R2 ^ 2) / 16); %!标准幅度,需要根据模式进行修改
        else
            R1 = 1; R2 = 2; R3 = 3;
            Standard_Amplitude = sqrt((4 * R1 ^ 2 + 12 * R2 ^ 2 + 16 * R3 ^ 2) / 32); %!标准幅度,需要根据模式进行修改
        end
        % * 归一化幅值
        Average_Recv = sqrt(mean(abs(RecvSignalDownComp) .^ 2));
        Recv_Signal_Normalization = Standard_Amplitude * RecvSignalDownComp / Average_Recv;
        C1 = 1/2 ^ 12; C2 = C1 ^ 2/2;
        theta = zeros(1, N);
        PLL_Detector = zeros(1, N); %鉴相器输出
        NCO = zeros(1, N); %NCO输出
        LOOP_OUT = zeros(1, N); %环路滤波器输出
        R = zeros(1, N); %幅度平方
        A = zeros(1, N);
        B = zeros(1, N);
        I_PLL = zeros(1, N);
        Q_PLL = zeros(1, N);
        PLL_OUT = [Recv_Signal_Normalization(1) zeros(1, N - 1)];

        for i = 2:N
            PLL_OUT(i) = Recv_Signal_Normalization(i) * exp(-1i * NCO(i));
            % R(i) = abs(PLL_OUT(i)) .^ 2;
            R(i) = imag(PLL_OUT(i)) .^ 2 + real(PLL_OUT(i)) .^ 2;
            %数字鉴相器
            I_PLL(i) = real(PLL_OUT(i));
            Q_PLL(i) = imag(PLL_OUT(i));
            theta(i) = atan2(Q_PLL(i), I_PLL(i)); % atan2 [-pi,pi]  atan [-pi/2,pi/2]
            % theta(i) = (mod(4 * theta(i), 2 * pi) - pi)/4;% ?

            phi = mod(4 * theta(i) - pi, 2 * pi);
            if (phi <= pi)
                theta(i) = phi / 4;
            else
                theta(i) = (phi - 2 * pi) / 4;
            end
            % * 鉴相器(目前16APSK和32APSK都用的是内环)
            % * 鉴相器(16APSK也可以用内外环加快收敛速度)
            if (R(i) <= (9/4))
                PLL_Detector(i) = theta(i);
            else
                PLL_Detector(i) = theta(i - 1);
            end

            % * 环路滤波器
            A(i) = PLL_Detector(i) * C2 + A(i - 1);
            B(i) = C1 * PLL_Detector(i);
            LOOP_OUT(i) = A(i) + B(i);
            % * NCO输出
            NCO(i + 1) = NCO(i) + LOOP_OUT(i);
            NCO(i + 1) = mod(NCO(i + 1), 2 * pi);

        end
    end
end

avatar
Poster
路漫漫其修远兮，吾将上下而求索
文章
30
标签
15
分类
3
Follow Me
公告
This is my Blog
系列文章
载波频偏粗补偿
载波频偏粗补偿
FPGA-DDC
FPGA-DDC
CMA盲均衡原理及应用
CMA盲均衡原理及应用
滤波器设计实例
滤波器设计实例
RFSOC47DR AD-DA配置
RFSOC47DR AD-DA配置
generate用法
generate用法
时钟分频设计
时钟分频设计
通信发端设计
通信发端设计
通信发端设计
通信发端设计
成形滤波器设计
成形滤波器设计
最新文章
【温故而知新】滤波器结构设计
【温故而知新】滤波器结构设计
拓展学习记录
拓展学习记录
通信FPGA相关合集
通信FPGA相关合集
DFT-S-OFDM 波形设计
DFT-S-OFDM 波形设计
技能总结
技能总结