使用MATLAB实现平方倍频法对DSSS/BPSK信号进行载频估计
DSSS/BPSK信号模型与平方倍频法原理
信号模型
DSSS/BPSK信号可以表示为：
s(t) = A·d(t)·c(t)·cos(2πf_c t + φ) + n(t)
其中：
d(t)：数据序列（±1）c(t)：扩频码序列（±1）f_c：载波频率（待估计）φ：载波相位n(t)：加性高斯白噪声
平方倍频法原理
通过信号平方运算消除调制信息，产生载频的二次谐波分量：
s²(t) = A²·d²(t)·c²(t)·cos²(2πf_c t + φ) + 噪声项
= A²/2 · [1 + cos(4πf_c t + 2φ)] + 噪声项
由于 d²(t) = 1，c²(t) = 1，平方后的信号包含：
- 直流分量
- 2f_c 频率分量 ← 这是我们用于载频估计的关键分量
MATLAB实现：平方倍频法载频估计
1. 生成DSSS/BPSK测试信号
function [signal, t, fs, fc_true] = generate_dsss_bpsk_signal()
% 生成DSSS/BPSK测试信号
% 信号参数
fs = 100e6; % 采样频率 100MHz
T = 1e-3; % 信号时长 1ms
fc_true = 10.7e6; % 真实载频 10.7MHz
% 扩频参数
chip_rate = 1e6; % 码片速率 1MHz
code_length = 1023; % 扩频码长度
% 生成时间序列
t = 0:1/fs:T-1/fs;
N = length(t);
% 生成数据序列
data_rate = chip_rate / code_length; % 数据速率
data_bits = randi([0 1], 1, ceil(T * data_rate));
data_wave = repelem(2*data_bits-1, ceil(fs/data_rate));
data_wave = data_wave(1:N);
% 生成扩频码（Gold码）
gold_code = generate_gold_code(code_length);
spreading_wave = repelem(2*gold_code-1, ceil(fs/chip_rate));
spreading_wave = spreading_wave(1:N);
% 生成DSSS/BPSK信号
carrier = cos(2*pi*fc_true*t + pi/4); % 载波
signal = data_wave .* spreading_wave .* carrier;
% 添加高斯白噪声
SNR_dB = 10; % 信噪比
signal = awgn(signal, SNR_dB, 'measured');
fprintf('生成DSSS/BPSK信号:\n');
fprintf(' 真实载频: %.3f MHz\n', fc_true/1e6);
fprintf(' 采样频率: %.1f MHz\n', fs/1e6);
fprintf(' 码片速率: %.1f MHz\n', chip_rate/1e6);
fprintf(' 信噪比: %d dB\n', SNR_dB);
end
function gold_seq = generate_gold_code(length)
% 生成Gold序列
% 简化的Gold码生成器
r1 = randi([0 1], 1, 10); % 寄存器1初始状态
r2 = randi([0 1], 1, 10); % 寄存器2初始状态
gold_seq = zeros(1, length);
for i = 1:length
% 线性反馈移位寄存器
feedback1 = mod(r1(3) + r1(10), 2);
feedback2 = mod(r2(2) + r2(3) + r2(6) + r2(8) + r2(9) + r2(10), 2);
r1 = [feedback1, r1(1:end-1)];
r2 = [feedback2, r2(1:end-1)];
gold_seq(i) = mod(r1(end) + r2(end), 2);
end
end
2. 平方倍频法载频估计
function fc_estimated = square_law_frequency_estimation(signal, fs)
% 平方倍频法载频估计
% 信号平方运算
signal_squared = signal .^ 2;
% 计算功率谱密度
N = length(signal_squared);
nfft = 2^nextpow2(N);
[Pxx, f] = pwelch(signal_squared, hamming(nfft/4), nfft/8, nfft, fs);
% 转换为dB尺度
Pxx_db = 10*log10(Pxx);
% 寻找峰值对应的频率
% 注意：平方后的峰值在2*fc处
[~, peak_idx] = max(Pxx_db);
f_peak = f(peak_idx);
% 载频估计值 = 峰值频率 / 2
fc_estimated = f_peak / 2;
fprintf('平方倍频法估计结果:\n');
fprintf(' 平方后峰值频率: %.3f MHz\n', f_peak/1e6);
fprintf(' 估计载频: %.3f MHz\n', fc_estimated/1e6);
end
3. 完整的载频估计系统
function dsss_carrier_estimation_demo()
% DSSS/BPSK信号载频估计演示
close all; clc;
% 生成测试信号
[signal, t, fs, fc_true] = generate_dsss_bpsk_signal();
% 方法1：平方倍频法
fc_square = square_law_frequency_estimation(signal, fs);
% 方法2：循环谱分析法（对比方法）
fc_cyclic = cyclic_spectrum_estimation(signal, fs);
% 方法3：基于自相关的载频估计
fc_corr = correlation_frequency_estimation(signal, fs);
% 显示比较结果
fprintf('\n=== 载频估计结果比较 ===\n');
fprintf('真实载频: %.6f MHz\n', fc_true/1e6);
fprintf('平方倍频法: %.6f MHz (误差: %.3f kHz)\n', ...
fc_square/1e6, abs(fc_square-fc_true)/1e3);
fprintf('循环谱分析法: %.6f MHz (误差: %.3f kHz)\n', ...
fc_cyclic/1e6, abs(fc_cyclic-fc_true)/1e3);
fprintf('自相关法: %.6f MHz (误差: %.3f kHz)\n', ...
fc_corr/1e6, abs(fc_corr-fc_true)/1e3);
% 绘制结果图形
plot_estimation_results(signal, fs, fc_true, fc_square, t);
end
4. 对比方法：循环谱分析
function fc_estimated = cyclic_spectrum_estimation(signal, fs)
% 循环谱分析法载频估计
N = length(signal);
nfft = 2^nextpow2(sqrt(N));
% 计算循环谱（简化的实现）
% 在实际应用中可以使用专门的循环谱分析工具箱
[S, alpha, f] = simplified_cyclic_spectrum(signal, fs, nfft);
% 在循环频率alpha=2fc处寻找峰值
alpha_idx = find(alpha > 0, 1); % 寻找正循环频率
if ~isempty(alpha_idx)
[~, f_peak_idx] = max(abs(S(alpha_idx, :)));
fc_estimated = f(f_peak_idx);
else
fc_estimated = 0;
end
end
function [S, alpha, f] = simplified_cyclic_spectrum(signal, fs, nfft)
% 简化的循环谱计算
N = length(signal);
window = hamming(nfft);
noverlap = nfft/2;
% 计算频谱相关密度
[S, f] = pwelch(signal, window, noverlap, nfft, fs);
% 简化的循环频率轴
alpha = (-fs/2:fs/nfft:fs/2-fs/nfft);
% 这里应该实现完整的循环谱计算
% 为简化演示，返回标准功率谱
S = repmat(S', length(alpha), 1);
end
5. 基于自相关的载频估计
function fc_estimated = correlation_frequency_estimation(signal, fs)
% 基于自相关的载频估计
% 计算信号的自相关函数
[R, lags] = xcorr(signal, 'coeff');
center_idx = find(lags == 0);
% 寻找自相关函数的第一个峰值（对应载波周期）
search_range = center_idx+1:center_idx+round(fs/(10e6)); % 假设载频>10MHz
[~, peak_idx] = max(abs(R(search_range)));
lag_peak = lags(search_range(peak_idx));
% 载频估计 = 采样频率 / 峰值滞后
if lag_peak ~= 0
fc_estimated = fs / abs(lag_peak);
else
fc_estimated = 0;
end
end
6. 结果可视化函数
function plot_estimation_results(signal, fs, fc_true, fc_estimated, t)
% 绘制载频估计结果
figure('Position', [100, 100, 1400, 1000]);
% 1. 原始信号时域波形（前100个采样点）
subplot(3, 3, 1);
plot(t(1:100)*1e6, real(signal(1:100)));
title('DSSS/BPSK信号时域波形 (前100个采样点)');
xlabel('时间 (\mus)'); ylabel('幅度');
grid on;
% 2. 原始信号功率谱
subplot(3, 3, 2);
[P_orig, f_orig] = pwelch(signal, hamming(1024), 512, 1024, fs);
plot(f_orig/1e6, 10*log10(P_orig));
hold on;
plot([fc_true/1e6, fc_true/1e6], ylim, 'r--', 'LineWidth', 2);
title('原始信号功率谱密度');
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
legend('信号PSD', '真实载频', 'Location', 'best');
grid on;
% 3. 平方后信号功率谱
subplot(3, 3, 3);
signal_squared = signal .^ 2;
[P_sq, f_sq] = pwelch(signal_squared, hamming(1024), 512, 1024, fs);
plot(f_sq/1e6, 10*log10(P_sq));
hold on;
plot([2*fc_true/1e6, 2*fc_true/1e6], ylim, 'r--', 'LineWidth', 2, ...
'DisplayName', '真实2倍载频');
plot([2*fc_estimated/1e6, 2*fc_estimated/1e6], ylim, 'g:', ...
'LineWidth', 2, 'DisplayName', '估计2倍载频');
title('平方后信号功率谱密度');
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
legend('平方信号PSD', '真实2f_c', '估计2f_c', 'Location', 'best');
grid on;
% 4. 平方倍频法详细分析
subplot(3, 3, 4);
N = length(signal_squared);
nfft = 2^nextpow2(N);
[Pxx_detailed, f_detailed] = pwelch(signal_squared, hamming(nfft/4), ...
nfft/8, nfft, fs);
Pxx_db = 10*log10(Pxx_detailed);
% 寻找峰值区域
peak_region = (f_detailed > 0.8*2*fc_true) & (f_detailed < 1.2*2*fc_true);
f_peak_region = f_detailed(peak_region);
Pxx_peak_region = Pxx_db(peak_region);
plot(f_peak_region/1e6, Pxx_peak_region, 'b-', 'LineWidth', 2);
hold on;
[max_val, max_idx] = max(Pxx_peak_region);
f_max = f_peak_region(max_idx);
plot(f_max/1e6, max_val, 'ro', 'MarkerSize', 8, 'LineWidth', 2);
plot([2*fc_true/1e6, 2*fc_true/1e6], ylim, 'r--', 'LineWidth', 2);
title('平方倍频法峰值检测');
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
legend('功率谱', '检测峰值', '真实2f_c', 'Location', 'best');
grid on;
% 5. 不同信噪比下的性能分析
subplot(3, 3, 5);
snr_range = -5:2:15;
estimation_errors = zeros(size(snr_range));
for i = 1:length(snr_range)
signal_noisy = awgn(signal, snr_range(i), 'measured');
fc_est = square_law_frequency_estimation(signal_noisy, fs);
estimation_errors(i) = abs(fc_est - fc_true) / 1e3; % kHz
end
plot(snr_range, estimation_errors, 'o-', 'LineWidth', 2);
xlabel('信噪比 (dB)'); ylabel('估计误差 (kHz)');
title('平方倍频法在不同信噪比下的性能');
grid on;
% 6. 信号星座图（解调后）
subplot(3, 3, 6);
% 使用估计的载频进行解调
t_local = 0:1/fs:(length(signal)-1)/fs;
local_carrier = cos(2*pi*fc_estimated*t_local);
demod_signal = signal .* local_carrier;
% 低通滤波
[b, a] = butter(6, 2e6/(fs/2)); % 2MHz低通
filtered_signal = filter(b, a, demod_signal);
plot(real(filtered_signal(1:1000)), imag(filtered_signal(1:1000)), '.');
title('解调后信号星座图');
xlabel('同相分量'); ylabel('正交分量');
axis equal; grid on;
% 7. 估计误差统计
subplot(3, 3, [7, 8, 9]);
% 多次蒙特卡洛仿真
num_trials = 100;
errors_khz = zeros(1, num_trials);
for trial = 1:num_trials
[test_signal, ~, ~, test_fc_true] = generate_dsss_bpsk_signal();
test_fc_est = square_law_frequency_estimation(test_signal, fs);
errors_khz(trial) = abs(test_fc_est - test_fc_true) / 1e3;
end
histogram(errors_khz, 20);
xlabel('估计误差 (kHz)'); ylabel('出现次数');
title('平方倍频法载频估计误差分布 (100次蒙特卡洛仿真)');
grid on;
mean_error = mean(errors_khz);
std_error = std(errors_khz);
text(0.7, 0.9, sprintf('平均误差: %.2f kHz\n标准差: %.2f kHz', ...
mean_error, std_error), 'Units', 'normalized', ...
'BackgroundColor', 'white');
sgtitle('DSSS/BPSK信号载频估计 - 平方倍频法分析', 'FontSize', 14);
end
参考代码 平方倍频法 www.youwenfan.com/contentcnl/79819.html
性能分析与改进建议
1. 平方倍频法的优势
- 简单易实现：只需平方运算和FFT
- 抗调制能力强：有效消除BPSK调制影响
- 计算效率高：适合实时处理
2. 局限性及改进方法
function fc_improved = improved_square_law_estimation(signal, fs)
% 改进的平方倍频法
% 1. 带通滤波预处理
fc_rough = initial_coarse_estimation(signal, fs); % 粗估计
[b, a] = butter(4, [0.8*fc_rough, 1.2*fc_rough]/(fs/2));
filtered_signal = filter(b, a, signal);
% 2. 多次平方增强效果
signal_squared = filtered_signal .^ 2;
% 3. 使用高分辨率频谱估计
N = length(signal_squared);
nfft = 4 * 2^nextpow2(N); % 增加FFT点数
% 4. 抛物线插值提高频率分辨率
[Pxx, f] = pwelch(signal_squared, hamming(nfft/4), nfft/8, nfft, fs);
[~, peak_idx] = max(Pxx);
% 抛物线插值
if peak_idx > 1 && peak_idx < length(Pxx)
alpha = 0.5 * (Pxx(peak_idx-1) - Pxx(peak_idx+1)) / ...
(Pxx(peak_idx-1) - 2*Pxx(peak_idx) + Pxx(peak_idx+1));
f_peak = f(peak_idx) + alpha * (f(2)-f(1));
else
f_peak = f(peak_idx);
end
fc_improved = f_peak / 2;
end
function fc_coarse = initial_coarse_estimation(signal, fs)
% 初始粗估计
[Pxx, f] = pwelch(signal, hamming(1024), 512, 1024, fs);
[~, idx] = max(Pxx);
fc_coarse = f(idx);
end
3. 应用建议
- 信噪比要求：平方倍频法在SNR > 0 dB时性能较好
- 信号长度：建议信号持续时间包含至少100个载波周期
- 频率分辨率：频率估计精度受FFT分辨率限制，可通过插值改善
- 多信号平均：对多个信号段分别估计后取平均，提高稳定性
浙公网安备 33010602011771号