OFDM 自适应功率/比特分配

——Hughes-Hartogs + 闭环功率重排（完全向量化/闭式）

---

1）问题模型

每路子载波（子信道）在接收端等效为：

\(\gamma_i=\frac{|H_i|^2 P_i}{\sigma_n^2},\quad\)
\(\text{SNR}_i=\frac{\gamma_i}{G_i},\quad G_i=\frac{\Gamma}{1-\alpha}\)

• \(P_i\)：分配给子载波 \(i\) 的功率
• \(|H_i|^2\)：信道增益（含路径损耗/阴影/小尺度都可并入）
• \(\sigma_n^2\)：噪声功率（IFFT周期 T，带宽 B：\(\sigma_n^2=N_0 B/N\)常见）
• \(\Gamma\)：SNR gap（QAM 常写作 \(\Gamma=-\frac{\ln(5\mathrm{BER})}{1.5}\)
• \(\alpha\)：编码/滚-off 余量（常见 0.0~0.1）

单子载波每增加 1 bit 所需“最小增量功率”（Hughes-Hartogs 核心）

\(\Delta P_i^{\min}= \frac{ 2^{b_i}-1 }{ \dfrac{|H_i|^2}{\sigma_n^2 G_i} } \quad\Rightarrow\quad \text{下一bit代价}= \boxed{ \frac{1}{\dfrac{|H_i|^2}{\sigma_n^2 G_i}} \;2^{\,b_i} }\)

也就是说：谁“每多传 1 bit 需要的增量功率最小”，就把这 1 bit 给它。

---

2）算法：Hughes-Hartogs（经典自适应比特加载）

2.1 输入/参数约定

N       : 可用子载波数
Ptot    : 总发射功率预算 (W)
N0      : PSD 单边 (W/Hz) 或噪声功率 per bin
H       : channel gains |H_i|^2  (size N×1, ≥0)
targetR : 目标总比特 R=sum b_i (整数)
Gamma   : SNR gap (≥1)
alpha   : margin (0~0.1)

2.2 主函数：alloc_hugheshartogs.m

function [b, P, Ractual, usedP] = alloc_hugheshartogs(H, Ptot, N0, targetR, Gamma, alpha, varargin)
% H      : |H_i|^2 (linear power gain)
% Ptot   : total power budget
% N0     : noise PSD (W/Hz) 或 noise power per tone (看你单位统一)
% targetR: required total bits over N tones
% Gamma  : SNR gap (e.g. 9.8 for BER~1e-7 uncoded approx)
% alpha  : margin (e.g. 0.05)
%
% 若你想“功率自适应到目标R” 且不限制Ptot，可设 Ptot=inf

if nargin>=7, binit=varargin{1}; else binit=zeros(size(H)); end

G = Gamma/(1-alpha);              % 等效gap
Nt = numel(H);
assert(Nt>=2);

% ------ 单位统一：最常见做法 ------
% 假设 noise power per tone = N0 (如果N0是PSD) 则乘 delta_f
% 若你已直接给出 sigma2_i=noise power per tone，就写 sigma2=N0*ones(Nt,1);
delta_f = 1;                       % 设为 1 或一个真实 bin width
sigma2  = N0*delta_f*ones(Nt,1);   % per-tone noise power

% 避免除0
H  = max(H, 1e-12);
sigma2 = max(sigma2, 1e-12);

% -------- 1) 起点（常可设全0或注水起点）--------
b = binit(:);          % integer bits per tone
P = zeros(Nt,1);       % 功率 per tone

% -------- 2) 增量功率代价函数 --------
% cost(i) = 下一bit需要多少功率
% 当前 b(i) 已分配，下1bit → 相当于从 b(i) 到 b(i)+1
cost = @(ii)  ( (2.^(b(ii)+1)-1) ) ./ ( H(ii)./(sigma2(ii)*G) );

% -------- 3) 贪心填充 --------
Rnow = sum(b);
step = 0;
while Rnow < targetR
    % 已分配功率（按当前b）
    for i=1:Nt
        P(i) = ( 2^b(i)-1 ) ./ max(H(i)/(sigma2(i)*G), 1e-12);
    end
    usedP_now = sum(P);

    % 若硬功率约束且已经超过，不能再塞bit（除非你允许过载）
    if isfinite(Ptot) && usedP_now > Ptot
        warning('Power budget exceeded before reaching targetR.');
        break;
    end

    [~, idx] = min( cost(1:Nt) );   % 确定性：谁“增量功率最小”给谁
    b(idx) = b(idx) + 1;
    Rnow = sum(b);

    step = step+1;
    if step>Nt*targetR*2, break; end % 防呆
end

% -------- 4) 按最终b重算功率（确定性闭式）--------
for i=1:Nt
    P(i) = max( (2^b(i)-1) ./ max(H(i)/(sigma2(i)*G), 1e-12) , 0);
end
usedP = sum(P);

% 若你要“把预算Ptot正好用完”，可做等比缩放（不改变bit分配顺序）：
if isfinite(Ptot) && usedP>1e-12
    P = P * (Ptot / usedP);
    usedP = Ptot;
end

Ractual = sum(b);
end

---

3）可选：Chow截断注水（802.20 style）

这是所谓 Water-filling with discrete bits：先算注水水位，再截断到整数bit与功率上限。仍然只用到排序/求根，不靠fmincon等“搜索工具”。

3.1 关键公式

注水：

\(P_i^*=\Big(\nu-\frac{\sigma_n^2}{|H_i|^2}\Big)^+, \quad b_i=\log_2\Big(1+\frac{P_i^*|H_i|^2}{\sigma_n^2}\Big)/\log_2(1+1/\Gamma)\)

工程上用 截断位：

\(b_i=\max\Big(0,\;\Big\lfloor \frac{1}{2} \log_2\Big( \frac{|H_i|^2} {\nu\,\sigma_n^2\,\Gamma} \Big) \Big\rfloor \Big)\)

并且功率回到：

\(P_i=\frac{(2^{b_i}-1)\,\sigma_n^2\,\Gamma}{|H_i|^2}\)

3.2 闭式求水位 (\nu)

水位 (\nu) 由功率约束决定：

\(\sum_{i\in\mathcal{S}(\nu)}\Big(\nu-\frac{\sigma_n^2}{|H_i|^2}\Big)=P_{\text{tot}}\)

对 (\nu) 单调增 ⇒ 用 二分法：

function [b,P,nu] = alloc_chow_waterfilling(H, Ptot, N0, Gamma, varargin)
% 这里用 per-tone noise = N0
Nt=numel(H); H=max(H(:),1e-12); sigma2=N0*ones(Nt,1);
if nargin>=5, margin=varargin{1}; else margin=0; end
G = Gamma/(1-margin);

% 排序辅助
invH = sigma2./H;   % = sigma2/|H|^2

% 可行集随nu扩大：S(nu)=find(nu>invH)
% 二分 nu ∈ [min(invH)+eps, Ptot/Nt + max(invH)]
nu_l = max(invH)+1e-6;
nu_r = Ptot + max(invH) + 1;
for it=1:60
    nu_m = (nu_l+nu_r)/2;
    Psum = sum(max(nu_m - invH,0));
    if Psum > Ptot, nu_l=nu_m; else nu_r=nu_m; end
    if abs(Psum-Ptot)/Ptot<1e-5, break; end
end
nu = (nu_l+nu_r)/2;

active = nu > invH;
P = zeros(Nt,1); b = zeros(Nt,1);
for i=1:Nt
    if active(i)
        P(i) = nu - invH(i);
        b(i) = floor( 0.5*log2( max(nu*H(i)/(sigma2(i)*G),1) ) );
        b(i) = max(b(i),0);
        % 转回确定功率
        P(i) = (2^b(i)-1)*sigma2(i)*G / H(i);
    end
end
nu=nu;
end

参考代码 OFDM系统的自适应功率和比特分配算法 www.youwenfan.com/contentcnv/79421.html

4）最小可跑 Demo

%% demo：频率选择性衰落 + 自适应 loading
N=32; Ptot=1; N0=1e-3; Gamma=9.8; targetR=120;

% 简单指数衰减速衰落信道
H = rand(N,1).*(0.8.^((0:N-1)')./max(0.8.^((0:N-1)'),1));
H = H + 1e-3; % avoid dead carriers

[b,P,Ract,usedP]=alloc_hugheshartogs(H,Ptot,N0,targetR,Gamma,0.05);

figure;
subplot(3,1,1); stem(1:N,b,'filled'); ylabel('b_i (bits)'); title('Bit Loading');
subplot(3,1,2); stem(1:N,P,'filled'); ylabel('P_i'); title('Power Allocation');
subplot(3,1,3); semilogy(1:N,H,'o-'); ylabel('|H_i|^2'); title('Channel gains');
sgtitle(sprintf('R=%d/%d bits, Pused=%.4f/%.2f',Ract,targetR,usedP,Ptot));