optimization in Cognitive radio network (spectrum sensing, total error rate)

clc; 
close all;
clear all;
%This program is for optimization of spectrum sensing in 
%Cognitive radio network.
N=20;
j=1;
tt=[];
err2=[];
Pmi=[];
Pdc=[];
error=[];
err1=[];
K=10;
snr=10;
Qd=0;
Qf=0;
tt=10:0.5:60;
vec=['-+','-o','-v','-d','->','-x','-s','-<','-*','-^'];

for n=1:1:10
    
s=ones(1,N);
w=randn(1,N); 

u=N/2;                    %Time-delay bandwidth product

for t=10:0.5:60 
  
Qd=0;
Qf=0;
SNR=10^(snr/10);   %for linear scale 

a=sqrt(2*SNR);
b=sqrt(t);



Pd = marcumq(a,b,u );         % AVG. PROB OF DETECTION(computes the generalized Marcum Q)
                 
Pf = gammainc((t/2),u,'upper');% AVG. PROB OF FALSE ALARM(compute incompelete gamma function)

Pm=1-Pd;                %AVG. PROB OF MISSED DETECTION OVER AWGN
for l=n:1:K
Qd=Qd+(factorial(K)*(Pd^l)*((1-Pd)^(K-l))/(factorial(l)*factorial(K-l)));
Qf=Qf+(factorial(K)*(Pf^l)*((1-Pf)^(K-l))/(factorial(l)*factorial(K-l)));
end

Qm=1-Qd;
err=Qf+Qm;
err1=[err1 err];
end

end
l=1;
i=1;
for j=1:1:10
semilogy(tt,err1(i:i+100),vec(l:l+1),'LineWidth',1.5)
i=i+101;
l=l+2;
hold on;
end
grid on;
ylabel('Total Error rate');
xlabel('Threshold');


%----------------------Energy Detection----------------------------------------
n=5;
rel=10000;
tt1=10:0.5:60;
er1=[];
for t=10:0.5:60
Pdc=0;
Pfc=0;
Qd=0;
Qf=0;
Qm=0;

for i=1:1:rel
SNR=10;

snr=10^(SNR/10);

    
s=ones(1,N);
w=randn(1,N); 
vari=var(w);               %variance of noise

Es=sum(s.^2);

N02=(Es)/(2*snr);   

x1=s+w; 
x2=w;
W=1;                    %Time-delay bandwidth product

E0=(sum(x2.^2))/((W*N02));

E1=(sum(x1.^2))/((W*N02));


if E1>t
    Pdc=Pdc+1;
else
end
if E0>t
    Pfc=Pfc+1;
else
end
end
Pd=Pdc/rel;
Pf=Pfc/rel;
   
  
for l=n:1:K
Qd=Qd+(factorial(K)*(Pd^l)*((1-Pd)^(K-l))/(factorial(l)*factorial(K-l)));
Qf=Qf+(factorial(K)*(Pf^l)*((1-Pf)^(K-l))/(factorial(l)*factorial(K-l)));
end
Qm=1-Qd;
er=Qf+Qm;
er1=[er1 er];
end
hold on;
semilogy(tt1,er1,'*r')
grid on;

ylabel('Total Error rate');
xlabel('Threshold');
legend('n=1','n=2','n=3','n=4','n=5','n=6','n=7','n=8','n=9','n=10','n=5 by modelling');

Low Energy Adaptive Clustering Hierarchy protocol (LEACH) in WSN & path_delay_calculation in Wireless Sensor Networks GUI

LEACH is a hierarchical protocol in which most nodes transmit to cluster heads, and the cluster heads aggregate and compress the data and forward it to the base station(sink). Each node uses a stochastic algorithm at each round to determine whether it will become a cluster head in this round. LEACH assumes that each node has a radio powerful enough to directly reach the base station or the nearest cluster head, but that using this radio at full power all the time would waste energy. 
Nodes that have been cluster heads cannot become cluster heads again for P rounds, where P is the desired percentage of cluster heads. Thereafter, each node has a 1/P probability of becoming a cluster head in each round. At the end of each round, each node that is not a cluster head selects the closest cluster head and joins that cluster. The cluster head then creates a schedule for each node in its cluster to transmit its data. 
All nodes that are not cluster heads only communicate with the cluster head in a TDMA fashion, according to the schedule created by the cluster head. They do so using the minimum energy needed to reach the cluster head, and only need to keep their radios on during their time slot. 
LEACH also uses CDMA so that each cluster uses a different set of CDMA codes, to minimize interference between clusters.

 MATLAB release MATLAB 8.0 (R2012b)

Other requirements Windows x64

clear;
%%%%%%%%%%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Field Dimensions - x and y maximum (in meters)/ Tarife Size mohite shabake(m)
xm=100;
ym=100;
%x and y Coordinates of the Sink  / 
sink.x=0.5*xm;
sink.y=0.5*ym;

%Number of Nodes in the field  / Tedade Node haye shabake
n=100

%Optimal Election Probability of a node to become cluster head/ Ehtemale Entekhab Node be onvane Cluster Head
p=0.1;
%Energy Model (all values in Joules)/ Energy ha bar hasbe Joule
%Initial Energy / Energy Avaliye
Eo=0.5;
%Eelec=Etx=Erx 
ETX=50*0.000000001;
ERX=50*0.000000001;
%Transmit Amplifier types / Ghodrate Ersal
Efs=10*0.000000000001;
Emp=0.0013*0.000000000001;
%Data Aggregation Energy/ Energy Masrafi Tajmie Dade
EDA=5*0.000000001;
%Values for Hetereogeneity 
%Percentage of nodes than are advanced 
m=0.1;
%\alpha
a=1;
%maximum number of rounds/ Max tedade round ha
rmax=200
%%%%%%%%%%%%%%%%%%%%%%%%% END OF PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%
%Computation of do/ 
do=sqrt(Efs/Emp);
%Creation of the random Sensor Network/ Tolide Randome shabake
figure(1);
for i=1:1:n
    S(i).xd=rand(1,1)*xm;
    XR(i)=S(i).xd;
    S(i).yd=rand(1,1)*ym;
    YR(i)=S(i).yd;
    S(i).G=0;
    %initially there are no cluster heads only nodes/ Dar ebteda hich Cluster Head i mojud nist
    S(i).type='N';
   
    temp_rnd0=i;
    %Random Election of Normal Nodes/ Entekhabe Tasadofi Node ha 
    if (temp_rnd0>=m*n+1) 
        S(i).E=Eo;
        S(i).ENERGY=0;
        plot(S(i).xd,S(i).yd,'o');
        hold on;
    end
    %Random Election of Advanced Nodes/ Entekhab Tasadofie CH ha
    if (temp_rnd0<m*n+1)  
        S(i).E=Eo*(1+a)
        S(i).ENERGY=1;
        plot(S(i).xd,S(i).yd,'+');
        hold on;
    end
end
S(n+1).xd=sink.x;
S(n+1).yd=sink.y;
plot(S(n+1).xd,S(n+1).yd,'x');
  
      
%First Iteration
figure(1);

%counter for CHs/ Tedade Cluster Head ha 
countCHs=0;
%counter for CHs per round/ Tedade CH haye har Round
rcountCHs=0;
cluster=1;

countCHs;
rcountCHs=rcountCHs+countCHs;
flag_first_dead=0;

for r=0:1:rmax
    r

  %Operation for epoch/ Formule entekhabe CH
  if(mod(r, round(1/p) )==0)
    for i=1:1:n
        S(i).G=0;
        S(i).cl=0;
    end
  end

hold off;

%Number of dead nodes/ Tedade Node haye morde dar kol
dead=0;
%Number of dead Advanced Nodes/ Tedade Node haye CH morde
dead_a=0;
%Number of dead Normal Nodes/ Tedade Node haye morde mamuli
dead_n=0;

%counter for bit transmitted to Bases Station and to Cluster Heads/ Tedade packet haye ersali be BS 
packets_TO_BS=0;
packets_TO_CH=0;
%counter for bit transmitted to Bases Station and to Cluster Heads /Tedade packet haye Ersali be BS dar har round
%per round
PACKETS_TO_CH(r+1)=0;
PACKETS_TO_BS(r+1)=0;

figure(1);

for i=1:1:n
    %checking if there is a dead node/ Check kardane zende budane Node ha 
    if (S(i).E<=0)
        plot(S(i).xd,S(i).yd,'red .');
        dead=dead+1;
        if(S(i).ENERGY==1)
            dead_a=dead_a+1;
        end
        if(S(i).ENERGY==0)
            dead_n=dead_n+1;
        end
        hold on;    
    end
    if S(i).E>0
        S(i).type='N';
        if (S(i).ENERGY==0)  
        plot(S(i).xd,S(i).yd,'o');
        end
        if (S(i).ENERGY==1)  
        plot(S(i).xd,S(i).yd,'+');
        end
        hold on;
    end
end
plot(S(n+1).xd,S(n+1).yd,'x');


STATISTICS(r+1).DEAD=dead;
DEAD(r+1)=dead;
DEAD_N(r+1)=dead_n;
DEAD_A(r+1)=dead_a;

%When the first node dies/ Pas az marge avalin Node
if (dead==1)
    if(flag_first_dead==0)
        first_dead=r
        flag_first_dead=1;
    end
end

countCHs=0;
cluster=1;
for i=1:1:n
   if(S(i).E>0)
   temp_rand=rand;     
   if ( (S(i).G)<=0)

 %Election of Cluster Heads/ Entekhabe CH
 if(temp_rand<= (p/(1-p*mod(r,round(1/p)))))
            countCHs=countCHs+1;
            packets_TO_BS=packets_TO_BS+1;
            PACKETS_TO_BS(r+1)=packets_TO_BS;
            
            S(i).type='C';
            S(i).G=round(1/p)-1;
            C(cluster).xd=S(i).xd;
            C(cluster).yd=S(i).yd;
            plot(S(i).xd,S(i).yd,'k*');
            
            distance=sqrt( (S(i).xd-(S(n+1).xd) )^2 + (S(i).yd-(S(n+1).yd) )^2 );
            C(cluster).distance=distance;
            C(cluster).id=i;
            X(cluster)=S(i).xd;
            Y(cluster)=S(i).yd;
            cluster=cluster+1;
            
            %Calculation of Energy dissipated/  Mohasebe energy masrafi
            distance;
            if (distance>do)
                S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            end
            if (distance<=do)
                S(i).E=S(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
        end     
    
    end
  end 
end

STATISTICS(r+1).CLUSTERHEADS=cluster-1;
CLUSTERHS(r+1)=cluster-1;

%Election of Associated Cluster Head for Normal Nodes/ Entekhabe CH marbute baraye Node haye mamuli
for i=1:1:n
   if ( S(i).type=='N' && S(i).E>0 )
     if(cluster-1>=1)
       min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );
       min_dis_cluster=1;
       for c=1:1:cluster-1
           temp=min(min_dis,sqrt( (S(i).xd-C(c).xd)^2 + (S(i).yd-C(c).yd)^2 ) );
           if ( temp<min_dis )
               min_dis=temp;
               min_dis_cluster=c;
           end
       end
       
       %Energy dissipated by associated Cluster Head/ Mohasebe energy masrafi CH ha 
            min_dis;
            if (min_dis>do)
                S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S(i).E=S(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
        %Energy dissipated/ Masrafe energy kol
        if(min_dis>0)
          S(C(min_dis_cluster).id).E = S(C(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
         PACKETS_TO_CH(r+1)=n-dead-cluster+1; 
        end

       S(i).min_dis=min_dis;
       S(i).min_dis_cluster=min_dis_cluster;
           
   end
 end
end
hold on;

countCHs;
rcountCHs=rcountCHs+countCHs;
end

%  DEAD  : a rmax x 1 array of number of dead nodes/round          %
%  DEAD_A : a rmax x 1 array of number of dead Advanced nodes/round       %
%  DEAD_N : a rmax x 1 array of number of dead Normal nodes/round                     %
%  CLUSTERHS : a rmax x 1 array of number of Cluster Heads/round                      %
%  PACKETS_TO_BS : a rmax x 1 array of number packets send to Base Station/round      %
%  PACKETS_TO_CH : a rmax x 1 array of number of packets send to ClusterHeads/round   %
%  first_dead: the round where the first node died  

WSN protocol implementation leach program GUI

close all;
clear;
clc;

%-------------------------------
%Number of Nodes in the field
 n=200;
%n=input('Enter the number of nodes in the space : '); 
%Energy Model (all values in Joules)
%Initial Energy 
Eo=0.1;
%Eo=input('Enter the initial energy of sensor nJ : ');
%Field Dimensions - x and y maximum (in meters)
% xm=input('Enter x value for area plot : ');
% ym=input('Enter y value for area plot : ');  
xm=100;
ym=100;

%x and y Coordinates of the Sink
sink.x=1.5*xm;
sink.y=0.5*ym;

%Optimal Election Probability of a node
%to become cluster head
p=0.2;

%Eelec=Etx=Erx
ETX=50*0.000000001;
ERX=50*0.000000001;
%Transmit Amplifier types
Efs=10*0.000000000001;
Emp=0.0013*0.000000000001;
%Data Aggregation Energy
EDA=5*0.000000001;

%Values for Hetereogeneity
%Percentage of nodes than are advanced
m=0.5;
%\alpha
a=1;

%maximum number of rounds
%rmax=input('enter the number of iterations you want to run : '); 
rmax=50;

%------------------

%Computation of do
do=sqrt(Efs/Emp);

%Creation of the random Sensor Network
figure(1);
hold off;
for i=1:1:n
    S(i).xd=rand(1,1)*xm;
    XR(i)=S(i).xd;
    S(i).yd=rand(1,1)*ym;
    YR(i)=S(i).yd;
    S(i).G=0;
    %initially there are no cluster heads only nodes
    S(i).type='N';
   
    temp_rnd0=i;
    %Random Election of Normal Nodes
    if (temp_rnd0>=m*n+1) 
        S(i).E=Eo;
        S(i).ENERGY=0;
       plot(S(i).xd,S(i).yd,'o-r');
        hold on;
    end
    %Random Election of Advanced Nodes
    if (temp_rnd0<m*n+1)  
        S(i).E=Eo*(1+a);
        S(i).ENERGY=1;
       plot(S(i).xd,S(i).yd,'+');
        hold on;
    end
end

S(n+1).xd=sink.x;
S(n+1).yd=sink.y;
plot(S(n+1).xd,S(n+1).yd,'o', 'MarkerSize', 12, 'MarkerFaceColor', 'r');
figure(1);  
% figure(1)
%  plot(o1,o2,'^','LineWidth',1, 'MarkerEdgeColor','k', 'MarkerFaceColor','y', 'MarkerSize',12);
%    hold on
%First Iteration
%counter for CHs
countCHs=0;
%counter for CHs per round
rcountCHs=0;
cluster=1;

countCHs;
rcountCHs=rcountCHs+countCHs;
flag_first_dead=0;

for r=0:1:rmax
    r;

  %Operation for epoch
  if(mod(r, round(1/p) )==0)
    for i=1:1:n
        S(i).G=0;
        S(i).cl=0;
    end
  end

hold off;

%Number of dead nodes
dead=0;
%Number of dead Advanced Nodes
dead_a=0;
%Number of dead Normal Nodes
dead_n=0;

%counter for bit transmitted to Bases Station and to Cluster Heads
packets_TO_BS=0;
packets_TO_CH=0;
%counter for bit transmitted to Bases Station and to Cluster Heads 
%per round
PACKETS_TO_CH(r+1)=0;
PACKETS_TO_BS(r+1)=0;

figure(1);

for i=1:1:n
    %checking if there is a dead node
    if (S(i).E<=0)
       plot(S(i).xd,S(i).yd,'^','LineWidth',1, 'MarkerEdgeColor','k', 'MarkerFaceColor','y', 'MarkerSize',8);
        dead=dead+1;
        if(S(i).ENERGY==1)
            dead_a=dead_a+1;
        end
        if(S(i).ENERGY==0)
            dead_n=dead_n+1;
        end
        hold on;    
    end
    if S(i).E>0
        S(i).type='N';
        if (S(i).ENERGY==0)  
     plot(S(i).xd,S(i).yd,'o','LineWidth',1, 'MarkerEdgeColor','k', 'MarkerFaceColor','g', 'MarkerSize',8);
        end
        if (S(i).ENERGY==1)  
        plot(S(i).xd,S(i).yd,'+','LineWidth',3, 'MarkerEdgeColor','k', 'MarkerFaceColor','r', 'MarkerSize',8);
        end
        hold on;
    end
end
plot(S(n+1).xd,S(n+1).yd,'x','LineWidth',1, 'MarkerEdgeColor','k', 'MarkerFaceColor','r', 'MarkerSize',8); 


STATISTICS(r+1).DEAD=dead;
DEAD(r+1)=dead;
DEAD_N(r+1)=dead_n;
DEAD_A(r+1)=dead_a;
%          plot(S(n+1).xd,S(n+1).yd,'o', 'MarkerSize', 12, 'MarkerFaceColor', 'r');
%          plot(S(n+1).xd,S(n+1).yd,'x','LineWidth',1, 'MarkerEdgeColor','k', 'MarkerFaceColor','r', 'MarkerSize',8);  
%When the first node dies
if (dead==1)
    if(flag_first_dead==0)
        first_dead=r;
        flag_first_dead=1;
    end
end

countCHs=0;
cluster=1;
for i=1:1:n
   if(S(i).E>0)
   temp_rand=rand;     
   if ( (S(i).G)<=0)

 %Election of Cluster Heads
 if(temp_rand<= (p/(1-p*mod(r,round(1/p)))))
            countCHs=countCHs+1;
            packets_TO_BS=packets_TO_BS+1;
            PACKETS_TO_BS(r+1)=packets_TO_BS;
            
            S(i).type='C';
            S(i).G=round(1/p)-1;
            C(cluster).xd=S(i).xd;
            C(cluster).yd=S(i).yd;
             plot(S(i).xd,S(i).yd,'k*');
            
            distance=sqrt( (S(i).xd-(S(n+1).xd) )^2 + (S(i).yd-(S(n+1).yd) )^2 );
            C(cluster).distance=distance;
            C(cluster).id=i;
            X(cluster)=S(i).xd;
            Y(cluster)=S(i).yd;
            cluster=cluster+1;
            
            %Calculation of Energy dissipated
            distance;
            if (distance>do)
                S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            %S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            end
            if (distance<=do)
                S(i).E=S(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            %S(i).E=S(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
            Energy_disp(r+1) =  S(i).E;
        end     
    
    end
  end 
end

STATISTICS(r+1).CLUSTERHEADS=cluster-1;
CLUSTERHS(r+1)=cluster-1;

%Election of Associated Cluster Head for Normal Nodes
for i=1:1:n
   if ( S(i).type=='N' && S(i).E>0 )
     if(cluster-1>=1)
       min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );
       min_dis_cluster=1;
       for c=1:1:cluster-1
           temp=min(min_dis,sqrt( (S(i).xd-C(c).xd)^2 + (S(i).yd-C(c).yd)^2 ) );
           if ( temp<min_dis )
               min_dis=temp;
               min_dis_cluster=c;
           end
       end
       
       %Energy dissipated by associated Cluster Head
            min_dis;
            if (min_dis>do)
                S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S(i).E=S(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
        %Energy dissipated
        if(min_dis>0)
         distance=sqrt( (S(C(min_dis_cluster).id).xd-(S(n+1).xd) )^2 + (S(C(min_dis_cluster).id).yd-(S(n+1).yd) )^2 );
          S(C(min_dis_cluster).id).E = S(C(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
if (distance>do)
                S(C(min_dis_cluster).id).E=S(C(min_dis_cluster).id).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            end
            if (distance<=do)
                S(C(min_dis_cluster).id).E=S(C(min_dis_cluster).id).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
          PACKETS_TO_CH(r+1)=n-dead-cluster+1; 
        end

       S(i).min_dis=min_dis;
       S(i).min_dis_cluster=min_dis_cluster;
           
   end
 end
end
hold on;

countCHs;
rcountCHs=rcountCHs+countCHs;
sum=0;
for i=1:1:n
if(S(i).E>0)
    sum=sum+S(i).E;
end
end
avg=sum/n;
STATISTICS(r+1).AVG=avg;
sum;


%Code for Voronoi Cells
%Unfortynately if there is a small
%number of cells, Matlab's voronoi
%procedure has some problems
warning('OFF');
[vx,vy]=voronoi(X(:),Y(:));
plot(X,Y,'g+',vx,vy,'m-');
hold on;
voronoi(X,Y);
axis([10 xm 0 ym]);

end
% figure1 = figure11;
% % Create axes
% axes1 = axes('Parent',figure1,'YGrid','on','XGrid','on','GridLineStyle','--');
% box(axes1,'on');
% hold(axes1,'all');
figure(2);
for r=0:1:24
    ylabel('Average Energy of Each Node');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).AVG STATISTICS(r+2).AVG],'red');
    hold on;
end
figure(3);
for r=0:1:49
    ylabel('Average Energy of Each Node');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).AVG STATISTICS(r+2).AVG],'red');
    hold on;
end
figure(4);
for r=0:1:74
    ylabel('Average Energy of Each Node');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).AVG STATISTICS(r+2).AVG],'red');
    hold on;
end
figure(5);
for r=0:1:99
    ylabel('Average Energy of Each Node');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).AVG STATISTICS(r+2).AVG],'red');
    hold on;
end
figure(6);
for r=0:1:24
    ylabel('Number of Dead Nodes');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).DEAD STATISTICS(r+2).DEAD],'red');
    hold on;
end
figure(7);
for r=0:1:49
        ylabel('Number of Dead Nodes');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).DEAD STATISTICS(r+2).DEAD],'red');
    hold on;
end
figure(8);
for r=0:1:74
        ylabel('Number of Dead Nodes');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).DEAD STATISTICS(r+2).DEAD],'red');
    hold on;
end
figure(9);
for r=0:1:99
        ylabel('Number of Dead Nodes');
    xlabel('Round Number');
    plot([r r+1],[STATISTICS(r+1).DEAD STATISTICS(r+2).DEAD],'red');
    hold on;
end