Matlab：繪制簡單能量的接收機工作特性曲線(Energy_detection_simulation_ok)

輸出結果

實作代碼

% This code is to plot receiver operating characteristic curve for simple energy

%繪制簡單能量的接收機工作特性曲線

% detection, when the primary signal is real Gaussian signal and noise is

% addive white real Gaussian. Here, the threshold is available

% analytically.

% Code written by: Sanket Kalamkar, Indian Institute of Technology Kanpur,

% India.

%% 以下代碼繪制在虛警機率一定時，檢測機率和信噪比之間的關系曲線稱為檢測器的檢測性能曲線

clc

close all

clear all

L = 1000;                                             % The number of samples

snr = 0.01:0.01:10;

Pf = 10e-4;                                           % Pf = Probability of False Alarm 虛警機率确定

%% Simulation to plot SNR vs.Probability of Detection (Pd)

for m = 1:length(snr)

   i = 0;

   thresh = (qfuncinv(Pf)./sqrt(L))+ 1;              % Theoretical value of Threshold, refer, Sensing-Throughput Tradeoff for Cognitive Radio Networks, Y. C. Liang

   for kk = 1:5000                                   % Number of Monte Carlo Simulations(https://cn.mathworks.com/discovery/monte-carlo-simulation.html)

       n = randn(1,L);                               % AWGN noise with mean 0 and variance 1

       s = sqrt(snr(m)).*randn(1,L);                 % Real valued Gaussina Primary User Signal

       y = s + n;                                    % Received signal at SU(認知使用者接收到的信号)

       energy = abs(y).^2;                           % Energy of received signal over N samples

       energy_fin =(1/L).*sum(energy);               % Test Statistic for the energy detection

       if(energy_fin >= thresh)                      % Check whether the received energy is greater than threshold, if so, increment Pd (Probability of detection) counter by 1

           i = i+1;

       end

   end

   Pd(m) = i/kk;

end

plot(10*log(snr), Pd,  'r')

xlabel('信噪比，機關db');

ylabel('檢測機率');

title('能量感覺檢測性能曲線');

grid on

hold on

%% Theroretical expression of Probability of Detection; refer above reference.

thresh = (qfuncinv(Pf)./sqrt(L))+ 1;

%Pd_the = qfunc(((thresh - (snr + 1)).*sqrt(L))./(sqrt(2).*(snr + 1))); % 原來代碼中的表達與論文中不一緻

for k = 1:length(snr)

   Pd_the(k) = qfunc(((thresh - (snr(k) + 1)).*sqrt(L))./(sqrt(2).*snr(k) + 1)); % 與論文中的方程式保持一緻

plot(10*log(snr), Pd_the, 'b')

legend('實際檢測機率', '理論檢測機率', 'Location', 'SouthEast');