QPSK / 8QAM（8PSK-like 8-QAM）— SNR(Eb/N0)–BER 曲线（仿真 + 理论对照）

说明：这里的 8QAM 采用 方形 8QAM（2×4 星座），功率归一化，Gray 比特映射；QPSK 就是 4-QAM（功率归一化，相邻点距离一致）。
若你说的 8QAM 实际指 8PSK，最后我也给你替换成 8PSK（恒包络）版本。

---

1）核心：Gray 映射函数（可复用）

%% ========== 调制/解调 工具 ==========
function [sym, bits] = qpsk_mod(b)
% b: Nbits×1 (0/1), Nbits 偶数
% 输出 sym: Nsym×1 complex, 单位平均功率=1
    M = 4; Nsym = length(b)/2; assert(mod(length(b),2)==0);
    mapping = [1+1i, -1+1i, -1-1i, 1-1i]/sqrt(2); % Gray: 00,01,11,10
    idx = bi2de(reshape(b,2,[])','left-msb')+1;
    sym = mapping(idx).';
end

function [sym, bits] = qam8_mod(b)
% 8QAM 方形 (±1±j, ±3±j 去掉四角)
% 2×4 点集: 平均功率归一化
    pts = [ ...
      -1+1j,   1+1j, ...
      -1-1j,   1-1j, ...
      -3+1j,   3+1j, ...
      -3-1j,   3-1j ];
    % Gray 按 I,Q 各自 Gray 编码 (3bit/sym: b0b1|b2)
    % b0b1->I∈{-3,-1,1,3}  Gray:00:-3,01:-1,11:1,10:3  (顺序可调)
    Ilev = [-3,-1,1,3];  Jlev = [-1,1];
    % bit->idx: b0b1->Ii(0..3), b2->Jj(0/1)
    nb = length(b); Nsym = nb/3; assert(mod(nb,3)==0);
    bmat = reshape(b,3,[])';
    ii = bi2de(bmat(:,[1 2]),'left-msb')+1; % 0..3
    jj = bmat(:,3)+1;                        % 0..1
    I = Ilev(ii.').';
    J = Jlev(jj.').';
    sym_raw = complex(I,J);
    % 功率归一化
    Es = mean(abs(sym_raw).^2);
    sym = sym_raw / sqrt(Es);   % 现在 avg Es=1
end

function b_rec = qpsk_demod(sym)
    % 硬判决 -> 2bit/sym
    s = sym(:);
    b1 = (real(s)<0);                 % MSB
    b0 = (imag(s)<0);                 % LSB（取决于你Gray）
    b_rec = reshape([b1(:).'; b0(:).'], [],1);
end

function b_rec = qam8_demod(sym)
    s = sym(:);
    % 逐符号最小距离
    Ilev =[-3,-1,1,3]; Jlev=[-1,1];
    Es = (mean(abs(complex(Ilev(2),Jlev(2))).^2)*2 ...
          + mean(abs(complex(Ilev(3),Jlev(2))).^2)*2)/4; % 与发端一致
    pts = [complex(Ilev(1),Jlev(1)),complex(Ilev(2),Jlev(1)), ...
           complex(Ilev(1),Jlev(2)),complex(Ilev(2),Jlev(2)), ...
           complex(Ilev(3),Jlev(1)),complex(Ilev(4),Jlev(1)), ...
           complex(Ilev(3),Jlev(2)),complex(Ilev(4),Jlev(2))];
    pts = pts / sqrt(mean(abs(pts).^2));
    b_rec = zeros(3*length(s),1);
    for k=1:length(s)
        [~,idx]=min(abs(s(k)-pts));
        % 逆映射：idx->bits (按发端固定顺序)
        % 这里用查表最稳：
        lut_bits = [0 0 0;
                    1 0 0;
                    0 1 0;
                    1 1 0;
                    0 0 1;
                    1 0 1;
                    0 1 1;
                    1 1 1]; % 你得按你的Gray规则严格对齐
        b_rec((k-1)*3+(1:3)) = lut_bits(idx,:).';
    end
end

你若嫌手写 8QAM 查表麻烦，最稳是用 qammod(...,'UnitAveragePower',true,'InputType','bit')（需要 Comm Toolbox）。

---

2）Monte Carlo 主循环（Eb/N0 → BER）

下面代码 不依赖工具箱，只靠 AWGN 生成，EbN0 从 0 到 12 dB。

%% ========== 主仿真 ==========
clear; clc;

Nbits_per_run = 2e5;     % 每EbN0 发送比特数（可调大换精度）
EbN0_dB = 0:2:12;
Nsnr = length(EbN0_dB);

BER_QPSK = nan(1,Nsnr);
BER_8QAM = nan(1,Nsnr);

rng(42);

for s = 1:Nsnr
    EbN0 = 10^(EbN0_dB(s)/10);

    %% ===== QPSK =====
    k = 2; M = 4;
    % QPSK 归一化：Es=1，故 Eb=Es/k=0.5 → N0 = Eb/EbN0 =0.5/EbN0
    Eb = 1/k; N0 = Eb / EbN0;
    sigma = sqrt(N0/2);   % AWGN 方差 per dim

    b = randi([0 1], Nbits_per_run,1);
    sym = qpsk_mod(b);    % Es=1

    noise = sigma*(randn(size(sym)) + 1i*randn(size(sym)));
    r = sym + noise;

    b_hat = qpsk_demod(r);
    BER_QPSK(s) = sum(b~=b_hat)/length(b);

    %% ===== 8QAM (方形 8QAM) =====
    k = 3; M = 8;
    Eb = 1/k; N0 = Eb / EbN0;
    sigma = sqrt(N0/2);

    b = randi([0 1], Nbits_per_run,1);
    sym = qam8_mod(b);    % 发端已经归一化到 Es=1

    noise = sigma*(randn(size(sym)) + 1i*randn(size(sym)));
    r = sym + noise;

    b_hat = qam8_demod(r);
    BER_8QAM(s) = sum(b~=b_hat)/length(b);

    fprintf('Eb/N0=%ddB  QPSK BER=%.3e  8QAM BER=%.3e\n', ...
        EbN0_dB(s), BER_QPSK(s), BER_8QAM(s));
end

---

3）理论近似曲线

% 理论：QPSK = BPSK 等价（Gray）
BER_QPSK_th = 0.5*erfc( sqrt(10.^(EbN0_dB/10)) );

% 8QAM（方形）近似：可拆成两路PAM（4电平+2电平）再Gray近似
% 常用工程近似（AWGN,Gray）：
gamma = 10.^(EbN0_dB/10);  % Eb/N0 linear
% 方形 8QAM：Es=1 => d_min^2 ≈ (2/Es_norm_factor)，但用“等效 2bit路”近似更直观：
% 这里给经典上界/近似：
Pb_4PAM = (2/4)*0.5*erfc( sqrt(gamma * (3/(M^2-1)) ) ); % M=4 PAM
Pb_2PAM = 0.5 * erfc( sqrt(gamma) );              % 2PAM（I路0/1简化）
BER_8QAM_app = (2/3)*Pb_4PAM + (1/3)*Pb_2PAM;      % 加权（保守可用 Pb_4PAM直接）

参考代码 QPSK编码，8QAM编码下的多种求解信噪比和误码率的曲线 www.youwenfan.com/contentcnv/81156.html

4）画图

figure('Color','w');
semilogy(EbN0_dB, BER_QPSK, 'o-','LineWidth',2,'DisplayName','QPSK 仿真');
hold on; semilogy(EbN0_dB, BER_QPSK_th,'--','DisplayName','QPSK 理论');
semilogy(EbN0_dB, BER_8QAM, 's-','LineWidth',2,'DisplayName','方形8QAM 仿真');
semilogy(EbN0_dB, BER_8QAM_app,':','DisplayName','8QAM 近似上界');

grid on; xlabel('E_b/N_0 (dB)'); ylabel('BER');
title('QPSK & 方形8QAM — Eb/N0 对 BER 曲线');
legend('Location','southwest');
ylim([1e-5 1])