北斗/GPS接收机系统与信号仿真
北斗和GPS接收机是现代导航定位系统的核心设备，能够接收、处理卫星信号并计算用户位置。信号仿真是接收机设计、算法验证和性能评估的关键环节。
一、北斗/GPS接收机系统架构
1. 接收机基本组成
┌─────────────────────────────────────────────────────┐
│ 北斗/GPS接收机系统架构 │
├─────────────────────────────────────────────────────┤
│ 1. 射频前端 (RF Front-end) │
│ ├── 天线 (Antenna) │
│ ├── 低噪声放大器 (LNA) │
│ ├── 下变频器 (Down-converter) │
│ └── 模数转换器 (ADC) │
│ │
│ 2. 数字信号处理 (Digital Signal Processing) │
│ ├── 信号捕获 (Acquisition) │
│ ├── 信号跟踪 (Tracking) │
│ ├── 位同步与帧同步 (Bit/Frame Sync) │
│ └── 导航电文解调 (Navigation Data Demodulation) │
│ │
│ 3. 定位解算 (Positioning Solution) │
│ ├── 伪距测量 (Pseudorange Measurement) │
│ ├── 载波相位测量 (Carrier Phase Measurement) │
│ ├── 多普勒测量 (Doppler Measurement) │
│ └── 定位算法 (Positioning Algorithm) │
│ │
│ 4. 辅助系统 (Auxiliary Systems) │
│ ├── 惯性导航单元 (IMU) │
│ ├── 时钟系统 (Clock System) │
│ └── 数据接口 (Data Interface) │
└─────────────────────────────────────────────────────┘
2. 北斗与GPS信号特性对比
| 参数 | GPS L1 C/A | 北斗 B1I | 北斗 B1C | 北斗 B2a |
|---|---|---|---|---|
| 载波频率 | 1575.42 MHz | 1561.098 MHz | 1575.42 MHz | 1176.45 MHz |
| 码类型 | C/A码 | B1I码 | B1C数据/导频 | B2a数据/导频 |
| 码速率 | 1.023 MHz | 2.046 MHz | 1.023 MHz | 10.23 MHz |
| 码长 | 1023 chips | 2046 chips | 10230 chips | 10230 chips |
| 调制方式 | BPSK | BPSK | BOC(1,1) | AltBOC(15,10) |
| 数据速率 | 50 bps | 50 bps | 50/100 bps | 50/100 bps |
| 频宽 | 2.046 MHz | 4.092 MHz | 32.736 MHz | 32.736 MHz |
二、北斗/GPS信号仿真方法
1. 信号生成数学模型
北斗/GPS信号可以表示为：
s(t) = A·D(t)·C(t)·cos(2πf_c t + φ(t)) + n(t)
其中：
A：信号幅度D(t)：导航数据（±1）C(t)：伪随机码（±1）f_c：载波频率φ(t)：载波相位（包含多普勒频移）n(t)：加性高斯白噪声
2. MATLAB信号仿真实现
%% 北斗B1I信号仿真
classdef BeiDouB1ISignalGenerator
properties
% 信号参数
fc = 1561.098e6; % 载波频率 (Hz)
fs = 16.368e6; % 采样频率 (Hz)
code_rate = 2.046e6; % 码速率 (chips/s)
code_length = 2046; % 码长
data_rate = 50; % 数据速率 (bps)
prn = 1; % 卫星PRN号
cn0 = 45; % 载噪比 (dB-Hz)
% 仿真参数
duration = 0.001; % 信号时长 (s)
doppler = 5000; % 多普勒频移 (Hz)
code_phase = 100; % 码相位偏移 (chips)
carrier_phase = 0; % 载波相位 (rad)
end
methods
function [signal, t] = generate_signal(obj)
% 生成北斗B1I信号
% 1. 生成时间序列
num_samples = round(obj.duration * obj.fs);
t = (0:num_samples-1)' / obj.fs;
% 2. 生成导航数据
data_bits = randi([0 1], ceil(obj.duration * obj.data_rate), 1);
data_bits(data_bits == 0) = -1; % BPSK调制：0->-1, 1->1
data_signal = repelem(data_bits, round(obj.fs/obj.data_rate));
data_signal = data_signal(1:num_samples);
% 3. 生成伪随机码（简化版Gold码）
code = obj.generate_gold_code(obj.prn);
code_signal = repmat(code, ceil(num_samples/(obj.fs/obj.code_rate)), 1);
code_signal = code_signal(1:num_samples);
% 4. 生成载波
carrier_freq = obj.fc + obj.doppler;
carrier = cos(2*pi*carrier_freq*t + obj.carrier_phase);
% 5. 组合信号
signal = data_signal .* code_signal .* carrier;
% 6. 添加码相位偏移
code_phase_samples = round(obj.code_phase * obj.fs / obj.code_rate);
signal = circshift(signal, code_phase_samples);
% 7. 添加噪声
signal_power = mean(abs(signal).^2);
noise_power = signal_power / (10^(obj.cn0/10) * obj.duration);
noise = sqrt(noise_power/2) * (randn(size(signal)) + 1i*randn(size(signal)));
signal = signal + real(noise);
% 8. 归一化
signal = signal / max(abs(signal));
end
function code = generate_gold_code(obj, prn)
% 生成Gold码（简化实现）
% 实际北斗B1I使用截短的Gold码
% 生成两个m序列
g1 = obj.generate_m_sequence([10, 3], 1); % 多项式: x^10 + x^3 + 1
g2 = obj.generate_m_sequence([10, 9, 8, 6, 3, 2], 1); % 多项式: x^10 + x^9 + x^8 + x^6 + x^3 + x^2 + 1
% 生成Gold码
code = mod(g1 + circshift(g2, prn-1), 2);
code(code == 0) = -1; % 转换为±1
code = code(1:obj.code_length);
end
function seq = generate_m_sequence(obj, taps, init_state)
% 生成m序列
% taps: 反馈抽头位置
% init_state: 初始状态
n = max(taps); % 寄存器长度
register = de2bi(init_state, n, 'left-msb'); % 初始状态
seq_length = 2^n - 1;
seq = zeros(1, seq_length);
for i = 1:seq_length
% 输出位
seq(i) = register(end);
% 计算反馈位
feedback = mod(sum(register(taps)), 2);
% 移位
register = [feedback, register(1:end-1)];
end
end
function plot_signal_characteristics(obj, signal, t)
% 绘制信号特性
figure('Position', [100 100 1200 800]);
% 时域波形
subplot(3, 2, 1);
plot(t(1:1000), real(signal(1:1000)));
xlabel('时间 (s)'); ylabel('幅度');
title('信号时域波形（前1000点）');
grid on;
% 功率谱密度
subplot(3, 2, 2);
[pxx, f] = pwelch(signal, 1024, 512, 1024, obj.fs, 'centered');
plot(f/1e6, 10*log10(pxx));
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
title('信号功率谱密度');
grid on;
xlim([obj.fc/1e6-2, obj.fc/1e6+2]);
% 自相关函数
subplot(3, 2, 3);
[corr, lags] = xcorr(signal(1:10000), 'normalized');
plot(lags/obj.fs*obj.code_rate, corr);
xlabel('码片偏移'); ylabel('自相关系数');
title('信号自相关函数');
grid on;
xlim([-10, 10]);
% 星座图
subplot(3, 2, 4);
scatter(real(signal(1:1000)), imag(signal(1:1000)), 10, 'filled');
xlabel('同相分量'); ylabel('正交分量');
title('信号星座图');
axis equal; grid on;
% 直方图
subplot(3, 2, 5);
histogram(real(signal), 50, 'Normalization', 'pdf');
xlabel('幅度'); ylabel('概率密度');
title('信号幅度分布');
grid on;
% 频谱图
subplot(3, 2, 6);
spectrogram(signal(1:min(100000, length(signal))), 256, 250, 256, obj.fs, 'yaxis');
title('信号频谱图');
colorbar;
end
end
end
3. GPS L1 C/A信号仿真
%% GPS L1 C/A信号仿真
classdef GPSL1CASignalGenerator
properties
fc = 1575.42e6; % 载波频率 (Hz)
fs = 16.368e6; % 采样频率 (Hz)
code_rate = 1.023e6; % 码速率 (chips/s)
code_length = 1023; % 码长
data_rate = 50; % 数据速率 (bps)
prn = 1; % 卫星PRN号
cn0 = 45; % 载噪比 (dB-Hz)
% 仿真参数
duration = 0.001; % 信号时长 (s)
doppler = 5000; % 多普勒频移 (Hz)
code_phase = 100; % 码相位偏移 (chips)
carrier_phase = 0; % 载波相位 (rad)
end
methods
function [signal, t] = generate_signal(obj)
% 生成GPS L1 C/A信号
% 1. 生成时间序列
num_samples = round(obj.duration * obj.fs);
t = (0:num_samples-1)' / obj.fs;
% 2. 生成导航数据
data_bits = randi([0 1], ceil(obj.duration * obj.data_rate), 1);
data_bits(data_bits == 0) = -1;
data_signal = repelem(data_bits, round(obj.fs/obj.data_rate));
data_signal = data_signal(1:num_samples);
% 3. 生成C/A码
ca_code = obj.generate_ca_code(obj.prn);
code_signal = repmat(ca_code, ceil(num_samples/(obj.fs/obj.code_rate)), 1);
code_signal = code_signal(1:num_samples);
% 4. 生成载波
carrier_freq = obj.fc + obj.doppler;
carrier = cos(2*pi*carrier_freq*t + obj.carrier_phase);
% 5. 组合信号
signal = data_signal .* code_signal .* carrier;
% 6. 添加码相位偏移
code_phase_samples = round(obj.code_phase * obj.fs / obj.code_rate);
signal = circshift(signal, code_phase_samples);
% 7. 添加噪声
signal_power = mean(abs(signal).^2);
noise_power = signal_power / (10^(obj.cn0/10) * obj.duration);
noise = sqrt(noise_power/2) * (randn(size(signal)) + 1i*randn(size(signal)));
signal = signal + real(noise);
% 8. 归一化
signal = signal / max(abs(signal));
end
function ca_code = generate_ca_code(obj, prn)
% 生成GPS C/A码
% 基于Gold码，PRN 1-32对应不同的延迟组合
% G1和G2寄存器初始化
g1 = ones(1, 10); % 全1初始化
g2 = ones(1, 10);
ca_code = zeros(1, obj.code_length);
for i = 1:obj.code_length
% 输出C/A码
ca_code(i) = mod(g1(10) + g2(obj.get_g2_tap(prn)), 2);
% 更新G1寄存器
g1_feedback = mod(g1(3) + g1(10), 2);
g1 = [g1_feedback, g1(1:9)];
% 更新G2寄存器
g2_feedback = mod(g2(2) + g2(3) + g2(6) + g2(8) + g2(9) + g2(10), 2);
g2 = [g2_feedback, g2(1:9)];
end
ca_code(ca_code == 0) = -1; % 转换为±1
end
function tap = get_g2_tap(~, prn)
% 获取G2寄存器的抽头位置
% 不同PRN对应不同的G2延迟组合
g2_taps = [
2, 6; % PRN 1
3, 7; % PRN 2
4, 8; % PRN 3
5, 9; % PRN 4
1, 9; % PRN 5
2, 10; % PRN 6
1, 8; % PRN 7
2, 9; % PRN 8
3, 10; % PRN 9
2, 3; % PRN 10
% ... 其他PRN
];
if prn <= size(g2_taps, 1)
tap = g2_taps(prn, 2); % 使用第二个抽头
else
tap = 2; % 默认值
end
end
end
end
三、接收机信号处理仿真
1. 信号捕获算法
%% 北斗/GPS信号捕获仿真
classdef GNSSAcquisition
properties
fs = 16.368e6; % 采样频率 (Hz)
fc = 1561.098e6; % 载波频率 (Hz)
code_rate = 2.046e6; % 码速率 (chips/s)
code_length = 2046; % 码长
doppler_search = -10000:500:10000; % 多普勒搜索范围 (Hz)
end
methods
function [doppler_est, code_phase_est, correlation] = acquire_signal(obj, signal, prn)
% 并行码相位搜索捕获算法
% 1. 生成本地码
local_code = obj.generate_local_code(prn);
local_code = repmat(local_code, 1, ceil(length(signal)/length(local_code)));
local_code = local_code(1:length(signal));
% 2. 准备搜索矩阵
num_doppler_bins = length(obj.doppler_search);
correlation = zeros(num_doppler_bins, obj.code_length);
% 3. 并行多普勒搜索
for d_idx = 1:num_doppler_bins
doppler = obj.doppler_search(d_idx);
% 生成本地载波
t = (0:length(signal)-1)' / obj.fs;
local_carrier = exp(-1i*2*pi*(obj.fc + doppler)*t);
% 下变频
baseband_signal = signal .* local_carrier;
% 并行码相位搜索
for phase = 1:obj.code_length
% 码相位对齐
shifted_code = circshift(local_code, phase-1);
% 相关计算
correlation(d_idx, phase) = abs(sum(baseband_signal .* shifted_code));
end
end
% 4. 寻找峰值
[max_val, max_idx] = max(correlation(:));
[doppler_idx, phase_idx] = ind2sub(size(correlation), max_idx);
doppler_est = obj.doppler_search(doppler_idx);
code_phase_est = phase_idx;
% 5. 计算检测统计量
mean_corr = mean(correlation(:));
std_corr = std(correlation(:));
detection_metric = (max_val - mean_corr) / std_corr;
fprintf('捕获结果:\n');
fprintf(' 估计多普勒: %.2f Hz\n', doppler_est);
fprintf(' 估计码相位: %d chips\n', code_phase_est);
fprintf(' 检测统计量: %.2f\n', detection_metric);
fprintf(' 峰值相关值: %.2f\n', max_val);
if detection_metric > 3.0
fprintf(' 信号检测: 成功\n');
else
fprintf(' 信号检测: 失败\n');
end
end
function plot_acquisition_results(obj, correlation, doppler_est, code_phase_est)
% 绘制捕获结果
figure('Position', [100 100 1200 400]);
% 3D相关平面
subplot(1, 3, 1);
[X, Y] = meshgrid(1:obj.code_length, obj.doppler_search);
surf(X, Y/1000, correlation, 'EdgeColor', 'none');
xlabel('码相位 (chips)'); ylabel('多普勒频率 (kHz)'); zlabel('相关值');
title('捕获相关平面');
colorbar; view(45, 30);
% 多普勒剖面
subplot(1, 3, 2);
doppler_idx = find(obj.doppler_search == doppler_est);
plot(obj.doppler_search/1000, correlation(:, code_phase_est), 'b-', 'LineWidth', 2);
hold on;
plot(doppler_est/1000, correlation(doppler_idx, code_phase_est), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('多普勒频率 (kHz)'); ylabel('相关值');
title('多普勒频率剖面');
grid on;
legend('相关曲线', '峰值点');
% 码相位剖面
subplot(1, 3, 3);
plot(1:obj.code_length, correlation(doppler_idx, :), 'b-', 'LineWidth', 2);
hold on;
plot(code_phase_est, correlation(doppler_idx, code_phase_est), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('码相位 (chips)'); ylabel('相关值');
title('码相位剖面');
grid on;
legend('相关曲线', '峰值点');
end
function local_code = generate_local_code(~, prn)
% 生成本地码（简化版）
% 实际实现需要根据具体信号类型生成
code_length = 2046; % 北斗B1I码长
local_code = 2*randi([0 1], 1, code_length) - 1;
end
end
end
2. 信号跟踪环路
%% 北斗/GPS信号跟踪仿真
classdef GNSSTracking
properties
fs = 16.368e6; % 采样频率 (Hz)
code_rate = 2.046e6; % 码速率 (chips/s)
integration_time = 0.001; % 积分时间 (s)
% 锁相环参数
pll_bandwidth = 15; % PLL带宽 (Hz)
pll_damping = 0.707; % PLL阻尼系数
% 延迟锁定环参数
dll_bandwidth = 2; % DLL带宽 (Hz)
dll_damping = 0.707; % DLL阻尼系数
early_late_spacing = 0.5; % 早晚码间距 (chips)
end
methods
function [carrier_phase, code_phase, measurements] = track_signal(obj, signal, init_doppler, init_code_phase, prn)
% 信号跟踪主函数
% 初始化
num_samples = length(signal);
samples_per_integration = round(obj.fs * obj.integration_time);
num_integrations = floor(num_samples / samples_per_integration);
carrier_phase = zeros(num_integrations, 1);
code_phase = zeros(num_integrations, 1);
measurements.doppler = zeros(num_integrations, 1);
measurements.code_freq = zeros(num_integrations, 1);
measurements.I = zeros(num_integrations, 1);
measurements.Q = zeros(num_integrations, 1);
measurements.cn0 = zeros(num_integrations, 1);
% 初始状态
current_carrier_phase = 0;
current_code_phase = init_code_phase;
current_doppler = init_doppler;
current_code_freq = obj.code_rate * (1 + current_doppler/1575.42e6);
% 生成本地码
local_code = obj.generate_local_code(prn);
% 跟踪循环
for k = 1:num_integrations
% 获取当前积分段数据
start_idx = (k-1)*samples_per_integration + 1;
end_idx = k*samples_per_integration;
signal_segment = signal(start_idx:end_idx);
% 生成本地信号
[I, Q, prompt_code, early_code, late_code] = ...
obj.generate_local_signals(signal_segment, current_carrier_phase, ...
current_code_phase, current_code_freq, local_code);
% 计算鉴别器输出
pll_error = obj.pll_discriminator(I, Q);
dll_error = obj.dll_discriminator(I, Q, early_code, late_code, prompt_code);
% 更新环路滤波器
[current_doppler, current_carrier_phase] = ...
obj.pll_update(pll_error, current_doppler, current_carrier_phase);
[current_code_freq, current_code_phase] = ...
obj.dll_update(dll_error, current_code_freq, current_code_phase);
% 保存结果
carrier_phase(k) = current_carrier_phase;
code_phase(k) = current_code_phase;
measurements.doppler(k) = current_doppler;
measurements.code_freq(k) = current_code_freq;
measurements.I(k) = I;
measurements.Q(k) = Q;
measurements.cn0(k) = obj.estimate_cn0(I, Q);
% 更新载波相位
current_carrier_phase = mod(current_carrier_phase + ...
2*pi*current_doppler*obj.integration_time, 2*pi);
% 更新码相位
current_code_phase = mod(current_code_phase + ...
current_code_freq*obj.integration_time, length(local_code));
end
% 绘制跟踪结果
obj.plot_tracking_results(carrier_phase, code_phase, measurements);
end
function [I, Q, prompt_code, early_code, late_code] = generate_local_signals(...
obj, signal_segment, carrier_phase, code_phase, code_freq, local_code)
% 生成本地参考信号
num_samples = length(signal_segment);
t = (0:num_samples-1)' / obj.fs;
% 生成载波
local_carrier_i = cos(2*pi*carrier_phase*t);
local_carrier_q = -sin(2*pi*carrier_phase*t);
% 生成码序列
code_samples = round(obj.fs / code_freq);
code_indices = mod(round(code_phase + (0:num_samples-1)' * code_freq / obj.fs), ...
length(local_code)) + 1;
prompt_code = local_code(code_indices);
early_code = local_code(mod(code_indices - round(obj.early_late_spacing*code_samples), ...
length(local_code)) + 1);
late_code = local_code(mod(code_indices + round(obj.early_late_spacing*code_samples), ...
length(local_code)) + 1);
% 相关计算
I = sum(signal_segment .* prompt_code .* local_carrier_i);
Q = sum(signal_segment .* prompt_code .* local_carrier_q);
end
function error = pll_discriminator(~, I, Q)
% PLL鉴别器（二象限反正切）
error = atan2(Q, I);
end
function error = dll_discriminator(~, I, Q, early_code, late_code, prompt_code)
% DLL鉴别器（非相干早迟包络）
early_power = I^2 + Q^2; % 简化计算
late_power = I^2 + Q^2; % 实际需要分别计算早码和迟码的相关
error = (early_power - late_power) / (early_power + late_power);
end
function [new_doppler, new_phase] = pll_update(obj, error, current_doppler, current_phase)
% PLL环路滤波器更新
% 简化的一阶环路滤波器
pll_gain = 4 * obj.pll_bandwidth / (obj.pll_damping + 1/(4*obj.pll_damping));
new_doppler = current_doppler + pll_gain * error;
new_phase = current_phase;
end
function [new_code_freq, new_code_phase] = dll_update(obj, error, current_code_freq, current_code_phase)
% DLL环路滤波器更新
dll_gain = 4 * obj.dll_bandwidth / (obj.dll_damping + 1/(4*obj.dll_damping));
new_code_freq = current_code_freq + dll_gain * error;
new_code_phase = current_code_phase;
end
function cn0 = estimate_cn0(~, I, Q)
% 估计载噪比
signal_power = I^2 + Q^2;
noise_power = (I^2 + Q^2) / 2; % 简化估计
cn0 = 10 * log10(signal_power / noise_power);
end
function plot_tracking_results(obj, carrier_phase, code_phase, measurements)
% 绘制跟踪结果
time = (1:length(carrier_phase))' * obj.integration_time;
figure('Position', [100 100 1200 800]);
% 载波相位跟踪
subplot(3, 2, 1);
plot(time, unwrap(carrier_phase), 'b-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('载波相位 (rad)');
title('载波相位跟踪');
grid on;
% 码相位跟踪
subplot(3, 2, 2);
plot(time, code_phase, 'r-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('码相位 (chips)');
title('码相位跟踪');
grid on;
% 多普勒频率
subplot(3, 2, 3);
plot(time, measurements.doppler/1000, 'g-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('多普勒频率 (kHz)');
title('多普勒频率估计');
grid on;
% I/Q支路
subplot(3, 2, 4);
plot(time, measurements.I, 'b-', 'LineWidth', 1.5);
hold on;
plot(time, measurements.Q, 'r-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('幅度');
title('I/Q支路输出');
legend('I支路', 'Q支路');
grid on;
% 星座图
subplot(3, 2, 5);
scatter(measurements.I, measurements.Q, 20, time, 'filled');
xlabel('I支路'); ylabel('Q支路');
title('跟踪星座图');
axis equal; grid on;
colorbar;
% 载噪比估计
subplot(3, 2, 6);
plot(time, measurements.cn0, 'm-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('C/N0 (dB-Hz)');
title('载噪比估计');
grid on;
ylim([30 50]);
end
end
end
参考代码 北斗GPS接收机，包含信号仿真 www.youwenfan.com/contentcst/63146.html
四、完整仿真系统集成
%% 北斗/GPS接收机完整仿真系统
classdef GNSSReceiverSimulator
properties
% 仿真参数
fs = 16.368e6; % 采样频率
duration = 0.1; % 仿真时长 (s)
num_satellites = 6; % 卫星数量
% 卫星参数
satellites = struct();
% 接收机模块
acquisition_module;
tracking_module;
navigation_module;
end
methods
function obj = GNSSReceiverSimulator()
% 初始化仿真系统
% 初始化卫星参数
obj.initialize_satellites();
% 初始化接收机模块
obj.acquisition_module = GNSSAcquisition();
obj.tracking_module = GNSSTracking();
obj.navigation_module = GNSSNavigation();
end
function initialize_satellites(obj)
% 初始化卫星参数
for i = 1:obj.num_satellites
obj.satellites(i).prn = i;
obj.satellites(i).cn0 = 45 + randn() * 2; % 载噪比
obj.satellites(i).doppler = 5000 + randn() * 1000; % 多普勒
obj.satellites(i).code_phase = randi([1 2046]); % 码相位
obj.satellites(i).carrier_phase = rand() * 2*pi; % 载波相位
obj.satellites(i).elevation = 30 + rand() * 60; % 仰角
obj.satellites(i).azimuth = rand() * 360; % 方位角
end
end
function [composite_signal, satellite_signals] = generate_composite_signal(obj)
% 生成复合信号（多颗卫星）
num_samples = round(obj.fs * obj.duration);
composite_signal = zeros(num_samples, 1);
satellite_signals = cell(obj.num_satellites, 1);
% 生成每颗卫星的信号
for i = 1:obj.num_satellites
sat = obj.satellites(i);
% 创建信号生成器
if mod(i, 2) == 0
% 偶数PRN使用北斗
generator = BeiDouB1ISignalGenerator();
generator.prn = sat.prn;
generator.cn0 = sat.cn0;
generator.doppler = sat.doppler;
generator.code_phase = sat.code_phase;
generator.carrier_phase = sat.carrier_phase;
else
% 奇数PRN使用GPS
generator = GPSL1CASignalGenerator();
generator.prn = sat.prn;
generator.cn0 = sat.cn0;
generator.doppler = sat.doppler;
generator.code_phase = sat.code_phase;
generator.carrier_phase = sat.carrier_phase;
end
generator.fs = obj.fs;
generator.duration = obj.duration;
% 生成单星信号
[signal, ~] = generator.generate_signal();
satellite_signals{i} = signal;
% 叠加到复合信号
composite_signal = composite_signal + signal;
end
% 添加噪声
noise_power = var(composite_signal) / (10^(40/10)); % 40 dB SNR
noise = sqrt(noise_power) * randn(size(composite_signal));
composite_signal = composite_signal + noise;
% 归一化
composite_signal = composite_signal / max(abs(composite_signal));
end
function run_simulation(obj)
% 运行完整仿真
fprintf('=== 北斗/GPS接收机仿真开始 ===\n\n');
% 1. 生成复合信号
fprintf('1. 生成复合信号...\n');
[composite_signal, satellite_signals] = obj.generate_composite_signal();
% 2. 信号捕获
fprintf('\n2. 信号捕获...\n');
acquisition_results = cell(obj.num_satellites, 1);
for i = 1:obj.num_satellites
fprintf(' 捕获卫星 PRN%d...\n', i);
[doppler_est, code_phase_est, correlation] = ...
obj.acquisition_module.acquire_signal(composite_signal, i);
acquisition_results{i}.doppler_est = doppler_est;
acquisition_results{i}.code_phase_est = code_phase_est;
acquisition_results{i}.correlation = correlation;
% 绘制捕获结果
if i <= 2 % 只绘制前两颗卫星的结果
figure('Name', sprintf('卫星PRN%d捕获结果', i));
obj.acquisition_module.plot_acquisition_results(...
correlation, doppler_est, code_phase_est);
end
end
% 3. 信号跟踪
fprintf('\n3. 信号跟踪...\n');
tracking_results = cell(obj.num_satellites, 1);
for i = 1:min(4, obj.num_satellites) % 只跟踪前4颗卫星
fprintf(' 跟踪卫星 PRN%d...\n', i);
[carrier_phase, code_phase, measurements] = ...
obj.tracking_module.track_signal(...
composite_signal, ...
acquisition_results{i}.doppler_est, ...
acquisition_results{i}.code_phase_est, ...
i);
tracking_results{i}.carrier_phase = carrier_phase;
tracking_results{i}.code_phase = code_phase;
tracking_results{i}.measurements = measurements;
end
% 4. 定位解算
fprintf('\n4. 定位解算...\n');
position = obj.navigation_module.calculate_position(tracking_results);
fprintf('\n=== 仿真完成 ===\n');
fprintf('估计位置: 纬度 %.6f°, 经度 %.6f°, 高度 %.2f m\n', ...
position.latitude, position.longitude, position.altitude);
fprintf('位置误差: %.2f m\n', position.error);
end
end
end
%% 导航定位解算模块
classdef GNSSNavigation
methods
function position = calculate_position(~, tracking_results)
% 简化的定位解算
% 实际实现需要卫星星历、伪距测量等
% 这里使用模拟数据
position.latitude = 30.0000 + randn() * 0.0001;
position.longitude = 120.0000 + randn() * 0.0001;
position.altitude = 50 + randn() * 10;
position.error = 5 + rand() * 3;
% 实际定位算法应包括：
% 1. 伪距测量
% 2. 卫星位置计算
% 3. 最小二乘或卡尔曼滤波定位
% 4. 误差修正（电离层、对流层、钟差等）
end
end
end
%% 主仿真程序
% 运行完整仿真
clear; clc; close all;
% 创建仿真器
simulator = GNSSReceiverSimulator();
% 运行仿真
simulator.run_simulation();
% 信号分析
[composite_signal, satellite_signals] = simulator.generate_composite_signal();
% 绘制信号频谱
figure('Position', [100 100 800 600]);
subplot(2, 2, 1);
plot(real(composite_signal(1:1000)));
xlabel('采样点'); ylabel('幅度');
title('复合信号时域波形（前1000点）');
grid on;
subplot(2, 2, 2);
[pxx, f] = pwelch(composite_signal, 1024, 512, 1024, simulator.fs, 'centered');
plot(f/1e6, 10*log10(pxx));
xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB/Hz)');
title('复合信号功率谱');
grid on;
subplot(2, 2, 3);
histogram(real(composite_signal), 50, 'Normalization', 'pdf');
xlabel('幅度'); ylabel('概率密度');
title('信号幅度分布');
grid on;
subplot(2, 2, 4);
spectrogram(composite_signal(1:min(100000, length(composite_signal))), ...
256, 250, 256, simulator.fs, 'yaxis');
title('信号频谱图');
colorbar;
五、MATLAB工具与资源
1. 官方工具箱
- Communications Toolbox：提供信号生成、调制解调功能
- Signal Processing Toolbox：提供频谱分析、滤波器设计
- Phased Array System Toolbox：提供阵列信号处理
- Satellite Communications Toolbox：专门用于卫星通信仿真
2. 第三方工具包
- GNSS-SDR：开源的GNSS软件定义接收机
- GPSTk：GPS工具包
- RTKLIB：开源的GNSS定位库
3. 仿真建议
- 信号生成：从简单到复杂，先仿真单星信号，再仿真多星复合信号
- 算法验证：使用已知参数的信号验证捕获和跟踪算法
- 性能评估：在不同信噪比、动态场景下测试接收机性能
- 实时性考虑：优化算法以满足实时处理要求
4. 实际应用注意事项
- 信号特性：北斗和GPS信号有不同特性，需要分别处理
- 干扰抑制：实际环境中存在多径、干扰等问题
- 动态适应：高动态场景需要更宽的捕获范围和更快的跟踪环路
- 多系统兼容：现代接收机通常支持多系统（GPS、北斗、GLONASS、Galileo）
