注水法（Water-filling）解决MIMO信道功率优化

注水法（Water-filling）解决MIMO信道功率优化，通过智能地将更多功率分配给信道条件更好的子信道，实现信道容量的最大化。

一、注水法核心原理

1.1 物理直观理解

想象一系列高度不同的容器（代表各子信道的逆增益1/λᵢ），向其中注水（代表总功率P）。水位（代表注水电平μ）会均匀上升，每个容器中的水深（代表分配给该子信道的功率pᵢ）等于水位减去容器高度，且水深≥0。

graph TD
    A[总发射功率约束 P_total] --> B(计算信道增益 λ_i);
    B --> C{排序增益并求逆 h_i = 1/λ_i};
    C --> D[初始化注水电平 μ];
    D --> E[计算各信道功率 p_i = μ - h_i];
    E --> F{检查 p_i ≥ 0 且<br>Σp_i = P_total?};
    F -- 否 --> G[调整 μ 重新计算];
    F -- 是 --> H[获得最优功率分配];
    H --> I[计算最大信道容量];

1.2 优化问题数学模型

对于具有N个子信道的MIMO系统，优化问题为：

最大化：C = Σ log₂(1 + pᵢ·λᵢ/σ²)
约束条件：Σ pᵢ ≤ P_total，pᵢ ≥ 0

其中：

• λᵢ：第i个子信道的增益（信道矩阵奇异值的平方）
• σ²：噪声功率
• pᵢ：分配给第i个子信道的功率

二、完整MATLAB实现

2.1 基础注水算法

function [p_opt, mu, capacity] = waterfilling_basic(lambda, P_total, sigma2)
% 基础注水法功率分配
% 输入：
%   lambda - 信道增益向量 (N×1)，λ_i = svd(H)^2
%   P_total - 总发射功率
%   sigma2 - 噪声功率
% 输出：
%   p_opt - 最优功率分配向量
%   mu - 注水电平
%   capacity - 可达信道容量

N = length(lambda);
h = 1 ./ lambda;  % 容器高度（逆增益）

% 对高度进行排序
[h_sorted, idx] = sort(h, 'ascend');

% 迭代搜索注水电平μ
mu_low = 0;
mu_high = P_total + max(h_sorted);  % μ的上界

for iter = 1:1000  % 二分搜索
    mu = (mu_low + mu_high) / 2;
    
    % 计算当前μ下的功率分配
    p_temp = max(0, mu - h_sorted);
    
    % 计算总功率
    total_power = sum(p_temp);
    
    % 调整搜索边界
    if total_power < P_total - 1e-10
        mu_low = mu;
    elseif total_power > P_total + 1e-10
        mu_high = mu;
    else
        break;
    end
end

% 恢复原始顺序
p_opt_sorted = p_temp;
p_opt = zeros(N, 1);
p_opt(idx) = p_opt_sorted;

% 计算信道容量
capacity = sum(log2(1 + p_opt .* lambda / sigma2));

fprintf('注水算法结果：\n');
fprintf('注水电平 μ = %.4f\n', mu);
fprintf('总功率约束 = %.4f, 实际分配 = %.4f\n', P_total, sum(p_opt));
fprintf('信道容量 = %.4f bits/s/Hz\n\n', capacity);
end

2.2 更高效的向量化实现

function [p_opt, mu] = waterfilling_vectorized(lambda, P_total, sigma2)
% 向量化注水法实现（更高效）
% 基于闭式解：μ = (P_total + Σ_{i=1}^k h_i) / k

N = length(lambda);
h = 1 ./ lambda;
[h_sorted, idx] = sort(h, 'ascend');

% 累积高度和
cumsum_h = cumsum(h_sorted);

% 寻找最优的k（活跃信道数）
k_opt = N;
for k = N:-1:1
    mu_temp = (P_total + cumsum_h(k)) / k;
    
    % 检查第k个信道是否活跃
    if mu_temp > h_sorted(k)
        k_opt = k;
        mu = mu_temp;
        break;
    end
end

% 计算功率分配
p_sorted = max(0, mu - h_sorted(1:k_opt));
p_sorted = [p_sorted; zeros(N-k_opt, 1)];

% 恢复原始顺序
p_opt = zeros(N, 1);
p_opt(idx) = p_sorted;

% 验证结果
assert(abs(sum(p_opt) - P_total) < 1e-10, '功率分配不满足约束');
end

2.3 完整MIMO系统仿真

%% MIMO系统注水法功率分配完整仿真
clear; close all; clc;

% ========== 1. 系统参数设置 ==========
Nt = 4;                 % 发射天线数
Nr = 4;                 % 接收天线数
P_total_dB = 10;        % 总功率 (dB)
P_total = 10^(P_total_dB/10); % 线性功率
sigma2_dB = 0;          % 噪声功率 (dB)
sigma2 = 10^(sigma2_dB/10);

num_channels = 1000;    % 仿真信道数量

% ========== 2. 性能统计初始化 ==========
capacity_waterfill = zeros(num_channels, 1);
capacity_uniform = zeros(num_channels, 1);
power_allocation_history = zeros(Nt, num_channels);

% ========== 3. 主仿真循环 ==========
fprintf('开始MIMO注水法功率分配仿真...\n');
fprintf('系统配置: %dx%d MIMO, 总功率: %.1f dB, 噪声功率: %.1f dB\n\n', ...
        Nt, Nr, P_total_dB, sigma2_dB);

for ch_idx = 1:num_channels
    % 生成随机瑞利衰落信道
    H = (randn(Nr, Nt) + 1i*randn(Nr, Nt)) / sqrt(2);
    
    % 计算信道增益（奇异值分解）
    [~, S, ~] = svd(H);
    lambda = diag(S).^2;  % 信道增益λ_i = σ_i^2
    
    % 方法1：注水法功率分配
    [p_wf, mu, cap_wf] = waterfilling_basic(lambda, P_total, sigma2);
    capacity_waterfill(ch_idx) = cap_wf;
    
    % 方法2：均匀功率分配（作为对比）
    p_uniform = P_total / Nt * ones(Nt, 1);
    cap_uniform = sum(log2(1 + p_uniform .* lambda / sigma2));
    capacity_uniform(ch_idx) = cap_uniform;
    
    % 记录功率分配
    power_allocation_history(:, ch_idx) = p_wf;
    
    % 显示前几个信道的详细情况
    if ch_idx <= 3
        fprintf('=== 信道实例 %d ===\n', ch_idx);
        fprintf('信道增益 λ_i: %s\n', num2str(lambda', '%.3f  '));
        fprintf('注水法功率分配: %s\n', num2str(p_wf', '%.3f  '));
        fprintf('均匀分配: %s\n', num2str(p_uniform', '%.3f  '));
        fprintf('容量增益: %.2f%%\n\n', ...
                100*(cap_wf/cap_uniform - 1));
    end
end

% ========== 4. 结果分析与可视化 ==========
figure('Position', [100, 100, 1400, 800]);

% 子图1：容量累积分布函数
subplot(2, 3, 1);
[cdf_wf, x_wf] = ecdf(capacity_waterfill);
[cdf_uni, x_uni] = ecdf(capacity_uniform);
plot(x_wf, cdf_wf, 'b-', 'LineWidth', 2); hold on;
plot(x_uni, cdf_uni, 'r--', 'LineWidth', 2);
xlabel('信道容量 (bits/s/Hz)'); ylabel('累积概率');
title('信道容量CDF比较');
legend('注水法', '均匀分配', 'Location', 'southeast');
grid on;

% 子图2：容量增益直方图
subplot(2, 3, 2);
gain_percentage = 100 * (capacity_waterfill ./ capacity_uniform - 1);
histogram(gain_percentage, 30, 'FaceColor', [0.2, 0.6, 0.8]);
xlabel('容量增益 (%)'); ylabel('频数');
title(sprintf('容量增益分布\n平均增益: %.2f%%', mean(gain_percentage)));
grid on;

% 子图3：功率分配模式示例
subplot(2, 3, 3);
sample_idx = 1;
lambda_sample = lambda;
p_wf_sample = p_wf;

% 创建注水法可视化
bar_data = [lambda_sample/max(lambda_sample), p_wf_sample/max(p_wf_sample)]';
bar(1:Nt, bar_data', 'grouped');
xlabel('子信道索引'); ylabel('归一化值');
title('注水法功率分配示例');
legend('信道增益(归一化)', '分配功率(归一化)', 'Location', 'northwest');
grid on;

% 子图4：不同信噪比下的性能
subplot(2, 3, 4);
SNR_range = -10:2:20;  % dB
capacity_vs_snr = zeros(length(SNR_range), 2);

for snr_idx = 1:length(SNR_range)
    P_total_current = 10^(SNR_range(snr_idx)/10);
    
    % 使用平均统计计算
    cap_wf_sum = 0;
    cap_uni_sum = 0;
    
    % 简化计算：使用代表性信道
    for rep = 1:100
        H_rep = (randn(Nr, Nt) + 1i*randn(Nr, Nt)) / sqrt(2);
        [~, S_rep, ~] = svd(H_rep);
        lambda_rep = diag(S_rep).^2;
        
        [p_wf_rep, ~, cap_wf_rep] = waterfilling_basic(lambda_rep, P_total_current, sigma2);
        p_uni_rep = P_total_current / Nt * ones(Nt, 1);
        cap_uni_rep = sum(log2(1 + p_uni_rep .* lambda_rep / sigma2));
        
        cap_wf_sum = cap_wf_sum + cap_wf_rep;
        cap_uni_sum = cap_uni_sum + cap_uni_rep;
    end
    
    capacity_vs_snr(snr_idx, 1) = cap_wf_sum / 100;
    capacity_vs_snr(snr_idx, 2) = cap_uni_sum / 100;
end

plot(SNR_range, capacity_vs_snr(:, 1), 'b-o', 'LineWidth', 2, 'MarkerSize', 6);
hold on;
plot(SNR_range, capacity_vs_snr(:, 2), 'r--s', 'LineWidth', 2, 'MarkerSize', 6);
xlabel('SNR (dB)'); ylabel('平均信道容量 (bits/s/Hz)');
title('不同SNR下的容量比较');
legend('注水法', '均匀分配', 'Location', 'northwest');
grid on;

% 子图5：功率分配矩阵（热图）
subplot(2, 3, 5);
% 选择信道增益差异最大的几个信道展示
lambda_variation = std(power_allocation_history, 0, 1);
[~, top_indices] = sort(lambda_variation, 'descend');
top_channels = top_indices(1:min(10, num_channels));

power_matrix = power_allocation_history(:, top_indices(1:10));
imagesc(power_matrix);
colorbar; colormap('jet');
xlabel('信道实例'); ylabel('发射天线');
title('功率分配矩阵（前10个变化最大的信道）');
set(gca, 'XTick', 1:length(top_indices(1:10)));

% 子图6：注水法可视化示意图
subplot(2, 3, 6);
% 创建一个简单的注水法示意图
h_example = sort(1./lambda_sample, 'ascend');
mu_example = mu;

% 绘制"容器"
for i = 1:length(h_example)
    rectangle('Position', [i-0.4, 0, 0.8, h_example(i)], ...
              'FaceColor', [0.9, 0.9, 0.9], 'EdgeColor', 'k');
end

% 绘制"水位线"
line([0.5, length(h_example)+0.5], [mu_example, mu_example], ...
     'Color', 'b', 'LineWidth', 3, 'LineStyle', '-');

% 填充"水"
for i = 1:length(h_example)
    if mu_example > h_example(i)
        rectangle('Position', [i-0.4, h_example(i), 0.8, mu_example-h_example(i)], ...
                  'FaceColor', [0.2, 0.4, 0.8], 'EdgeColor', 'none', 'FaceAlpha', 0.6);
    end
end

xlabel('子信道（按增益排序）'); ylabel('高度 / 功率');
title(sprintf('注水法示意图\n水位 μ=%.3f', mu_example));
xlim([0.5, length(h_example)+0.5]);
grid on;

% ========== 5. 性能统计 ==========
fprintf('\n========== 性能统计结果 ==========\n');
fprintf('仿真信道数：%d\n', num_channels);
fprintf('平均信道容量：\n');
fprintf('  注水法：%.4f bits/s/Hz\n', mean(capacity_waterfill));
fprintf('  均匀分配：%.4f bits/s/Hz\n', mean(capacity_uniform));
fprintf('平均容量增益：%.2f%%\n', mean(gain_percentage));
fprintf('增益标准差：%.2f%%\n', std(gain_percentage));
fprintf('最大容量增益：%.2f%%\n', max(gain_percentage));
fprintf('最小容量增益：%.2f%%\n', min(gain_percentage));

% 计算容量增益与信道条件数的关系
cond_numbers = zeros(num_channels, 1);
for ch_idx = 1:num_channels
    % 生成信道矩阵并计算条件数
    H_temp = (randn(Nr, Nt) + 1i*randn(Nr, Nt)) / sqrt(2);
    cond_numbers(ch_idx) = cond(H_temp);
end

% 计算相关性
correlation = corrcoef(cond_numbers, gain_percentage);
fprintf('\n信道条件数与容量增益的相关系数：%.4f\n', correlation(1, 2));

三、注水法的高级变体

3.1 多用户MIMO注水法

function [P_alloc, capacity_total] = multiuser_waterfilling(H, P_total, sigma2, users)
% 多用户MIMO注水法功率分配
% 输入：
%   H - 多用户信道矩阵（每个用户的信道堆叠）
%   users - 用户数
% 输出：
%   P_alloc - 各用户功率分配矩阵
%   capacity_total - 系统和容量

[Nr_total, Nt] = size(H);
P_alloc = zeros(Nt, users);
capacity_users = zeros(users, 1);

% 总功率约束下的多用户注水
remaining_power = P_total;

while remaining_power > 1e-6
    % 为每个用户计算边际收益
    marginal_gains = zeros(users, 1);
    
    for u = 1:users
        % 提取当前用户信道
        start_idx = (u-1)*floor(Nr_total/users) + 1;
        end_idx = u*floor(Nr_total/users);
        H_u = H(start_idx:end_idx, :);
        
        % 计算当前功率下的信道容量
        [~, S_u, ~] = svd(H_u);
        lambda_u = diag(S_u).^2;
        current_cap = sum(log2(1 + P_alloc(:, u)' .* lambda_u / sigma2));
        
        % 增加少量功率计算边际增益
        delta_p = 1e-3;
        new_cap = sum(log2(1 + (P_alloc(:, u)' + delta_p) .* lambda_u / sigma2));
        marginal_gains(u) = (new_cap - current_cap) / delta_p;
    end
    
    % 将功率分配给边际增益最大的用户
    [~, best_user] = max(marginal_gains);
    
    % 为最佳用户分配一小部分功率
    delta_alloc = min(remaining_power, P_total * 0.01);
    P_alloc(:, best_user) = P_alloc(:, best_user) + delta_alloc;
    remaining_power = remaining_power - delta_alloc;
end

% 计算最终容量
for u = 1:users
    start_idx = (u-1)*floor(Nr_total/users) + 1;
    end_idx = u*floor(Nr_total/users);
    H_u = H(start_idx:end_idx, :);
    [~, S_u, ~] = svd(H_u);
    lambda_u = diag(S_u).^2;
    capacity_users(u) = sum(log2(1 + P_alloc(:, u)' .* lambda_u / sigma2));
end

capacity_total = sum(capacity_users);
end

3.2 频率选择性信道注水法

function [P_freq, capacity] = frequency_waterfilling(H_freq, P_total, sigma2, N_subcarriers)
% 频率选择性信道注水法
% H_freq: 每个子载波上的信道矩阵

P_freq = zeros(N_subcarriers, 1);
capacity_sub = zeros(N_subcarriers, 1);

% 收集所有子载波的增益
all_gains = [];
for f = 1:N_subcarriers
    [~, S_f, ~] = svd(H_freq{f});
    lambda_f = diag(S_f).^2;
    all_gains = [all_gains; lambda_f];
end

% 全局注水法
[p_all, mu] = waterfilling_basic(all_gains, P_total, sigma2);

% 重新分配到各子载波
idx = 1;
for f = 1:N_subcarriers
    [~, S_f, ~] = svd(H_freq{f});
    lambda_f = diag(S_f).^2;
    n_sub = length(lambda_f);
    
    p_f = p_all(idx:idx+n_sub-1);
    idx = idx + n_sub;
    
    P_freq(f) = sum(p_f);
    capacity_sub(f) = sum(log2(1 + p_f .* lambda_f / sigma2));
end

capacity = sum(capacity_sub);
end

四、性能分析与实际应用建议

4.1 注水法性能总结

场景	容量增益	计算复杂度	适用条件
高信噪比	显著（可达100%+）	O(N log N)	信道条件数大
低信噪比	有限（<10%）	O(N log N)	所有信道质量差
各向同性信道	可忽略	O(N log N)	奇异值接近相等
多用户系统	显著	O(KN log N)	用户间信道正交性强

4.2 实际系统调整建议

1. 不完全CSI（信道状态信息）情况：

function [p_robust] = robust_waterfilling(lambda_est, lambda_error, P_total, sigma2, confidence)
% 鲁棒注水法：考虑信道估计误差
% lambda_error: 估计误差的方差

% 最坏情况设计：考虑估计下界
lambda_min = max(0, lambda_est - confidence * sqrt(lambda_error));

% 在最坏情况下优化
[p_robust, ~] = waterfilling_basic(lambda_min, P_total, sigma2);

% 可选：增加功率余量
p_robust = p_robust * 0.95;  % 保留5%功率余量
end

2. 混合预编码系统中的注水法：

function [P_opt, F_RF, F_BB] = hybrid_precoding_waterfilling(H, P_total, N_RF)
% 混合预编码系统中的功率分配
% N_RF: RF链数量

[Nr, Nt] = size(H);

% 设计模拟预编码矩阵（基于信道主特征向量）
[U, ~, ~] = svd(H);
F_RF = U(:, 1:N_RF);

% 等效基带信道
H_eq = H * F_RF;

% 基带注水法
[~, S_eq, V_eq] = svd(H_eq);
lambda_eq = diag(S_eq).^2;

[P_opt, ~] = waterfilling_basic(lambda_eq, P_total, 1);

% 基带数字预编码
F_BB = V_eq(:, 1:N_RF) * diag(sqrt(P_opt));
end

4.3 算法局限性及改进方向

1. 复杂度问题：大规模MIMO（如128×128）需使用低复杂度近似：

function [p_approx] = approximate_waterfilling(lambda, P_total, sigma2, epsilon)
% 近似注水法：减少计算复杂度
% 基于阈值截断的小增益信道

% 忽略增益小于阈值的信道
threshold = epsilon * max(lambda);
active_idx = lambda > threshold;
lambda_active = lambda(active_idx);

% 仅对活跃信道应用注水法
[p_active, ~] = waterfilling_basic(lambda_active, P_total, sigma2);

% 重建完整功率向量
p_approx = zeros(length(lambda), 1);
p_approx(active_idx) = p_active;
end

2. 公平性考虑：多用户系统中可结合比例公平：

function [p_fair] = proportional_fair_waterfilling(lambda, P_total, user_weights)
% 比例公平注水法
% user_weights: 用户权重（基于历史吞吐量）

N = length(lambda);
h = 1 ./ lambda;

% 修改的注水法：考虑权重
for i = 1:100
    mu = (P_total + sum(h)) / N;  % 初始估计
    
    % 权重调整
    p_temp = max(0, mu - h ./ user_weights);
    
    if abs(sum(p_temp) - P_total) < 1e-6
        break;
    end
end

p_fair = p_temp;
end

参考代码 注水法优化解决MIMO信道功率优化 www.youwenfan.com/contentcso/96289.html

五、快速验证与测试

% 快速测试脚本
fprintf('=== 注水法快速测试 ===\n');

% 测试用例1：典型MIMO场景
H_test = [1.2+0.5i, 0.3-0.8i, -0.5+0.2i, 0.1+0.4i;
          -0.3+0.7i, 0.9-0.1i, 0.4+0.3i, -0.2+0.6i;
          0.5-0.2i, -0.4+0.3i, 1.1+0.9i, 0.3-0.1i;
          -0.1+0.4i, 0.2-0.5i, 0.6+0.3i, 0.8-0.2i];

[U_test, S_test, V_test] = svd(H_test);
lambda_test = diag(S_test).^2;

P_total_test = 10;  % 总功率
sigma2_test = 1;    % 噪声功率

[p_opt_test, mu_test, cap_test] = waterfilling_basic(lambda_test, P_total_test, sigma2_test);

fprintf('测试用例结果：\n');
fprintf('信道增益：%s\n', num2str(lambda_test', '%.3f  '));
fprintf('最优功率分配：%s\n', num2str(p_opt_test', '%.3f  '));
fprintf('注水电平：%.4f\n', mu_test);
fprintf('信道容量：%.4f bits/s/Hz\n', cap_test);

% 与均匀分配对比
p_uniform_test = P_total_test / length(lambda_test) * ones(size(lambda_test));
cap_uniform_test = sum(log2(1 + p_uniform_test .* lambda_test / sigma2_test));
gain_percentage = 100 * (cap_test / cap_uniform_test - 1);

fprintf('均匀分配容量：%.4f bits/s/Hz\n', cap_uniform_test);
fprintf('容量增益：%.2f%%\n\n', gain_percentage);

关键结论：注水法在MIMO信道功率优化中提供了理论最优解，特别适用于信道条件数大的场景。实际系统中需结合信道估计误差、硬件限制和公平性要求进行适当调整。