power allocation for mc cdma system using iterative water filling algorithm

My project is power allocation in MC CDMA using Iterative Water filling algorithm. i have the problem with plotting the graph for ber vs no users. Can u pls help me to plot. Here is my Code.

MC CDMA
clc;
close all;
clear all;
user=6;
nSubChannel = 6;
totalPower = 1e-5; % -20 dBm
channelStateInformation = random('rayleigh',1/0.6552,1,nSubChannel);
bandwidth = 1e6; % 1 MHz
noiseDensity = 1e-11; % -80 dBm
% total no of bits to be transmitted
% data=input('enter the datas which is need to transmit');
% n=length(data);
% no of sub carrier channels
c=6;
% bits per channels
bits=54;
n=324;
for i=1:n
data(i)= 2*round(rand)-1;
end
%first user data-------------------
[transmit carrier1 carrier2 carrier3 carrier4 carrier5 carrier6]=ch_ofdm(data,c,bits);
%second user data------------------
for i=1:n
data1(i)= 2*round(rand)-1;
end
[transmit1 carrier11 carrier12 carrier13 carrier14 carrier15 carrier16]=ch_ofdm(data1,c,bits);
% generating the noise
% p=rand(1,800)*2*pi;
p=rand*2*pi;
snr=10;
r=sqrt(-1*(1/snr*log(1 - rand)));
% no = 5*(r.* exp(j*p));
no = (r.* exp(j*p));
% value of alpha
al=rand+j*rand;
%al=1;
% Spreading channel with the alpha as the variable

for k=2:2:646
for l = 1:2
%al=round(rand)+j*round(rand)
rec(k+l)=transmit(k+l)+al*transmit(k-2+l);
end
end

rxdata=rec+ no ;
%-------
p=rand*2*pi;
snr=10;
r=sqrt(-1*(1/snr*log(1 - rand)));
% no = 5*(r.* exp(j*p));
no = (r.* exp(j*p));
% value of alpha
al=rand+j*rand;
%al=1;
% Spreading channel with the alpha as the variable

for k=2:2:646
for l = 1:2
%al=round(rand)+j*round(rand)
rec1(k+l)=transmit1(k+l)+al*transmit1(k-2+l);
end
end

rxdata1=rec1+ no ;
%--------------------------------------------------------------------------
[Capacity PowerAllocated] = ofdmwaterfilling(nSubChannel,totalPower,channelStateInformation,bandwidth,noiseDensity);
%-------------------------------------------------------------------------
%Reciever side 1st user
[demod,error]=rx_ofdm(rxdata,carrier1,carrier2,carrier3,carrier4,carrier5,carrier6,data)
ber=error/324
figure(3)
stem(data)
hold
stem(demod,'rx')
%-------------------------------------------------------------------------
%reciever side 2nd user
[demod1,error1]=rx_ofdm(rxdata1,carrier11,carrier12,carrier13,carrier14,carrier15,carrier16,data1)
ber=error1/324
figure(4)
stem(data1)
hold
stem(demod1,'rx')

OFDM WATER FILLING
function [shanonCapacity powerAllocated] = ofdmwaterfilling(nSubChannel,totalPower,channelStateInformation,bandwidth,noiseDensity)
%< Parameter Computation >

subchannelNoise = ...
noiseDensity*bandwidth/nSubChannel;
carrierToNoiseRatio = ...
channelStateInformation.^2/subchannelNoise;

initPowerAllo = ... (Formula 1)
(totalPower + sum(1./carrierToNoiseRatio))...
/nSubChannel - 1./carrierToNoiseRatio;



% < Iterative part of the algorithm >
while(length( find(initPowerAllo < 0 )) > 0 )
negIndex = find(initPowerAllo <= 0);
posIndex = find(initPowerAllo > 0);
nSubchannelRem = length(posIndex);
initPowerAllo(negIndex) = 0;
CnrRem = carrierToNoiseRatio(posIndex);
powerAlloTemp = (totalPower + sum(1./CnrRem))...
/nSubchannelRem - 1./CnrRem;
initPowerAllo(posIndex) = powerAlloTemp;
end

% < Output Computation >

% amount of power allocated to each subchannel
powerAllocated = initPowerAllo';
% total capacity of a channel based on shanon theory
shanonCapacity = bandwidth/nSubChannel * ...
sum(log2(1 + initPowerAllo.*carrierToNoiseRatio));

% <Graphical Observation>


f1 = figure(1);
clf;
set(f1,'Color',[1 1 1]);
bar((initPowerAllo + 1./carrierToNoiseRatio),1,'r')
hold on;
bar(1./carrierToNoiseRatio,1);
xlabel('subchannel indices');
title('Water filling algorithm')

legend('amount of power allocated to each subchannel',...
'Noise to Carrier Ratio')


CH OFDM
function [transmit carrier1 carrier2 carrier3 carrier4 carrier5 carrier6]=ch_ofdm(data,c,bits)
n=324;
% Converting the series into parallel for the channels
s = reshape(data,c,bits);
% BPSK modulation for a single channel
signal1 = s(1,;
%first expand the bit stream
exdata1=[];
exdata2=[];
exdata3=[];
exdata4=[];
exdata5=[];
exdata6=[];
for i=1:length(signal1)
for rep=1:2
exdata1= [exdata1 signal1(i)];
end
end

signal2 = s(2,;
%first expand the bit stream
exdata=[];
for i=1:length(signal2)
for rep=1:2
exdata2= [exdata2 signal2(i)];
end
end

signal3= s(3,;
%first expand the bit stream
exdata=[];
for i=1:length(signal3)
for rep=1:2
exdata3= [exdata3 signal3(i)];
end
end

signal4= s(4,;
%first expand the bit stream
exdata=[];
for i=1:length(signal4)
for rep=1:2
exdata4= [exdata4 signal4(i)];
end
end

signal5= s(5,;
%first expand the bit stream
exdata=[];
for i=1:length(signal5)
for rep=1:2
exdata5= [exdata5 signal5(i)];
end
end

signal6= s(6,;
%first expand the bit stream
exdata=[];
for i=1:length(signal6)
for rep=1:2
exdata6= [exdata6 signal6(i)];
end
end

% Bpsk modulation
m=10*n;
% Generating the carrier signal
ts=.1;
tp=1:ts:11.79;
carrier1=cos(2*pi*tp);
% Generating the modulated signal 1
bpsk_sig1=exdata1.*carrier1;
% Generating the modulated signal 2
carrier2=cos(4*pi*tp);
bpsk_sig2=exdata2.*carrier2;
% Generating the modulated signal 3
carrier3=cos(6*pi*tp);
bpsk_sig3=exdata3.*carrier3;
% Generating the modulated signal 4
carrier4=cos(8*pi*tp);
bpsk_sig4=exdata4.*carrier4;
% Generating the modulated signal 5
carrier5=cos(10*pi*tp);
bpsk_sig5=exdata5.*carrier5;
carrier6=cos(12*pi*tp);
% Generating the modulated signal
bpsk_sig6=exdata6.*carrier6;
% taking the iFFT of each of these signals
if_sig1=ifft(bpsk_sig1);
if_sig2=ifft(bpsk_sig2);
if_sig3=ifft(bpsk_sig3);
if_sig4=ifft(bpsk_sig4);
if_sig5=ifft(bpsk_sig5);
if_sig6=ifft(bpsk_sig6);

fin(1,=if_sig1;
fin(2,=if_sig2;
fin(3,=if_sig3;
fin(4,=if_sig4;
fin(5,=if_sig5;
fin(6,=if_sig6;

transmit=reshape(fin,1,648);


RECEIVER OFDM
function [demod,error]=rx_ofdm(rxdata,carrier1,carrier2,carrier3,carrier4,carrier5,carrier6,data)
% Converting from serial to parallel
myrec=reshape(rxdata,6,108);

rxdata1=fft(myrec(1,);
rxdata2=fft(myrec(2,);
rxdata3=fft(myrec(3,);
rxdata4=fft(myrec(4,);
rxdata5=fft(myrec(5,);
rxdata6=fft(myrec(6,);

% taking the FFT


%begin demodulation
%first multiply recieved bitstream by cosine wave with carrier frequency


uncarry1=rxdata1.*carrier1;
uncarry2=rxdata2.*carrier2;
uncarry3=rxdata3.*carrier3;
uncarry4=rxdata4.*carrier4;
uncarry5=rxdata5.*carrier5;
uncarry6=rxdata6.*carrier6;

%plot(uncarry)
%demodulate by integrating
dec1=[];
dec2=[];
dec3=[];
dec4=[];
dec5=[];
dec6=[];
for inc=1:2:length(uncarry1)
dec=trapz(inc:inc+1,uncarry1(inc:inc+1));
dec1=[dec1 dec];
end
%2
for inc=1:2:length(uncarry2)
dec=trapz(inc:inc+1,uncarry2(inc:inc+1));
dec2=[dec2 dec];
end
%3
for inc=1:2:length(uncarry3)
dec=trapz(inc:inc+1,uncarry3(inc:inc+1));
dec3=[dec3 dec];
end
%4
for inc=1:2:length(uncarry4)
dec=trapz(inc:inc+1,uncarry4(inc:inc+1));
dec4=[dec4 dec];
end
%5
for inc=1:2:length(uncarry5)
dec=trapz(inc:inc+1,uncarry5(inc:inc+1));
dec5=[dec5 dec];
end
%6
for inc=1:2:length(uncarry6)
dec=trapz(inc:inc+1,uncarry6(inc:inc+1));
dec6=[dec6 dec];
end
final_rec(1,=dec1;
final_rec(2,=dec2;
final_rec(3,=dec3;
final_rec(4,=dec4;
final_rec(5,=dec5;
final_rec(6,=dec6;

fin_rec_parallel=reshape(final_rec,1,324)
%make decision with a threshold of zero
demod=[];
for i=1:length(fin_rec_parallel)
if fin_rec_parallel(i)>0
demod=[demod 1];
else
demod=[demod -1];
end
end

%stem(demod)

%calculate errors
error=0;
for i=1:length(demod)
if data(i)~=demod(i)
error=error+1;
end
end
error;