（期末课程设计
基于VPI或matlab仿真软件实现单偏振
28GBaudQPSK信号传输2000km;
1.基于matlab实现色散补偿算法、教波相
位恢复算法：直接运行的模型代码（必须
含注释）
2.传输系统：单跨80km、EDFA噪声指数
5dB.）帮我完成这个课程设计，目前已有代码，请帮我看看有没有什么问题，哪里需要改成等等（%% =============================================================
% 期末课程设计（最终整合版，可直接运行，单文件提交）
% 单偏振 28 Gbaud QPSK 信号传输 2000 km（25×80 km）
%
% 要求：
% 1) MATLAB 实现色散补偿算法（EDC）+ 载波相位恢复算法（V&V）
% 2) 传输系统：单跨 80 km，EDFA 噪声指数 NF=5 dB
%
% 本代码包含：
% Tx: QPSK(Gray) + RRC 成形 + 激光相位噪声(线宽 linewidth)
% Link: 25×80 km SMF 色散 + 损耗 + EDFA补偿 + ASE噪声
% Rx: 频域色散补偿(EDC) + RRC 匹配滤波 + 理想抽样
% CPR: Viterbi&Viterbi (QPSK): 四次方 + unwrap + 平滑 + 去π/4偏置
% QPSK π/2 相位模糊：仿真中用“已知比特”选最优旋转（等效导频）
%
% 额外：可选 sweep（linewidth/NF）+ MonteCarlo + BER=0 的 UB95 输出
% =============================================================
clear; clc; close all;
rng(1);
%% --------------------- 开关区（按需改） ---------------------
doPlots_singlePoint = true; % 是否画“单点仿真”的星座/相位图
plot_noCD_constellation = true; % 是否画“未做CD补偿”的星座图（建议保留）
doSweep_linewidth = true; % 是否扫 linewidth
doSweep_NF = true; % 是否扫 NF
Nmc_sweep = 10; % sweep 时每个点做 Monte Carlo 次数（建议 >=5）
%% --------------------- 系统参数 -----------------------------
Rs = 28e9; % 符号率 28 Gbaud
sps = 8; % 每符号采样数
Fs = Rs*sps; % 采样率
Ts = 1/Fs;
Nsym_single = 2^17; % 单点仿真符号数
Nsym_sweep = 2^15; % sweep 时每次仿真符号数（为加速）
% RRC 滤波（成形/匹配）
rolloff = 0.2;
span_sym = 6;
rrc = rcosdesign(rolloff, span_sym, sps, 'sqrt');
rrc = rrc / sqrt(sum(rrc.^2)); % 单位能量归一化
rrcDelay = (length(rrc)-1)/2; % 群时延（采样）
%% --------------------- 光纤参数（只考虑线性：色散+损耗） -------
lambda = 1550e-9;
c = 3e8;
D_ps_nm_km = 16; % 典型 SMF
D_SI = D_ps_nm_km * 1e-12 / (1e-9*1e3); % s/m^2
beta2 = -(D_SI*lambda^2)/(2*pi*c); % s^2/m
L_span_km = 80;
L_span = L_span_km*1e3;
N_span = 2000 / L_span_km; % 25 跨
L_total = L_span*N_span;
alpha_dB_km = 0.2; % dB/km（常见）
lossSpan_dB = alpha_dB_km * L_span_km;
A_field = 10^(-lossSpan_dB/20); % 场衰减
G_span_dB = lossSpan_dB; % EDFA 每跨补偿损耗
G_span = 10^(G_span_dB/10);
%% --------------------- 噪声（ASE）相关 -----------------------
h = 6.626e-34;
nu = c/lambda;
% 课程要求：NF=5 dB（单点仿真固定用它；sweep 时会变）
NF0_dB = 5;
%% --------------------- 单点仿真参数 --------------------------
linewidth0 = 10e3; % 10 kHz
NF_dB0 = NF0_dB;
% V&V 平滑窗口：可“固定”也可“随 linewidth 自适应”
useAutoM = true; % linewidth 提高时自动减小窗口，避免跟踪不动
M_fixed = 501; % 不自适应时用这个
% 自适应窗口范围（防止过大/过小）
M_min = 51;
M_max = 2001;
%% =============================================================
% 单点仿真（用于展示：星座图/相位图/打印结果）
%% =============================================================
fprintf('\n=========== 单点仿真（用于画图/展示） ===========\n');
[ber0, numErr0, ub95_0, bestK0, M_used0, dbg0] = simulate_once( ...
Rs, sps, Fs, Ts, Nsym_single, rrc, rrcDelay, ...
beta2, L_span, N_span, A_field, G_span, ...
h, nu, linewidth0, NF_dB0, ...
useAutoM, M_fixed, M_min, M_max, ...
doPlots_singlePoint, plot_noCD_constellation);
fprintf('---------------- 仿真结果 ----------------\n');
fprintf('传输距离：%.0f km (%d × %.0f km)\n', L_total/1e3, N_span, L_span_km);
fprintf('QPSK 符号数：%d\n', Nsym_single);
fprintf('激光线宽：%.1f kHz\n', linewidth0/1e3);
fprintf('EDFA NF：%.1f dB\n', NF_dB0);
fprintf('V&V 窗口 M=%d（useAutoM=%d）\n', M_used0, useAutoM);
fprintf('相位模糊旋转：k=%d (%d°)\n', bestK0, bestK0*90);
fprintf('误比特数：%d\n', numErr0);
if numErr0==0
fprintf('BER ≈ 0.000e+00（若为0，看 95%%上界 UB95≈%.3e）\n', ub95_0);
else
fprintf('BER ≈ %.3e\n', ber0);
end
fprintf('【单点】err=%d, BER=%.3e, UB95=%.3e, bitsUsed=%d\n', ...
numErr0, ber0, ub95_0, dbg0.bitsPerRun);
%% =============================================================
% 扫 linewidth（固定 NF=5 dB）
%% =============================================================
if doSweep_linewidth
fprintf('\n=========== 扫 linewidth（固定 NF=%.1f dB, Nmc=%d） ===========\n', NF0_dB, Nmc_sweep);
linewidth_list = [ ...
1e3 2e3 5e3 10e3 20e3 50e3 100e3 200e3 500e3 1e6 ]; % Hz
ber_lw = zeros(size(linewidth_list));
ub95_lw = zeros(size(linewidth_list));
for i = 1:numel(linewidth_list)
lw = linewidth_list(i);
totalErr = 0;
totalBits = 0;
for mc = 1:Nmc_sweep
[~, err_i, ~, ~, ~, dbg_i] = simulate_once( ...
Rs, sps, Fs, Ts, Nsym_sweep, rrc, rrcDelay, ...
beta2, L_span, N_span, A_field, G_span, ...
h, nu, lw, NF0_dB, ...
useAutoM, M_fixed, M_min, M_max, ...
false, false);
totalErr = totalErr + err_i;
totalBits = totalBits + dbg_i.bitsPerRun;
end
if totalErr == 0
ber_lw(i) = 0;
ub95_lw(i) = -log(0.05) / totalBits;
else
ber_lw(i) = totalErr / totalBits;
ub95_lw(i) = ber_lw(i);
end
fprintf('lw=%7.1f kHz | totalErr=%6d | totalBits=%d | BER=%.3e | UB95=%.3e\n', ...
lw/1e3, totalErr, totalBits, ber_lw(i), ub95_lw(i));
end
figure;
semilogy(linewidth_list/1e3, max(ub95_lw, 1e-12), '-o');
grid on;
xlabel('Laser linewidth (kHz)');
ylabel('BER (0 -> UB95)');
title(sprintf('BER vs Linewidth (NF=%.1f dB, Nsym=%d, Nmc=%d)', NF0_dB, Nsym_sweep, Nmc_sweep));
end
%% =============================================================
% 扫 NF（固定 linewidth=10 kHz）
%% =============================================================
if doSweep_NF
fprintf('\n=========== 扫 NF（固定 linewidth=%.1f kHz, Nmc=%d） ===========\n', linewidth0/1e3, Nmc_sweep);
NF_list = 3:1:15;
ber_nf = zeros(size(NF_list));
ub95_nf = zeros(size(NF_list));
for i = 1:numel(NF_list)
NF_dB = NF_list(i);
totalErr = 0;
totalBits = 0;
for mc = 1:Nmc_sweep
[~, err_i, ~, ~, ~, dbg_i] = simulate_once( ...
Rs, sps, Fs, Ts, Nsym_sweep, rrc, rrcDelay, ...
beta2, L_span, N_span, A_field, G_span, ...
h, nu, linewidth0, NF_dB, ...
useAutoM, M_fixed, M_min, M_max, ...
false, false);
totalErr = totalErr + err_i;
totalBits = totalBits + dbg_i.bitsPerRun;
end
if totalErr == 0
ber_nf(i) = 0;
ub95_nf(i) = -log(0.05) / totalBits;
else
ber_nf(i) = totalErr / totalBits;
ub95_nf(i) = ber_nf(i);
end
fprintf('NF=%4.1f dB | totalErr=%6d | totalBits=%d | BER=%.3e | UB95=%.3e\n', ...
NF_dB, totalErr, totalBits, ber_nf(i), ub95_nf(i));
end
figure;
semilogy(NF_list, max(ub95_nf, 1e-12), '-o');
grid on;
xlabel('EDFA Noise Figure NF (dB)');
ylabel('BER (0 -> UB95)');
title(sprintf('BER vs NF (linewidth=%.1f kHz, Nsym=%d, Nmc=%d)', linewidth0/1e3, Nsym_sweep, Nmc_sweep));
end
%% =============================================================
% 子函数：单次仿真（写在同一个脚本末尾，便于“单文件提交”）
%% =============================================================
function [ber, numErr, ub95, bestK, M_used, dbg] = simulate_once( ...
Rs, sps, Fs, Ts, Nsym, rrc, rrcDelay, ...
beta2, L_span, N_span, A_field, G_span, ...
h, nu, linewidth, NF_dB, ...
useAutoM, M_fixed, M_min, M_max, ...
doPlots, plotNoCD)
%% ---------- EDFA ASE ----------
NF = 10^(NF_dB/10);
n_sp = NF/2; % 常用近似：NF ≈ 2*nsp
S_ASE = n_sp*h*nu*(G_span-1); % 单偏振、单边 PSD (W/Hz)
%% ---------- 生成 QPSK ----------
bits = randi([0 1], 2*Nsym, 1);
bitI = bits(1:2:end);
bitQ = bits(2:2:end);
txSym = ((1-2*bitI) + 1j*(1-2*bitQ))/sqrt(2); % Gray（I/Q独立映射）
txWave = upfirdn(txSym, rrc, sps, 1);
txWave = txWave(:);
txWave = txWave / rms(txWave);
%% ---------- 激光相位噪声（Wiener 过程） ----------
phi = sqrt(2*pi*linewidth*Ts) * cumsum(randn(size(txWave)));
txWave = txWave .* exp(1j*phi);
%% ---------- 频率轴 ----------
N = length(txWave);
df = Fs/N;
f = (-N/2:N/2-1).' * df;
%% ---------- 单跨色散 ----------
H_cd_span = exp(-1j*2*pi^2*beta2*L_span*(f.^2));
% 噪声离散化：对每个频点（带宽 df）加入复高斯噪声
% 若 PSD 为单边 S_ASE，则离散到每个 FFT bin 的功率 ~ S_ASE * df
noiseVar = S_ASE * df; % 复噪声功率 E{|n|^2}
noiseStd = sqrt(noiseVar/2); % 每个实/虚分量方差 = noiseVar/2
%% ---------- 链路：N_span 跨 ----------
rxWave = txWave;
for n = 1:N_span
RX = fftshift(fft(rxWave));
RX = RX .* H_cd_span;
rxWave = ifft(ifftshift(RX));
rxWave = rxWave * A_field;
rxWave = rxWave * sqrt(G_span);
rxWave = rxWave + noiseStd*(randn(size(rxWave))+1j*randn(size(rxWave)));
end
%% ---------- EDC（频域色散补偿） ----------
L_total = L_span*N_span;
H_cd_total = exp(-1j*2*pi^2*beta2*L_total*(f.^2));
H_eq = conj(H_cd_total);
RX = fftshift(fft(rxWave));
rxWave_cd = ifft(ifftshift(RX .* H_eq));
%% ---------- 匹配滤波 + 抽样 ----------
rxMF_cd = upfirdn(rxWave_cd, rrc, 1, 1);
rxMF_cd = rxMF_cd / rms(rxMF_cd);
if plotNoCD
rxMF_noCD = upfirdn(rxWave, rrc, 1, 1);
rxMF_noCD = rxMF_noCD / rms(rxMF_noCD);
end
startIdx = 2*rrcDelay + 1;
idx = startIdx : sps : (startIdx + (Nsym-1)*sps);
if idx(end) > length(rxMF_cd)
error('接收波形长度不足：请减小 Nsym 或增加零填充。');
end
rxSym = rxMF_cd(idx);
if plotNoCD
if idx(end) > length(rxMF_noCD)
error('接收波形长度不足（noCD）：请减小 Nsym 或增加零填充。');
end
rxSym_noCD = rxMF_noCD(idx);
end
%% ---------- V&V CPR ----------
if useAutoM
% 经验规则：linewidth 越大 -> M 越小（跟踪更快）
K = 50; % K 越大 -> M 越小
M_used = round(Rs/(K*max(linewidth,1)));
if mod(M_used,2)==0, M_used = M_used+1; end
M_used = min(max(M_used, M_min), M_max);
else
M_used = M_fixed;
if mod(M_used,2)==0, M_used = M_used+1; end
end
r4 = rxSym.^4;
phi4 = unwrap(angle(r4)); % = 4*phi_carrier + noise + 4*data_phase(被消掉)
if exist('movmean','file') == 2
phi4_s = movmean(phi4, M_used);
else
w = ones(M_used,1)/M_used;
phi4_s = conv(phi4, w, 'same');
end
phi_est = phi4_s/4 - pi/4; % 去掉QPSK固有 π/4 偏置
rxSym_cpr = rxSym .* exp(-1j*phi_est);
%% ---------- π/2 相位模糊：用“已知比特”选择最优旋转 ----------
valid = M_used:Nsym; % 丢弃前 M-1 个（V&V 边界不稳定）
txI = bitI(valid);
txQ = bitQ(valid);
y = rxSym_cpr(valid);
bestErr = inf;
bestK = 0;
bestY = y;
for k = 0:3
yk = y .* exp(-1j*k*pi/2);
rxI = real(yk) < 0;
rxQ = imag(yk) < 0;
errTotal = sum(rxI ~= txI) + sum(rxQ ~= txQ);
if errTotal < bestErr
bestErr = errTotal;
bestK = k;
bestY = yk;
end
end
rxSym_det = bestY;
%% ---------- BER ----------
rxI = real(rxSym_det) < 0;
rxQ = imag(rxSym_det) < 0;
txBits = zeros(2*length(valid),1);
rxBits = zeros(2*length(valid),1);
txBits(1:2:end) = txI; txBits(2:2:end) = txQ;
rxBits(1:2:end) = rxI; rxBits(2:2:end) = rxQ;
numErr = sum(txBits ~= rxBits);
nBits = length(txBits);
ber = numErr / nBits;
if numErr == 0
ub95 = -log(0.05) / nBits; % 95% one-sided upper bound
else
ub95 = ber;
end
dbg.bitsPerRun = nBits;
%% ---------- 画图 ----------
if doPlots
figure;
plot(real(txSym(1:1000)), imag(txSym(1:1000)), 'o');
axis square; grid on;
title('发送端 QPSK 星座（前 1000 符号）');
xlabel('I'); ylabel('Q');
if plotNoCD
figure;
plot(real(rxSym_noCD(1:2000)), imag(rxSym_noCD(1:2000)), '.');
axis square; grid on;
title('接收端星座（未做CD补偿，前 2000 符号）');
xlabel('I'); ylabel('Q');
end
figure;
plot(real(rxSym(1:2000)), imag(rxSym(1:2000)), '.');
axis square; grid on;
title('接收端星座（EDC 后，CPR 前）');
xlabel('I'); ylabel('Q');
figure;
plot(real(rxSym_det(1:min(2000,length(rxSym_det)))), imag(rxSym_det(1:min(2000,length(rxSym_det)))), '.');
axis square; grid on;
title('接收端星座（EDC + CPR 后）');
xlabel('I'); ylabel('Q');
% =======================
% Figure 5（更像教材的 V&V 相位跟踪图）
% =======================
Kplot = min(200, length(rxSym));
% “教材式”中通常展示：
% φ_raw = unwrap(angle(r^4))/4
% φ_smooth= movmean(φ_raw)
% φ_est = φ_smooth - π/4
% φ_res = unwrap(angle((r*e^{-jφ_est})^4))/4 (残余更小更平)
phi_raw = unwrap(angle(rxSym(1:Kplot).^4))/4;
if exist('movmean','file') == 2
phi_smooth = movmean(phi_raw, min(M_used,Kplot));
else
w = ones(min(M_used,Kplot),1)/min(M_used,Kplot);
phi_smooth = conv(phi_raw, w, 'same');
end
phi_est_p = phi_smooth - pi/4;
r_cpr_p = rxSym(1:Kplot).*exp(-1j*phi_est_p);
phi_res = unwrap(angle(r_cpr_p.^4))/4;
% 为了更直观：把曲线整体去均值（只看“变化量/抖动”）
phi_raw_z = phi_raw - mean(phi_raw);
phi_smooth_z = phi_smooth - mean(phi_smooth);
phi_est_z = phi_est_p - mean(phi_est_p);
phi_res_z = phi_res - mean(phi_res);
figure;
plot(1:Kplot, phi_raw_z, '-'); hold on;
plot(1:Kplot, phi_smooth_z, '-');
plot(1:Kplot, phi_est_z, '-');
plot(1:Kplot, phi_res_z, '-');
grid on;
xlabel('符号索引 k');
ylabel('相位 [rad]');
title('教材式载波相位跟踪示意（V&V：四次方去数据相位）');
legend( ...
'\phi_{raw} = unwrap(\angle(r^4))/4 (去数据后)', ...
'\phi_{smooth} = movmean(\phi_{raw})', ...
'\phi_{est} = \phi_{smooth} - \pi/4 (用于补偿)', ...
'\phi_{res} = unwrap(\angle((re^{-j\phi_{est}})^4))/4 (补偿后残余)', ...
'Location','best');
end
end
）已有运行结果为：
=========== 单点仿真（用于画图/展示） ===========
---------------- 仿真结果 ----------------
传输距离：2000 km (25 × 80 km)
QPSK 符号数：131072
激光线宽：10.0 kHz
EDFA NF：5.0 dB
V&V 窗口 M=2001（useAutoM=1）
相位模糊旋转：k=1 (90°)
误比特数：0
BER ≈ 0.000e+00（若为0，看 95%上界 UB95≈1.160e-05）
【单点】err=0, BER=0.000e+00, UB95=1.160e-05, bitsUsed=258144
=========== 扫 linewidth（固定 NF=5.0 dB, Nmc=10） ===========
lw= 1.0 kHz | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
lw= 2.0 kHz | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
lw= 5.0 kHz | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
lw= 10.0 kHz | totalErr= 2 | totalBits=615360 | BER=3.250e-06 | UB95=3.250e-06
lw= 20.0 kHz | totalErr= 3 | totalBits=615360 | BER=4.875e-06 | UB95=4.875e-06
lw= 50.0 kHz | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
lw= 100.0 kHz | totalErr= 840 | totalBits=615360 | BER=1.365e-03 | UB95=1.365e-03
lw= 200.0 kHz | totalErr= 1806 | totalBits=615360 | BER=2.935e-03 | UB95=2.935e-03
lw= 500.0 kHz | totalErr= 9 | totalBits=632960 | BER=1.422e-05 | UB95=1.422e-05
lw= 1000.0 kHz | totalErr= 766 | totalBits=644160 | BER=1.189e-03 | UB95=1.189e-03
=========== 扫 NF（固定 linewidth=10.0 kHz, Nmc=10） ===========
NF= 3.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
NF= 4.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
NF= 5.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
NF= 6.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
NF= 7.0 dB | totalErr= 6 | totalBits=615360 | BER=9.750e-06 | UB95=9.750e-06
NF= 8.0 dB | totalErr= 1 | totalBits=615360 | BER=1.625e-06 | UB95=1.625e-06
NF= 9.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
NF=10.0 dB | totalErr= 48 | totalBits=615360 | BER=7.800e-05 | UB95=7.800e-05
NF=11.0 dB | totalErr= 3 | totalBits=615360 | BER=4.875e-06 | UB95=4.875e-06
NF=12.0 dB | totalErr= 5 | totalBits=615360 | BER=8.125e-06 | UB95=8.125e-06
NF=13.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
NF=14.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
NF=15.0 dB | totalErr= 0 | totalBits=615360 | BER=0.000e+00 | UB95=4.868e-06
>>
扫一扫
85
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
