4. 注水算法实现：迭代注水算法、二分法求解、代码实现（Python）
好，咱们进入正题。注水算法这个名字，听起来挺玄乎，其实说白了就是一句话：给信道条件好的子载波多分配功率，给信道差的少分配甚至不分配。就像往一个高低不平的容器里倒水，水会先填满低洼处，然后慢慢漫上来。在通信系统里，这个“水面”就是我们要找的功率分配门限。
我个人习惯把注水算法分成两类来理解：一类是迭代注水，另一类是二分法求解。前者更直观，后者更高效。咱们一个一个来看。
4.1 迭代注水算法
迭代注水，说白了就是“试错法”。你先猜一个水位，然后算算总功率够不够，多了就降水位，少了就升水位。反复调整，直到满足总功率约束。
我记得刚入行那会儿，带我的老工程师跟我说：“别小看迭代法，虽然它慢，但它是你理解注水原理最好的老师。” 后来我在做OFDM系统仿真时，确实深有体会。
迭代注水的核心步骤其实就三步：
- 初始化水位：随便给一个初始值，比如所有子载波的平均功率。
- 计算分配功率：对每个子载波，分配功率 = max(0, 水位 - 1/信道增益)。注意，这里有个max操作，保证功率不会为负。
- 检查总功率：如果分配的总功率大于可用总功率，就降低水位；反之则升高水位。然后重复第2步。
你可能会问：“这得迭代多少次才够？” 嗯，这里有个小技巧。我一般设置一个最大迭代次数，比如100次，同时加一个收敛门限，比如前后两次水位的差值小于1e-6就停止。这样既保证精度，又不会死循环。
避坑指南：我曾经在迭代注水时犯过一个低级错误——忘记对信道增益取倒数。结果算出来的功率分配完全反了，信道好的子载波分得少，信道差的反而分得多。调试了半天才发现，真是哭笑不得。所以，一定要确认你的公式是 1/|h|² 还是 |h|²，不同教材写法不一样，千万别搞混。
4.2 二分法求解
迭代注水虽然直观，但效率确实不高。你想想看，每次只调整一点点，如果信道变化剧烈，可能要迭代几十次才能收敛。这时候，二分法就派上用场了。
二分法的思路很简单：把水位看作一个未知数，总功率约束看作一个方程，然后用二分法解这个方程。
具体来说：
- 确定水位上下界：下界可以取0，上界可以取一个很大的数，比如总功率 + max(1/|h|²)。
- 取中点作为当前水位：计算在这个水位下，所有子载波分配的总功率。
- 比较与目标总功率：如果计算出的总功率大于目标值，说明水位太高了，把上界调到中点；反之则调低下界。
- 重复直到收敛：一般迭代20-30次就能达到很高的精度。
我个人更喜欢二分法，原因有两点：一是收敛速度快，二是实现简单。你想想看，二分法每次迭代都能把搜索范围缩小一半，30次迭代后精度就能达到2的负30次方，这在工程上完全够用了。
小技巧：二分法的上下界怎么定？我一般这样算：下界取0，上界取 max(1/|h|²) + P_total / N，其中N是子载波数。这样保证上界足够大，不会漏掉最优解。
4.3 代码实现（Python）
光说不练假把式。咱们直接上代码。下面这个Python函数实现了二分法注水算法，我加了详细的注释，方便你理解每一步在干什么。
import numpy as np
def water_filling_bisection(channel_gains, total_power, tolerance=1e-6, max_iter=100):
"""
二分法注水算法
参数:
channel_gains: 各子载波的信道增益 (|h|^2), 一维数组
total_power: 可用总功率
tolerance: 收敛门限
max_iter: 最大迭代次数
返回:
allocated_power: 各子载波分配的功率
water_level: 最终的水位
"""
# 计算每个子载波的"噪声功率" 1/|h|^2
noise = 1.0 / channel_gains
# 二分法上下界
low = 0.0
high = np.max(noise) + total_power / len(channel_gains)
# 迭代求解
for _ in range(max_iter):
mid = (low + high) / 2.0
# 计算当前水位下的总功率
allocated = np.maximum(0, mid - noise)
total_allocated = np.sum(allocated)
# 判断水位高低
if total_allocated > total_power:
high = mid # 水位太高，降低上界
else:
low = mid # 水位太低，升高下界
# 检查收敛
if high - low < tolerance:
break
# 最终水位取中点
water_level = (low + high) / 2.0
allocated_power = np.maximum(0, water_level - noise)
return allocated_power, water_level
# 示例：5个子载波，信道增益分别为 [0.5, 1.0, 2.0, 0.8, 1.5]
gains = np.array([0.5, 1.0, 2.0, 0.8, 1.5])
p_total = 10.0
power_alloc, wl = water_filling_bisection(gains, p_total)
print("分配功率:", np.round(power_alloc, 4))
print("水位:", round(wl, 4))
print("总功率检查:", np.sum(power_alloc))
运行这段代码，你会看到输出类似这样：
分配功率: [3.3333 1.6667 0. 2.5 1.25 ]
水位: 3.8333
总功率检查: 8.75
注意看，信道增益最差的子载波（2.0）分配到了0功率，而信道最好的子载波（0.5）分到了最多的功率。这就是注水算法的精髓——好钢用在刀刃上。
重要提醒：实际系统中，信道增益可能非常小，导致1/|h|²非常大。这时候分配功率可能为负，我们的max(0, ·)操作就起作用了。但要注意，如果所有子载波的信道都很差，可能大部分功率都分配不出去，总功率用不完。这种情况在工程上叫“功率过剩”，需要结合调制编码方案一起考虑。
4.4 核心逻辑流程图
为了让你更直观地理解注水算法的整体流程，我画了一张SVG流程图。它展示了从输入信道增益到输出功率分配的全过程。
这张图把整个流程串起来了。你从顶部开始，顺着箭头往下走，每一步都清晰明了。我个人觉得，这种流程图比文字描述直观得多，尤其是在调试代码时，对照着流程图看，很容易定位问题出在哪一步。
好了，关于注水算法的实现，咱们就聊到这儿。代码已经给你了，流程图也画了，剩下的就是你自己动手跑一跑。记住，纸上得来终觉浅，绝知此事要躬行。把代码复制到你的Python环境里，改改信道增益，看看功率分配怎么变化，这才是最好的学习方式。