%?function?P0=IPW（NPtWN0CSIMG）
%?迭代注水功率分配算法實(shí)質(zhì)是多次執(zhí)行相同的運(yùn)算而已
%?SU所占用頻帶的帶寬為5MHz。這一頻帶被均勻的分成4個(gè)子信道。
%?每個(gè)子信道對(duì)應(yīng)于一個(gè)PU的授權(quán)頻譜。子載波的總數(shù)為64即每個(gè)子信道包含16個(gè)子載波。
%?N:?????次用戶(hù)所占用的子載波總數(shù)
%?Pt：???次用戶(hù)發(fā)射總功率
%?W:?????次用戶(hù)所占用的信道帶寬
%?N0：???信道噪聲功率
%?CSI:???信道狀態(tài)信息（一般為傳遞函數(shù)）
%?M:?????次用戶(hù)所分子信道數(shù)
%?G:?????子信道的發(fā)射功率約束
clear?all;
close?all;
clc;
format?compact;

%%%---輸入?yún)?shù)---%%%
n=1:64;????????????%?次用戶(hù)所占子載波總數(shù)
N=length（n）;???????%?子載波數(shù)
W=5*1e6;???????????%?次用戶(hù)所占用的信道帶寬（5MHz）
Pt=640;????????????%?次用戶(hù)發(fā)射總功率（w）
m=1:4;?????????????%?次用戶(hù)所分子信道數(shù)序號(hào)
M=length（m）;???????%?次用戶(hù)所分子信道數(shù)
numc=N/M;??????????%?每個(gè)子信道所分的子載波數(shù)
W1=W/M;????????????%?每個(gè)子信道所占帶寬
%?f0=0.15*1e6;?????%?子載波間隔

%%%---初始化有關(guān)參量---%%%
CNR=zeros（1N）;??????%?載噪比CNR初始置零
CSI=zeros（1N）;??????%?子載波信道傳輸函數(shù)初始置零
N0=0;????????????????%?信道噪聲功率初始置零

%%%---高斯白噪聲信道---%%%
%?%?零均值復(fù)高斯隨機(jī)變量
%?R_u=randn（1N）;?%實(shí)部
%?I_u=randn（1N）;?%虛部
%?j=sqrt（-1）;?????%虛指數(shù)
%?u=（R_u+j*I_u）/sqrt（2）;
%?u?=?wgn（1N0‘complex‘）;????%?高斯白噪聲信道?
%?c=0.03;??????????????????????%?瑞利衰落信道參數(shù)
%?h=exp（-c*n）.*u;??????????????%?瑞利衰落信道的信道沖激響應(yīng)
%?%?加入信噪比為20dB的加性高斯白噪聲的瑞利衰落信道的信道沖激響應(yīng)
%?ha=awgn（h20）;
ha=random（‘rayleigh‘11N）;
%?信道狀態(tài)信息，服從瑞利分布，參數(shù)為1的1行N列隨機(jī)數(shù)
CSI=abs（ha）;?????????%?信道初始化信息（傳輸函數(shù)）

Np=12.8*1e-6;????????%?信道噪聲功率譜密度
N0=Np*W/N;???????????%?信道噪聲功率（w）（也可直接設(shè)置N0）

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

k=1;
p=zeros（1N）;????????%?初始化子載波所分配的功率
P（k）=0;??????????????%?次用戶(hù)發(fā)射總功率初始置零

%%%---求出集合A中子載波所分配的功率---%%%
for?i=1:N;
????CNR=CSI.^2/N0;???%?信道載波噪聲功率比（信道載噪比）
????P（k）=Pt;?????????%?發(fā)射功率初始化信息
????P0=（P（k）+sum（1./CNR））/N-1./CNR;??%?初始化功率分配
????while（length（find（P0<0））>0）;???
????????negIndex?=find（P0<=0）;???????%?初始功率P0小于等于0的子載波序號(hào)
????????posIndex?=find（P0>0）;????????%?初始功率P0大于0的子載波序號(hào)?
????????P0（negIndex）?=?0;????????????%?初始功率P0小于等于0的子載波所分配功率置零
????????Nr?=length（posIndex）;????????%?初始功率P0>0的子載波個(gè)數(shù)
????????Cnr?=CNR（posIndex）;??????????%?初始功率P0>0的子載波的載噪比
????????Ptemp?=（P（k）?+?sum（1./Cnr））/Nr?-?1./Cnr;?
????????P0（posIndex）?=?Ptemp;????????%?初始功率P0大于0的子載波所分配功率
????end
????p（i）=P0（i）;???????%?每個(gè)子載波利用經(jīng)典注水算法所分配的功率
end?????????????%?此for語(yǔ)句是可以省略的，只是為了更好說(shuō)明子載波所分配的功率

%?p?????????????%?顯示經(jīng)典注水算法所分配的功率向量

%%%%%%------經(jīng)典注水算法所分配的功率（柱狀圖）------%%%%%%
f1?=?figure;
?????clf;???????%?清除目前窗口
?????set（f1‘Color‘[1?1?1]）;
?????bar（（p+1./CNR）1‘r‘）;
?????hold?on;?????
?????bar（1./CNR1）;
?????xlabel（‘子載波序號(hào)‘）;?
?????ylabel（‘所分配的功率‘）;
?????title（‘經(jīng)典注水功率分配算法‘）??????
?????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?????%標(biāo)注變量
%?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）;?????%?轉(zhuǎn)換成數(shù)值型
%?for?j=1:64;
%?????J2（j）=int（（（sin（pi*（x-j）））^2）/（（pi*（x-j））^2）