DESIGN OF ENERGY DETECTOR, MATCHED FILTER, COVARIANCE AND EIGEN VALUE SPECTRUM SENSING TECHNIQUES

clc
clear all;
close all;
warning off;

disp('----------------------------------');
disp('Greetings from Verilog Course Team');
disp('Website:www.verilogcourseteam.com');
disp('Email:info@verilogcourseteam.com');
disp('WhatsApp @ +91 790 456 8 456');
disp('----------------------------------');

disp('*************************************************************************************************');
disp('DESIGN OF ENERGY DETECTOR, MATCHED FILTER, COVARIANCE AND EIGEN VALUE SPECTRUM SENSING TECHNIQUES');
disp('*************************************************************************************************');
pause(5);
detail_info;
pause(5);

disp('---------------------------------');
disp('EXECUTING OR RULE WITH AWGN NOISE');
disp('---------------------------------');
srt_node=1;
N=20;
M=300;
hval=1;
ind=1;
snr_range=1:1:20;

RULE_ERROR
pause;
%%
disp('---------------------------------------');
disp('EXECUTING MAJORITY RULE WITH AWGN NOISE');
disp('---------------------------------------');
srt_node=N/2;
N=20;
M=300;
hval=1;
ind=1;
snr_range=1:1:20;

RULE_ERROR
pause;

%%
disp('----------------------------------');
disp('EXECUTING AND RULE WITH AWGN NOISE');
disp('----------------------------------');

srt_node=N;
N=20;
M=300;
hval=1;
ind=1;
snr_range=1:1:20;
RULE_ERROR

 

 

 % N=20;
% M=300;
% snr=10;


nsval=wgn(N,M,snr);
rjih0=nsval;
rjih1=ones(N,M)+nsval;

pfrange=0.1:0.1:1;

ind=1;
for pf=pfrange
    
    thres=gaminv(1-pf,1);
    yjh0=(sum((rjih0.').^2))/10000;
    yjh1=(sum((rjih1.').^2))/10000;
    
    
    pd=length(find(yjh1>=thres));
    
    pfxx=length(find(yjh0>thres));
    
    pdprop(ind)=pd/N;
    
    pmprop(ind)=1-pdprop(ind);
    
    
    pfprop(ind)=pf;
    
    
    qdprop_res=0;
    qrprop_res=0;
    
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=qrprop_res+qmprop_res;
    
    ind=ind+1;
    
end

outqm=sort(outqm,'descend');
energy_res_x=pfprop;
energy_res_y=outqm;

%%

nsval=wgn(N,M,snr);
rjih0=nsval;
rjih1=ones(N,M)+nsval;
ind=1;
for pf=pfrange
    
    thres=gaminv(1-pf,1);
    
    hmf=phased.MatchedFilter;
    
    for k3=1:N
        rjih0_filt(k3,:)=step(hmf,rjih0(k3,:));
        rjih1_filt(k3,:)=step(hmf,rjih1(k3,:));
    end
    yjh0=(sum((rjih0_filt.').^2))/10000;
    yjh1=(sum((rjih1_filt.').^2))/10000;
    pd=length(find(yjh1>=thres));
    pfxx=length(find(yjh0>thres));
    pdprop(ind)=pd/N;
    pmprop(ind)=1-pdprop(ind);
    pfprop(ind)=pf;
    qdprop_res=0;
    qrprop_res=0;
    
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=qrprop_res+qmprop_res;
    ind=ind+1;
end
outqm=sort(outqm,'descend');

match_res_y=outqm;
match_res_x=pfprop;


%%


nsval=wgn(N,M,snr);
rjih0=nsval;
rjih1=ones(N,M)+nsval;


smooth_val=4;
ind=1;
for pf=pfrange
    kcov1=1+((smooth_val-1)*sqrt(2/(M*3.14)));
    kcov2=1-(qfuncinv(pf)*sqrt(2/(M)));
    kcov=(kcov1/kcov2)*2;
    
    

    for k3=1:N
        rnval=rjih0(k3,:);
        [outlm tcov]=auto_corr_function(rnval,smooth_val);
        tcov_finalh0(k3)=(tcov);
        rnval=rjih1(k3,:);
        [outlm tcov]=auto_corr_function(rnval,smooth_val);
        tcov_finalh1(k3)=tcov;
    end
    
    yjh0=tcov_finalh0;
    yjh1=tcov_finalh1;
    pd=length(find(yjh1>=kcov));
    pfxx=length(find(yjh0>kcov));
    pdprop(ind)=pd/N;
    pmprop(ind)=1-pdprop(ind);
    pfprop(ind)=pf;
    qdprop_res=0;
    qrprop_res=0;
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=qrprop_res+qmprop_res;
    
    ind=ind+1;
    
end
outqm=sort(outqm,'descend');
cov_res_x=pfprop;
cov_res_y=outqm;

%%

outqm=zeros(1,length(pfrange));
smooth_val=4;
ind=1;
for pf=pfrange(1:end-2)
    
    keig1=((sqrt(M+15*smooth_val)).^2)/((sqrt(M-15*smooth_val)).^2);
    keig2=1+(((sqrt(M+15*smooth_val)).^(-2/3))/((M*15*smooth_val).^(1/6)))*expinv(1-pf);
    keig=(keig1*keig2)*150;
    
    for k3=1:N
        rnval=rjih0(k3,:);
        [outlm teig]=eig_val_function(rnval,smooth_val);
        tcov1_finalh0(k3)=teig;
        rnval=rjih1(k3,:);
        [outlm teig]=eig_val_function(rnval,smooth_val);
        tcov1_finalh1(k3)=teig;
        
        
    end
       
    yjh0=tcov1_finalh0;
    yjh1=tcov1_finalh1;
      
    pd=length(find(yjh1>=keig));
    pfxx=length(find(yjh0>keig));
    pdprop(ind)=pd/N;
    pmprop(ind)=1-pdprop(ind);
      
    pfprop(ind)=pf;
    
    
    qdprop_res=0;
    qrprop_res=0;
    
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=qrprop_res+qmprop_res;
    ind=ind+1;
    
end
eig_res_x=pfprop;
eig_res_y=outqm;

%%
figure,plot(energy_res_x,energy_res_y,'r:s','linewidth',2);
hold on,plot(match_res_x,match_res_y,'b-^','linewidth',2);
hold on,plot(cov_res_x,cov_res_y,'g->','linewidth',2);
hold on,plot(eig_res_x,eig_res_y,'k-<','linewidth',2);
grid on;
xlabel('Propability of False Alarm');
ylabel('Propability of Missed');
legend('Energy Detector','Matched Filter','Covariance','Eigen Method');

x1=normrnd(0,0.5);
x2=normrnd(0,0.5);
hval=sqrt(x1^2+x2^2);


for snr=snr_range
    

    nsval=wgn(N,M,snr);
    rjih0=nsval;
    rjih1=(hval*ones(N,M))+nsval;

    pf=0.5;
    thres=gaminv(1-pf,1);
    yjh0=(sum((rjih0.').^2))/10000;
    yjh1=(sum((rjih1.').^2))/10000;
    
    
    pd=length(find(yjh1>=thres));
    
    pfxx=length(find(yjh0>thres));
    pdprop(ind)=pd/N;
    pmprop(ind)=1-pdprop(ind);
    pfprop(ind)=pf;
    pfprop1(ind)=snr;
       
    
    qdprop_res=0;
    qrprop_res=0;
    
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=(qrprop_res+qmprop_res)-1;
    
    ind=ind+1;
    
end
error_prop=sort(error_prop,'descend');
energy_res_x=pfprop1;
energy_res_y=error_prop;

%% 2.    Matched filter

ind=1;
for snr=snr_range
    
    nsval=wgn(N,M,snr);
    rjih0=nsval;
    rjih1=(hval*ones(N,M))+nsval;
    
    
    thres=gaminv(1-pf,1);
    hmf=phased.MatchedFilter;
    
    for k3=1:N
        rjih0_filt(k3,:)=step(hmf,rjih0(k3,:));
        rjih1_filt(k3,:)=step(hmf,rjih1(k3,:));
    end
    
    yjh0=(sum((rjih0_filt.').^2))/10000;
    yjh1=(sum((rjih1_filt.').^2))/10000;
    pd=length(find(yjh1>=thres));
    pfxx=length(find(yjh0>thres));
    pdprop(ind)=pd/N;
    pmprop(ind)=1-pdprop(ind);
    pfprop(ind)=pf;
    pfprop1(ind)=snr;
    
    
    qdprop_res=0;
    qrprop_res=0;
    
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=(qrprop_res+qmprop_res)-1;
    ind=ind+1;
end
error_prop=sort(error_prop,'descend');
match_res_x=pfprop1;
match_res_y=error_prop;


%% 3.    Covariance

smooth_val=4;
ind=1;

for snr=snr_range
    
    nsval=wgn(N,M,snr);
    rjih0=nsval;
    rjih1=(hval*ones(N,M))+nsval;
    
    
    kcov1=1+((smooth_val-1)*sqrt(2/(M*3.14)));
    kcov2=1-(qfuncinv(pf)*sqrt(2/(M)));
    kcov=kcov1/kcov2;
     

    for k3=1:N
        rnval=rjih0(k3,:);
        [outlm tcov]=auto_corr_function(rnval,smooth_val);
        tcov_finalh0(k3)=(tcov);
        rnval=rjih1(k3,:);
        [outlm tcov]=auto_corr_function(rnval,smooth_val);
        tcov_finalh1(k3)=tcov;
    end
    
    yjh0=tcov_finalh0;
    yjh1=tcov_finalh1;
    pd=length(find(yjh1>=kcov));
    pfxx=length(find(yjh0>kcov));
    pdprop(ind)=pd/N;
    pmprop(ind)=1-pdprop(ind);
    
    pfprop(ind)=pf;
    pfprop1(ind)=snr;
    
    qdprop_res=0;
    qrprop_res=0;
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=(qrprop_res+qmprop_res)-1;
    
    ind=ind+1;
    
end
error_prop=sort(error_prop,'descend');
cov_res_x=pfprop1;
cov_res_y=error_prop;

%% eigen method

smooth_val=6;
ind=1;
for snr=snr_range
    
    nsval=wgn(N,M,snr);
    rjih0=nsval;
    rjih1=(hval*ones(N,M))+nsval;
    
    
    keig1=((sqrt(M+15*smooth_val)).^2)/((sqrt(M-15*smooth_val)).^2);
    keig2=1+(((sqrt(M+15*smooth_val)).^(-2/3))/((M*15*smooth_val).^(1/6)))*expinv(1-pf);
    keig=(keig1*keig2)*20;
    
    for k3=1:N
        rnval=rjih0(k3,:);
        [outlm teig]=eig_val_function(rnval,smooth_val);
        tcov1_finalh0(k3)=teig;
        rnval=rjih1(k3,:);
        [outlm teig]=eig_val_function(rnval,smooth_val);
        tcov1_finalh1(k3)=teig;
        
        
    end
    
    
    yjh0=tcov1_finalh0;
    yjh1=tcov1_finalh1;
        
    pd=length(find(yjh1>=keig));
    pfxx=length(find(yjh0>keig));
    pdprop(ind)=pd/N;
    pmprop(ind)=1-pdprop(ind);
    pfprop(ind)=pf;
    pfprop1(ind)=snr;
    
    qdprop_res=0;
    qrprop_res=0;
    pdprop_res=pdprop(ind);
    pfprop_res=pfprop(ind);
    
    for l=srt_node:1:N
        qdprop_res=qdprop_res+(factorial(N)*(pdprop_res^l)*((1-pdprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
        qrprop_res=qrprop_res+(factorial(N)*(pfprop_res^l)*((1-pfprop_res)^(N-l))/(factorial(l)*factorial(N-l)));
    end
    
    outqd(ind)=qdprop_res;
    qmprop_res=1-qdprop_res;
    outqm(ind)=qmprop_res;
    error_prop(ind)=(qrprop_res+qmprop_res)-1;
    
    ind=ind+1;
    
end
error_prop=sort(error_prop,'descend');
eig_res_x=pfprop1;
eig_res_y=error_prop;

%%

figure,plot(energy_res_x,energy_res_y,'m:s','linewidth',2);
hold on,plot(match_res_x,match_res_y,'g-^','linewidth',2);
hold on,plot(cov_res_x,cov_res_y,'c->','linewidth',2);
hold on,plot(eig_res_x,eig_res_y,'y-<','linewidth',2);
grid on;
xlabel('SNR');
ylabel('Total Error Propability ');
legend('Energy Detector','Matched Filter','Covariance','Eigen Method');