%?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）