简介适用于认知无线电功率分配研究的注水算法MATLAB实现包面向无线通信、信号处理方向的初学者与科研人员解决多信道环境下总发射功率受限时的速率优化问题。压缩包共3个文件包含2个可直接运行的M源码文件和1个说明性文本整体大小仅2KB代码精炼适合快速理解算法核心流程。已有1226人学习下载可见其作为入门参考的实用价值。文件中提供了注水算法的完整实现覆盖信道参数读取、水位阈值设定与迭代分配等关键逻辑可帮助读者将理论公式转化为可运行程序进一步分析不同信道条件下的功率分配结果并为认知无线电场景中的主用户干扰约束提供基础示例。配套文本或为相关链接与说明便于获取更多背景资料。1. 注水算法为什么是认知无线电功率分配的首选模型在多信道无线通信里频谱资源不是均匀可用的有些子信道噪声低、增益高有些则被主用户占用或受到强干扰。认知无线电的次用户要在不干扰主用户的前提下最大化自己的传输速率本质上就是一个带约束的凸优化问题。注水算法Water-Filling从香农信道容量公式出发把总功率像水一样注入“信道容量凹槽”噪声低、增益高的信道先分到功率直到水位线统一。这个结论不是经验近似而是拉格朗日乘子法直接导出的闭式解所以它在理论上是严密的在工程上又比线性规划简单得多。本文用 MATLAB 源码zhushuixian.m和powerallo.m拆解两种常见实现并给出认知无线电干扰温度约束下的完整可运行代码适合正在做无线资源管理课题、需要快速验证算法性能的研究生和工程师。2. 注水问题的数学建模与 KKT 条件推导2.1 多信道下的 Shannon 容量与目标函数考虑一个认知无线电次用户系统的下行链路OFDM 子载波或感知得到的空闲频段被划分为 (N) 个并行信道。第 (i) 个信道的信道增益为 (h_i)噪声功率谱密度为 (N_{0,i})发射功率为 (p_i)。在总功率约束下最大化总容量的标准形式为[ \max_{p_i} \sum_{i1}^{N} B \log_2\left(1 \frac{|h_i|^2 p_i}{N_{0,i}B}\right) ][ \text{s.t.} \quad \sum_{i1}^{N} p_i \le P_{total}, \quad p_i \ge 0 ]其中 (B) 是子信道带宽。令信道噪声等效系数 (n_i N_{0,i}B / |h_i|^2)则目标函数变成 (\sum \log_2(1 p_i / n_i))。这里 (n_i) 的物理含义是“让该信道达到单位信噪比所需的最小发射功率”它越小代表信道质量越好。认知无线电场景下如果第 (i) 个信道靠近主用户接收机还需要额外叠加一个干扰权重但最基础的注水问题先不考虑这一层后面第 4 章再展开。用拉格朗日乘子法处理约束。构造[ L \sum_{i1}^{N} \log_2\left(1 \frac{p_i}{n_i}\right) - \lambda \left(\sum_{i1}^{N} p_i - P_{total}\right) ]对 (p_i) 求偏导并令其为零得到条件 (p_i n_i 1/(\lambda \ln 2))。右侧是个常数记为 (\mu)因此最优功率满足 (p_i \mu - n_i)。这就是“注水”的由来(\mu) 是统一的水位线所有信道的 (p_i n_i) 都被抬到同一高度。如果某个信道的噪声系数 (n_i) 高于水位线则 (p_i) 为负不符合约束此时该信道应当关闭即 (p_i 0)。2.2 拉格朗日乘子与“水位线”的物理含义水位线 (\mu) 不是随便选的它由总功率约束唯一确定。闭式表达为[ \mu \frac{P_{total} \sum_{i \in \mathcal{A}} n_i}{|\mathcal{A}|} ]其中 (\mathcal{A}) 是参与功率分配的信道集合即满足 (\mu n_i) 的信道。这个公式在 MATLAB 里非常有用因为你可以直接根据当前集合计算 (\mu)而不需要像梯度下降那样迭代很多轮。2.2.1 为什么它是最优解从 KKT 条件看该问题满足 Slater 条件因此强对偶成立。注水解满足三个条件主可行性功率非负且总和不超过上限、对偶可行性(\lambda \ge 0)、互补松弛(\lambda(\sum p_i - P_{total}) 0)。只要按上述公式求解这三个条件自然成立。实际实现里常见的错误是忘记剔除 (p_i 0) 的信道导致水位线偏高、总功率溢出。正确做法是循环剔除先按所有信道计算去掉不合格信道后重新算直到所有 (p_i) 非负。3. MATLAB 实现注水算法从单循环到向量化3.1 核心函数 zhushuixian.m 逐行拆解压缩包里的zhushuixian.m是经典实现我见过很多教材里的版本都长这样。它接受信道噪声系数向量和总功率返回分配功率向量。function p zhushuixian(n, P_total) % n: 1*N 向量每个信道的噪声系数已归一化 % P_total: 标量总可用功率 % p: 1*N 向量分配功率 N length(n); p zeros(1, N); idx 1:N; % 初始认为所有信道都参与 mu (P_total sum(n)) / N; % 初始水位线 while true p_sel mu - n(idx); % 候选功率 if all(p_sel 0) break; end idx idx(p_sel 0); % 剔除不合格信道 mu (P_total sum(n(idx))) / length(idx); % 重算水位线 end p(idx) mu - n(idx); end逻辑说明idx保存当前参与分配的信道下标第一次假设所有信道都分到功率如果某信道噪声系数太大导致p_sel为负就把它从参与集合中剔除。剔除后总和公式里的 (n_i) 项减少分母也减少水位线变化因此需要循环重算。这个做法和理论推导完全一致复杂度最坏是 (O(N^2))但实际循环次数很少。参数说明n必须预先除以信道增益平方即n noise_power ./ (abs(h).^2)P_total的单位必须和n一致否则结果尺度全错。3.2 经典迭代分配与排序注入法对比powerallo.m在工程上更常用它先对信道质量排序然后从最好的信道开始逐个“注水”。这种方法的好处是可以提前知道哪些信道被激活也便于和等功率分配做比较。function p powerallo(n, P_total) % 排序注入法先按噪声系数升序排序再计算水位线 [n_sorted, idx] sort(n); N length(n); p_sorted zeros(1, N); % 从最好的信道开始累加找到最后一个被激活的信道 active 0; for k 1:N mu (P_total sum(n_sorted(1:k))) / k; if k N mu n_sorted(k1) active k; break; % 当前水位线无法覆盖下一个信道停止 end active k; end % 按最终水位线分配 mu (P_total sum(n_sorted(1:active))) / active; p_sorted(1:active) mu - n_sorted(1:active); p zeros(1, N); p(idx) p_sorted; % 映射回原始信道顺序 end逻辑说明sort按噪声系数从小到大排n_sorted(1)是质量最好的信道。循环里计算“前 k 个信道的公共水位线”如果这个水位线仍然小于第 k1 个信道的噪声系数说明第 k1 个信道分不到功率循环终止。这种方法的优势是只需一次排序和一次扫描复杂度主要是排序的 (O(N \log N))。参数说明idx用来记录排序前后的映射关系分配结果必须还原到原始信道序号否则后续叠加干扰约束时会张冠李戴。3.3 向量化实现与运行效率对比当信道数达到上千个比如 OFDMA 系统的子载波级别上面的循环版本依然很快。但如果你在蒙特卡洛仿真里跑几万次向量化版本能把时间再压一个量级。我一般用批量矩阵运算一次处理多组随机信道实现而不是逐次循环。% 批量模式下 n_matrix 是 M*N每行是一组信道噪声系数 % P_total_vec 是 M*1每行对应一组总功率 M size(n_matrix, 1); N size(n_matrix, 2); mu_all (P_total_vec sum(n_matrix, 2)) ./ N; % M*1 初始水位线 p_matrix mu_all - n_matrix; % M*N 初始功率 bad p_matrix 0; % 剔除不合格信道后重算这里用逻辑索引矩阵但需循环直到稳定 while any(bad(:)) valid ~bad; active_count sum(valid, 2); mu_all (P_total_vec sum(n_matrix .* valid, 2)) ./ active_count; p_matrix mu_all - n_matrix; new_bad p_matrix 0; if isequal(new_bad, bad) break; end bad new_bad; end p_matrix(~bad) p_matrix(~bad); p_matrix(bad) 0;逻辑说明valid是布尔矩阵sum(n_matrix .* valid, 2)只累加未被剔除的信道。每一次循环把不合格位置置零并重算水位线直到前后两次的剔除集合不再变化。注意 MATLAB 的布尔索引比用find快尤其是在大矩阵下。参数说明P_total_vec必须是列向量n_matrix每行一组这样mu_all的广播才能正确匹配。运行效率对比结果通常如此循环实现约 0.3 ms排序注入法约 0.4 msN1024向量化批量处理一万组时平均每组时间降到 0.02 ms。4. 认知无线电场景下的干扰温度约束与 powerallo.m 实现4.1 主用户干扰限制怎么进入目标函数认知无线电和普通多信道系统最大的区别在于次用户发射功率不仅受自身总功率限制还要保证在主用户接收机处产生的干扰功率低于干扰温度。假设主用户接收机在第 (i) 个信道上感受到的干扰系数为 (g_i)干扰门限为 (Q)那么额外约束变为[ \sum_{i1}^{N} g_i p_i \le Q ]此时问题变成两个线性约束的凸优化。如果 (g_i) 对所有信道都一样那么干扰约束等价于把总功率上限改成 (\min(P_{total}, Q / g))直接套用第 3 章的代码即可。但实际中 (g_i) 和信道增益 (h_i) 没有固定比例关系因为主用户和次用户的空间位置不同路径损耗、阴影衰落都不一样。此时需要把干扰约束也拉进拉格朗日函数引入第二个乘子 (\nu)最优解变成分段形式[ p_i \max\left(0, \frac{1}{\lambda \nu g_i} - n_i\right) ]这是一个广义注水水位线不再是单一常数而是随 (g_i) 变化的“倾斜水位线”。4.2 带干扰线约束的功率分配代码powerallo.m的完整版应该处理这种广义注水。常见做法是外循环找 (\lambda)内循环找 (\nu)因为 (\lambda) 和 (\nu) 相互影响。这里的实现采用二分搜索嵌套。function p powerallo_cr(n, g, P_total, Q) % n: 1*N 噪声系数 % g: 1*N 主用户干扰系数 % P_total: 总功率上限 % Q: 干扰温度上限 lambda_low 0; lambda_high 1e10; for iteration 1:200 % 外层二分lambda lambda (lambda_low lambda_high) / 2; % 内层找nu满足干扰约束 nu_low 0; nu_high 1e10; for inner 1:200 nu (nu_low nu_high) / 2; p max(0, 1 ./ (lambda nu * g) - n); interference sum(g .* p); if interference Q nu_high nu; else nu_low nu; end end % 检查总功率约束 used_power sum(p); if used_power P_total lambda_high lambda; % 用不完功率说明lambda太大 else lambda_low lambda; end if abs(used_power - P_total) 1e-9 break; end end end逻辑说明内层循环固定 (\lambda)用二分法调整 (\nu)让干扰量收敛到 (Q) 附近外层循环再调整 (\lambda) 让总功率收敛到 (P_{total})。二分边界设置很关键功率和干扰都有物理上限(1e10) 是安全值。如果某个信道1 ./ (lambda nu * g)算出来是 NaN多半是lambda nu * g出现了 0 或负数需要检查输入 (g) 是否为严格的正常数。4.3 参数调节与失败模式实际执行时最常遇到的三种失败模式如下。第一Q设置得太小导致干扰约束比总功率约束更紧外层二分搜不到可行解。此时可以先用无约束注水算一遍看干扰量sum(g .* p_no_constraint)是否超过Q如果超过说明这个场景下次用户需要降低总功率来避开主用户P_total实际有效值要缩小。第二g里有零元素导致干扰约束对某些信道无约束力公式里出现除以零风险常见做法是给g加一个极小量eps避免数值病态。第三外层二分循环次数不足当 (N) 很大时总功率随 (\lambda) 的变化在高维下可能很陡200 次二分通常够但如果你把P_total设置得非常大lambda会非常小此时要改用对数刻度搜索即lambda 10^(log10(lambda_low) ...)否则低水位线区域精度不足。运行完成后可以用下面的代码验证 KKT 条件是否满足% 验证总功率和干扰约束 assert(abs(sum(p) - min(P_total, sum(p))) 1e-6); assert(sum(g .* p) Q 1e-6); % 验证拉格朗日条件激活信道满足 1/(lambdanu*g) p n active p 1e-10; residual max(abs(1 ./ (lambda nu * g(active)) - p(active) - n(active)));5. 验证技巧收敛性检查与网格搜索基准5.1 用穷举法验证注水结果对于 N 较小时比如 2 到 4 个信道最稳妥的验证方式是用 MATLAB优化工具箱的fmincon做对照fun (p) -sum(log2(1 p ./ n)); x0 ones(1, N) * P_total / N; Aeq ones(1, N); beq P_total; lb zeros(1, N); ub []; opt optimoptions(fmincon, Display, off, Algorithm, sqp); p_fmincon fmincon(fun, x0, [], [], Aeq, beq, lb, ub, [], opt);这里把总功率约束写成等式因为注水解一定用完全部功率除非干扰约束先碰到。比较p_fmincon和zhushuixian的输出差值应小于 (10^{-6})。如果对不上优先检查噪声系数是否漏了abs(h).^2归一化。5.2 常见坑信道为零、水位线为负、浮点误差第一信道增益h为零时n变成无穷大排序后该信道永远排在最后不会参与分配但代码里sum(n(idx))可能计算出 Inf导致水位线 NaN。解决办法是先把n中大于某个阈值例如1e10 * max(n(n Inf))的元素剔除。第二P_total远小于所有n时注水解会只给最好的信道分配极少量功率而此时水位线接近min(n)p中会出现数量级为 (10^{-15}) 的负值用max(0, ...)截断即可但不该截断后参与求和否则总功率偏低。第三浮点误差会让p的总和差1e-12在后续叠加干扰约束时可能被判违规建议在分配完成后做一次缩放归一化p p * (P_total / sum(p))只对激活信道缩放保持零功率信道不变。最后给一个实用技巧在循环仿真对比算法性能时不需要每次调用二分函数。把注水结果缓存下来只有当信道实现变化时才重算如果信道变化不大可以用上一次的水位线作为初始值二分搜索的第一步就直接到达收敛区通常一次迭代就够了。对于认知无线电研究来说真正耗时的是上千次蒙特卡洛下的主用户位置和信道快照更新注水算法本身只占不到 1% 的运行时间。本文还有配套的精品资源点击获取