clc
close all
clear all
L = 10; % Number of sensing samples
iter =10^4; % Number of iterations for Monte Carlo Simulation
Pf =0.1:0.1:1; % Probability of False Alarm
snr_db = -1; % Signal-to-noise ratio (SNR) in dB
snr = 10.^(snr_db./10); % SNR in linear scale
for tt = 1:length(Pf) % Calculating threshold for each value of Pf
energy_fin = []; % Initialization
n = [];
y = [];
energy = [];
energy_fin = [];
n=randn(iter,L); % Gaussian noise, mean 0, variance 1
y = n; % Received signal at the secondary user under the hypothesis H0
energy = abs(y).^2; % Energy of received signal over L sensing samples under the hypothesis H0
for kk=1:iter % Start simulation loop to calculate the threhsold
energy_fin(kk,:) =sum(energy(kk,:)); % Test Statistic of the energy detection
end
energy_fin = energy_fin'; % Taking transpose to arrage values in descending order
energy_desc = sort(energy_fin,'descend')'; % Arrange values in descending order
thresh_sim(tt) = energy_desc(ceil(Pf(tt)*iter)); % Threshold obtained by simulations; the first 'Pf' fraction of values lie above the threshold
tt
end
%% Simulated probability of detection for awgn channel
for tt = 1:length(Pf)
s = []; % Initializtion
h = []; % Initializtion
mes = randi([0 1],iter,L); % Generating 0 and 1 with equal probability for BPSK
s = (2.*(mes)-1); % BPSK modulation
y = sqrt(snr).*s + n; % Received signal y at the secondary user
energy = (abs(y).^2); % Received energy under the hypotheis H1
for kk=1:iter % Number of Monte Carlo Simulations
energy_fin(kk) =sum(energy(kk,:)); % Test Statistic for the energy detection undet the hypothesis H1
end
upper = (energy_fin >= thresh_sim(tt)); % Checking whether the received energy above the threshold
i = sum(upper); % Count how many times out of 'iter', the received energy is above the threshold.
Pd_sim(tt) = i/kk; % Calculation of the probability of detection (simulated)
pm_sim(tt)=1-Pd_sim(tt);
tt
end
plot(Pf,pm_sim,'r')
hold on
%% Probability of detection (Theory):awgn
thresh_theory = 2*gammaincinv(1-Pf,L/2);% Theory threhsold, from the formula of Pf = gammainc(thresh_theory./(2), L/2,'upper')
u = L./2; % Time-Bandwidth product
for pp = 1:length(Pf)
n = 0:1:u-2;
Pd_theory_awgn = marcumq(sqrt(L*snr),sqrt(thresh_theory),L./2);
pm=1-Pd_theory_awgn;
end
plot(Pf,pm,'b*') % ROC curve
hold on
grid on
title(' Pf vs Pm for SNR = -1 dB in AWGN channel');
legend('Simulated P_d','Theory P_d')
xlabel('Probability of False Alarm')
ylabel('Probability of Missed Detection')
