OFDM迭代注水算法

% function P0=IPW（NPtWN0CSIMG）
% 迭代注水功率分配算法实质是多次执行相同的运算而已
% SU所占用频带的带宽为5MHz。这一频带被均匀的分成4个子信道。
% 每个子信道对应于一个PU的授权频谱。子载波的总数为64即每个子信道包含16个子载波。
% N:     次用户所占用的子载波总数
% Pt：   次用户发射总功率
% W:     次用户所占用的信道带宽
% N0：   信道噪声功率
% CSI:   信道状态信息（一般为传递函数）
% M:     次用户所分子信道数
% G:     子信道的发射功率约束
clear all;
close all;
clc;
format compact;

%%%---输入参数---%%%
n=1:64;            % 次用户所占子载波总数
N=length（n）;       % 子载波数
W=5*1e6;           % 次用户所占用的信道带宽（5MHz）
Pt=640;            % 次用户发射总功率（w）
m=1:4;             % 次用户所分子信道数序号
M=length（m）;       % 次用户所分子信道数
numc=N/M;          % 每个子信道所分的子载波数
W1=W/M;            % 每个子信道所占带宽
% f0=0.15*1e6;     % 子载波间隔

%%%---初始化有关参量---%%%
CNR=zeros（1N）;      % 载噪比CNR初始置零
CSI=zeros（1N）;      % 子载波信道传输函数初始置零
N0=0;                % 信道噪声功率初始置零

%%%---高斯白噪声信道---%%%
% % 零均值复高斯随机变量
% R_u=randn（1N）; %实部
% I_u=randn（1N）; %虚部
% j=sqrt（-1）;     %虚指数
% u=（R_u+j*I_u）/sqrt（2）;
% u = wgn（1N0‘complex‘）;    % 高斯白噪声信道 
% c=0.03;                      % 瑞利衰落信道参数
% h=exp（-c*n）.*u;              % 瑞利衰落信道的信道冲激响应
% % 加入信噪比为20dB的加性高斯白噪声的瑞利衰落信道的信道冲激响应
% ha=awgn（h20）;
ha=random（‘rayleigh‘11N）;
% 信道状态信息，服从瑞利分布，参数为1的1行N列随机数
CSI=abs（ha）;         % 信道初始化信息（传输函数）

Np=12.8*1e-6;        % 信道噪声功率谱密度
N0=Np*W/N;           % 信道噪声功率（w）（也可直接设置N0）

%%%%%%%%%%%%------------算法第一步-------------%%%%%%%%%%%%

k=1;
p=zeros（1N）;        % 初始化子载波所分配的功率
P（k）=0;              % 次用户发射总功率初始置零

%%%---求出集合A中子载波所分配的功率---%%%
for i=1:N;
    CNR=CSI.^2/N0;   % 信道载波噪声功率比（信道载噪比）
    P（k）=Pt;         % 发射功率初始化信息
    P0=（P（k）+sum（1./CNR））/N-1./CNR;  % 初始化功率分配
    while（length（find（P0<0））>0）;   
        negIndex =find（P0<=0）;       % 初始功率P0小于等于0的子载波序号
        posIndex =find（P0>0）;        % 初始功率P0大于0的子载波序号 
        P0（negIndex） = 0;            % 初始功率P0小于等于0的子载波所分配功率置零
        Nr =length（posIndex）;        % 初始功率P0>0的子载波个数
        Cnr =CNR（posIndex）;          % 初始功率P0>0的子载波的载噪比
        Ptemp =（P（k） + sum（1./Cnr））/Nr - 1./Cnr; 
        P0（posIndex） = Ptemp;        % 初始功率P0大于0的子载波所分配功率
    end
    p（i）=P0（i）;       % 每个子载波利用经典注水算法所分配的功率
end             % 此for语句是可以省略的，只是为了更好说明子载波所分配的功率

% p             % 显示经典注水算法所分配的功率向量

%%%%%%------经典注水算法所分配的功率（柱状图）------%%%%%%
f1 = figure;
     clf;       % 清除目前窗口
     set（f1‘Color‘[1 1 1]）;
     bar（（p+1./CNR）1‘r‘）;
     hold on;     
     bar（1./CNR1）;
     xlabel（‘子载波序号‘）; 
     ylabel（‘所分配的功率‘）;
     title（‘经典注水功率分配算法‘）      
     legend（‘分配功率Pi‘ ‘信道特性1/hi‘）;
     hold off

%%%%%%%%%%%%------------算法第二步-------------%%%%%%%%%%%%

J=zeros（MN）;             % 初始化子信道外子载波的功率泄露矩阵
for i=1:M;
    for j=numc*（i-1）+1:numc*i;
        J（ij）=1;         % 不考虑子信道外子载波的功率泄露矩阵
    end
end

%%%---考虑子信道外子载波的功率泄露矩阵---%%%
% for i=1:5;
%     m（i）=16*（i-1）+1;
% end
% 
% syms x     %标注变量
% for j=1:64;
%     J1（j）=int（（（sin（pi*（x-j）））^2）/（（pi*（x-j））^2）xm（1）-0.5m（2）-0.5）;
% end
%     J1=eval（J1）;     % 转换成数值型
% for j=1:64;
%     J2（j）=int（（（sin（pi*（x-j）））^2）/（（pi*（x-j））^2）