多用户资源分配

clc;
clear all;
K=2;         %认知用户数
N=64;        %系统总的子载波数
cmax=6;      %每个子载波平均所能携带的最大比特数目
cmin=120;    %满足自身QoS需求所需要的最小信道容量
B=1e6;       %总的信道带宽
Pe=1e-7;     %误码率
P_av=1; 
SNR_av=0; 
Noise=P_av./10.^（SNR_av./10）; %噪声功率
N0=Noise./（B/N）;              %噪声功率谱密度
% N0=1e-6;
% PL=B/N*N0*（10.^（SNR_av./10））;
PL=0.0078;%保证认知用户能够正常通信的功率需求
Pt=1;                         %所有认知用户的总功率上限
PT=0.5;                       %各认知用户的总功率上限
l=1:cmax;                     %子载波每次分配增加的比特数

% 子载波分配
% h = rand（KN）;  %信道系数
% h = random（‘rayleigh‘1KN）; %信道增益
% H（k）= h（k:）.^2;              %信道增益的平方
Nk_initial=zeros（K1）;          %各用户初始分配的子载波数
h=zeros（KN）;                   %信道增益
for k=1:K
    Nk_initial（k）=floor（cmin/cmax）;%每个用户最少占用的子载波数目
    if sum（Nk_initial）>N           %如果所有用户需要的最少的子载波数目之和仍大于所有的子载波数目，则报错
        error（‘There are two many users!\n‘）;
    end
    for n=1:N
        h（kn）=1/sqrt（2）*（randn（1） + 1j*randn（1））;  %信道增益
    end
end
H=abs（h）;
% sort（H2‘descend‘）;%对各用户的信道增益进行降序排列
Htemp = H.^2;         %信道系数的临时记录值，若该子载波被占用，则信道系数设置为0

P=zeros（K1）;    %各用户的功率需求
omiga=cell（K1）; %前20个子载波组成的集合，初始化为空集
f=zeros（KN）;    %用于存储接收到c个比特时所需的功率
p=zeros（KN）;    %用于存储认知用户k在第n个子载波上消耗的功率   
Htemp1=Htemp;
for k=1:K
    m=1;
    while m<=Nk_initial（k）
        Hmax=max（max（Htemp1（k:）））;%找出最大的信道系数
        [kh index]=find（Htemp1（k:）==Hmax）;%找出该信道系数的坐标
        %此时并计算认知用户k的第n个子载波的功率
        f（km）=N0/3*（qfuncinv（Pe/4）*qfuncinv（Pe/4））*（2^（cmin./Nk_initial（k））-1）;%MQAM
        p（km）=f（km）/Htemp1（kindex）;
        omiga{k}=union（omiga{k}index）;
        Htemp1（kindex）=0;
        m=m+1;
    end
    P（k）=sum（p（k:）2）;
end
% P=sum（p2）;
Ptotal=sum（P）;          %总功率需求
Ntotal=sum（Nk_initial）; %分配的子载波数总和

Pp=zeros（K1）;
PP=zeros（K1）;
ff=zeros（KN）;   
pp=zeros（KN）;
CP=zeros（K1）;%记录各用户每多分配一个子载波时所减小的功率
NK=Nk_initial;
Htemp2=Htemp;
while Ptotal>Pt*（Ntotal/N）
%     Htemp2=Htemp;
    for k=1:K
        m=1;
        while m<=NK（k）
            Hmax2=max（max（Htemp2（k:）））;%找出最大的信道系数
            [kh1 index1]=find（Htemp2（k:）==Hmax2）;%找出该信道系数的坐标
%             [Hmax nm]=max（Htemp2（k:））;%找出最大的信道系数  
            %此时并计算认知用户k的第n个子载波的功率
            f（km）=N0/3*（qfuncinv（Pe/4）*qfuncinv（Pe/4））*（2^（cmin./NK（k））-1）;%MQAM
%             p（km）=f（km）/Hmax;
            p（km）=f（km）/Htemp2（kindex1）;
            ff（km）=N0/3*（qfuncinv（Pe/4）*qfuncinv（Pe/4））*（2^（cmin./（NK（k）+1））-1）;%MQAM
%             pp（km）=ff（km）/Hmax;
            pp（km）=ff（km）/Htemp2（kindex1）;
            omiga{k}=union（omiga{k}index1）;
            Htemp2（kindex1）=0;
%             omiga{k}=union（omiga{k}nm）;
%             Htemp2（knm）=0;
            m=m+1;
        end
        Pp（k）=sum（p（k:）2）;
        PP（k）=sum（pp（k:）2）;
    end
        Htemp4=Htemp2;
        for kk=1:K
            Hmax3=max（max（Htemp4（kk:）））;%找出最大的信道系数
            [kh2 index2]=find（Htemp4（kk:）==Hmax3）;%找出该信道系数的坐标
            %此时并计算认知用户k的在原来基础上增加一个子载波的功率