对 DSSS/BPSK 信号进行二次功率谱估计以提取载频信息

一、核心原理：为什么平方能载频？

1.1 DSSS/BPSK 信号模型

\[x(t) = A \cdot d(t) \cdot c(t) \cdot \cos(2\pi f_c t + \phi) + n(t) \]

其中：

• \(d(t)\)：信息码（±1），速率 \(R_b\)
• \(c(t)\)：扩频码（±1），码片速率 \(R_c = N \cdot R_b\)（\(N\) 为扩频增益）
• \(f_c\)：载波频率（待估计）
• \(n(t)\)：加性高斯白噪声

1.2 平方运算的数学推导

对接收信号进行平方：

\[\begin{aligned} x^2(t) &= A^2 \cdot d^2(t) \cdot c^2(t) \cdot \cos^2(2\pi f_c t + \phi) + \text{交叉项} \\ &= A^2 \cdot \cos^2(2\pi f_c t + \phi) + \text{噪声项} \quad (\text{因为 } d^2(t)=1, c^2(t)=1) \\ &= \frac{A^2}{2} + \frac{A^2}{2} \cos(4\pi f_c t + 2\phi) + \text{噪声项} \end{aligned} \]

关键结论：平方后，±1 的扩频调制和信息调制被消除，产生：

1. 直流分量 \(\frac{A^2}{2}\)
2. 二倍频分量 \(\frac{A^2}{2} \cos(4\pi f_c t + 2\phi)\)
3. 噪声与信号的交叉项

因此，对 \(x^2(t)\) 做功率谱分析，会在 \(2f_c\) 处出现谱线，从而可估计 \(f_c\)。

---

二、算法实现（MATLAB）

2.1 信号生成（用于验证）

%% DSSS/BPSK 信号生成
clear; clc; close all;

% 参数设置
fc = 10e6;          % 载频 10 MHz
fs = 50e6;          % 采样率 50 MHz
Rb = 100e3;         % 信息速率 100 kbps
Rc = 10e6;          % 码片速率 10 Mcps
N = Rc/Rb;          % 扩频增益 = 100
T = 1e-3;           % 信号时长 1 ms
SNR_dB = 10;        % 信噪比

% 生成时间轴
t = 0:1/fs:T-1/fs;
Ns = length(t);

% 生成信息码和扩频码
info_bits = randi([0 1], 1, floor(Rb*T));
pn_seq = randi([0 1], 1, floor(Rc*T));

% 扩频：信息码每个比特重复 N 次
info_expanded = repelem(info_bits, N);
% BPSK调制：0->1, 1->-1
info_bpsk = 2*info_expanded(1:length(pn_seq)) - 1;
pn_bpsk = 2*pn_seq - 1;

% DSSS信号（基带）
dsss_baseband = info_bpsk .* pn_bpsk;

% 上变频到载频（实信号）
dsss_passband = dsss_baseband .* cos(2*pi*fc*t(1:length(dsss_baseband)));

% 添加噪声
signal_power = mean(abs(dsss_passband).^2);
noise_power = signal_power / (10^(SNR_dB/10));
noise = sqrt(noise_power/2) * randn(size(dsss_passband));
received_signal = dsss_passband + noise;

% 截取到相同长度
if length(received_signal) < Ns
    received_signal = [received_signal, zeros(1, Ns-length(received_signal))];
else
    received_signal = received_signal(1:Ns);
end

2.2 平方倍频法估计载频

%% 方法1：直接平方法（平方律检测）
x_squared = received_signal .^ 2;

% Welch功率谱估计
nfft = 2^nextpow2(length(x_squared));
window = hamming(1024);
noverlap = 512;
[Pxx, f] = pwelch(x_squared, window, noverlap, nfft, fs);

% 寻找峰值（跳过直流附近）
skip_bins = round(0.01 * length(f));  % 跳过前1%的频率（避免直流影响）
[~, idx] = max(Pxx(skip_bins:end));
peak_idx = idx + skip_bins - 1;
f_peak = f(peak_idx);

% 载频估计
fc_est1 = f_peak / 2;

fprintf('方法1（直接平方法）:\n');
fprintf('  检测到的峰值频率: %.3f MHz\n', f_peak/1e6);
fprintf('  估计的载频: %.3f MHz\n', fc_est1/1e6);
fprintf('  真实载频: %.3f MHz\n', fc/1e6);
fprintf('  估计误差: %.3f kHz\n', abs(fc_est1 - fc)/1e3);

2.3 改进方法：延时相乘（避免直流影响）

%% 方法2：延时相乘法（更稳健）
tau = 1/(4*fc);  % 延时约1/4载波周期
delay_samples = round(tau * fs);
if delay_samples < 1
    delay_samples = 1;
end

% 延时相乘（相当于带通滤波的平方）
x_delayed = [zeros(1, delay_samples), received_signal(1:end-delay_samples)];
x_multiplied = received_signal .* x_delayed;

% 功率谱估计
[Pxx2, f2] = pwelch(x_multiplied, window, noverlap, nfft, fs);

% 寻找峰值
[~, idx2] = max(Pxx2);
f_peak2 = f2(idx2);

% 载频估计（注意：延时相乘产生的是f_c分量，不是2f_c）
fc_est2 = f_peak2;

fprintf('\n方法2（延时相乘法）:\n');
fprintf('  检测到的峰值频率: %.3f MHz\n', f_peak2/1e6);
fprintf('  估计的载频: %.3f MHz\n', fc_est2/1e6);
fprintf('  估计误差: %.3f kHz\n', abs(fc_est2 - fc)/1e3);

2.4 可视化结果

%% 可视化
figure('Position', [100, 100, 1200, 800]);

% 原始信号时域
subplot(3,2,1);
plot(t(1:2000)*1e6, received_signal(1:2000));
xlabel('时间 (\mu s)'); ylabel('幅度');
title('接收信号（时域）');
grid on;

% 原始信号功率谱
subplot(3,2,2);
[P_orig, f_orig] = pwelch(received_signal, window, noverlap, nfft, fs);
plot(f_orig/1e6, 10*log10(P_orig));
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
title('原始信号功率谱');
xlim([0, fs/2/1e6]); grid on;

% 平方后信号时域
subplot(3,2,3);
plot(t(1:2000)*1e6, x_squared(1:2000));
xlabel('时间 (\mu s)'); ylabel('幅度');
title('平方后信号（时域）');
grid on;

% 平方后信号功率谱（方法1）
subplot(3,2,4);
plot(f/1e6, 10*log10(Pxx));
hold on;
plot(f_peak/1e6, 10*log10(Pxx(peak_idx)), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
title(sprintf('平方后功率谱（峰值在 %.3f MHz）', f_peak/1e6));
xlim([0, fs/2/1e6]); grid on;
legend('功率谱', '检测峰值');

% 延时相乘信号时域
subplot(3,2,5);
plot(t(1:2000)*1e6, x_multiplied(1:2000));
xlabel('时间 (\mu s)'); ylabel('幅度');
title('延时相乘信号（时域）');
grid on;

% 延时相乘信号功率谱（方法2）
subplot(3,2,6);
plot(f2/1e6, 10*log10(Pxx2));
hold on;
plot(f_peak2/1e6, 10*log10(Pxx2(idx2)), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
title(sprintf('延时相乘功率谱（峰值在 %.3f MHz）', f_peak2/1e6));
xlim([0, fs/2/1e6]); grid on;
legend('功率谱', '检测峰值');

sgtitle('DSSS/BPSK信号载频估计 - 二次功率谱方法');

---

三、性能分析与改进

3.1 理论性能分析

平方倍频法的输出信噪比（SNR）与输入 SNR 的关系：

\[\text{SNR}_{out} \approx \frac{\text{SNR}_{in}^2}{1 + 2\text{SNR}_{in}} \]

当 \(\text{SNR}_{in} \ll 1\) 时，\(\text{SNR}_{out} \propto \text{SNR}_{in}^2\)，性能急剧下降。

检测门限：通常需要 \(\text{SNR}_{in} > -15 \text{ dB}\) 才能可靠检测。

3.2 影响估计精度的因素

因素	影响	缓解措施
低信噪比	谱峰被噪声淹没	增加积累时间、使用循环谱分析
频谱泄漏	谱峰展宽，定位不准	使用合适的窗函数（如 Kaiser、Blackman）
频率分辨率	估计精度受限于 \(\Delta f = f_s/N_{FFT}\)	使用 CZT（Chirp-Z变换）进行局部细化
多径效应	产生多个谱峰	使用多信号分类（MUSIC）等超分辨算法

3.3 改进方法：CZT 细化频谱

%% 方法3：CZT细化频谱（提高频率分辨率）
f_start = fc_est1 * 0.9;  % 细化区间起点
f_end = fc_est1 * 1.1;    % 细化区间终点
M = 4096;                 % 细化点数

% 计算CZT
w = exp(-1j*2*pi*(f_end-f_start)/(M*fs));
a = exp(1j*2*pi*f_start/fs);
czt_result = czt(x_squared, M, w, a);

% CZT对应的频率轴
f_czt = linspace(f_start, f_end, M);

% 寻找峰值
[~, idx_czt] = max(abs(czt_result).^2);
fc_est_czt = f_czt(idx_czt) / 2;

fprintf('\n方法3（CZT细化）:\n');
fprintf('  细化区间: %.3f ~ %.3f MHz\n', f_start/1e6, f_end/1e6);
fprintf('  估计的载频: %.3f MHz\n', fc_est_czt/1e6);
fprintf('  估计误差: %.3f Hz\n', abs(fc_est_czt - fc));

3.4 改进方法：循环谱分析（抗噪声性能更好）

%% 方法4：循环谱分析（需要通信工具箱）
if exist('commcyclefreq', 'file')
    % 计算循环谱
    [Sx, alpha, f] = commcyclefreq(received_signal, fs, 'Method', 'fft');
    
    % 寻找循环频率 alpha = 2fc 处的峰值
    alpha_idx = find(alpha > 2*fc*0.8 & alpha < 2*fc*1.2);
    [~, max_idx] = max(max(abs(Sx(:, alpha_idx)), [], 1));
    fc_est_cyclic = alpha(alpha_idx(max_idx)) / 2;
    
    fprintf('\n方法4（循环谱分析）:\n');
    fprintf('  估计的载频: %.3f MHz\n', fc_est_cyclic/1e6);
    fprintf('  估计误差: %.3f kHz\n', abs(fc_est_cyclic - fc)/1e3);
else
    fprintf('\n方法4（循环谱分析）: 需要通信工具箱\n');
end

参考代码 二次谱法 www.youwenfan.com/contentcnv/79091.html

四、工程应用注意事项

4.1 实际系统考虑

1. 采样率选择：需满足 \(f_s > 4f_c\)（避免二倍频分量混叠）
2. 抗混叠滤波：平方前需进行适当滤波，防止高频分量混叠
3. 动态范围：平方运算会扩大信号动态范围，需注意ADC量化效应
4. 实时处理：FPGA实现时可用CORDIC算法计算平方，流水线FFT进行谱估计

4.2 与其他方法的比较

方法	优点	缺点	适用场景
平方倍频法	简单、实时性好	低SNR性能差、有3dB损失	高SNR、快速粗估计
延时相乘法	避免直流影响、更稳健	需要精确延时控制	中等SNR、存在直流干扰
循环谱法	抗噪声能力强、可同时估计多个参数	计算复杂、需要先验信息	低SNR、非协作侦察
四阶累积量	完全抑制高斯噪声	计算量极大、实现复杂	极低SNR、特殊应用

4.3 性能评估指标

%% 性能评估：不同SNR下的估计误差
SNR_range = -20:5:20;  % SNR范围
num_trials = 100;      % 蒙特卡洛次数
errors = zeros(length(SNR_range), 1);

for i = 1:length(SNR_range)
    SNR = SNR_range(i);
    trial_errors = zeros(num_trials, 1);
    
    for trial = 1:num_trials
        % 生成带噪声信号（同前）
        % ...（信号生成代码）
        
        % 平方倍频法估计
        % ...（估计代码）
        
        trial_errors(trial) = abs(fc_est - fc);
    end
    
    errors(i) = mean(trial_errors);
end

figure;
semilogy(SNR_range, errors/1e3, 'b-o', 'LineWidth', 2);
xlabel('输入 SNR (dB)');
ylabel('平均估计误差 (kHz)');
title('平方倍频法载频估计性能 vs SNR');
grid on;

---

五、总结

核心结论：对 DSSS/BPSK 信号进行平方运算（或延时相乘）后做功率谱分析，确实可以估计载频信息，原理是平方消除了 ±1 调制，使被抑制的载波分量以二倍频形式显现。

关键要点：

1. 方法有效性：在 SNR > -10 dB 时通常有效，估计精度可达 \(\Delta f \approx f_s/(2N_{FFT})\)
2. 工程实现：需注意抗混叠、动态范围、实时性等问题
3. 改进方向：低 SNR 时可结合循环谱、高阶累积量等方法
4. 应用场景：非协作通信侦察、频谱监测、信号参数盲估计

推荐流程：

接收信号 → 带通滤波（防混叠） → 平方/延时相乘 → 
功率谱估计（Welch + 加窗） → 峰值检测 → 
CZT细化（可选） → 载频估计（f_peak/2）