From: Quan on
hello,every one!
As a novice of cognitive radio, I am now very interesting in the topic of energy detection (ED), especially the performance of ED over different fading channels.
Accoring to some papers, like On the energy detection of unknown signals over fading channels (which was included in IEEE Xplore, 2003). I have done some simulation about this in Matlab 2009b: Firstly, produce a real signal(single tone or BPSK signal), and then detect them after transmission over different channels with ED theory. However, the problem is: the complementary ROC curve (Pm-Pf) obtained by this simulation is very different from the theory results and analysis given in your papers. I have also referred some other high cited papers but their results are the same with yours, so I am very eager to wonder where the problem is? It seems that I have some misunderstanding of the detection process.
I haven’t found any solution for this problem after a long-time serious check. So I’d like to resort to your help. Part of my Mathlab script is attached as follows. Did you ever do such a simulation? Or If you have any program codes related to ED Simulation, could you please send me some for study?

Matlab Simulation Code: (Is there any problem ??)
close all;
clear all;
clc;

m=5;
N =2*m; %sample points N=T/Ts=T/1/2W=2TW
Base= 0.01:0.02:1;
Pf =Base.^2; %False alarm
snr_avgdB =5; %SNR=Ps/(N0*W)
snr_avg = power(10,snr_avgdB/10); %db to linear snr
for i=1:length(Pf)
Over_Num = 0;
for kk = 1:1000 %Monter-Carlo times
%If Transimmting Single tone
t = 1:N; %samples
x = sin(pi*t/2); %single tone, Fs=2F0
amp = sqrt(1/2/snr_avg); %amp. for noise, avg power of signal=1/2
noise = amp*randn(1,N); %Noise production
pn = (std(noise)).^2; %Noise power average

% %If trans. BPSK singnal
% Ns=10;%bits Number
% M=2;
% Fd=1;%code rate
% Fs=1;%sample freq.
% a=randint(1,Ns,M); %source
% x=dmodce(a,Fd,Fs,'psk',M); %BPSK baseband singal
% amp = sqrt(1/snr_avg); %amp. for noise,, avg power of signal=1
% noise = amp*randn(1,N); %Noise production
% pn = (std(noise)).^2; %Noise power average

signal = x(1:N);
ps = mean(abs(signal).^2); %signal average power

%awgn channels
Rev_sig = signal + noise; %received signal for detection

Th(i) = gaminv(1-Pf(i),m,1)*2; %Threshold for given Pf(i) according to Digham 2003
accum_power(i) = sum(abs(Rev_sig.^2))/pn; %accumulated received power(normalized to noise power)

if accum_power(i)>Th(i)
Over_Num = Over_Num +1; %decide 1 or 0, present or absent
end
end
Pd_sim(i) = Over_Num / kk; %Detection probability computation
Pd_theory(i) = marcumq(sqrt(snr_avg*2*m),sqrt(Th(i)),m);
%Theory results according Digham 2003 in awgn channel
Pd_appm(i) = 0.5*erfc((sqrt(2)*erfcinv(2*Pf(i))-snr_avg*sqrt(m))/sqrt(2+4*snr_avg));
%approximation when m>100, chi-squre appoximate gaussion
end

Pm_sim =1-Pd_sim;
Pm_theory = 1- Pd_theory;
Pm_appm =1-Pd_appm;

figure
loglog(Pf,Pm_sim,'*r',Pf,Pm_theory,'-.k',Pf,Pm_appm,'--b');%Complementary ROC Curve

title('Complementary ROC of ED under AWGN')
grid on
axis([0.0001,1,0.0001,1]);
xlabel('Pf');
ylabel('Pm');
legend('Simulation','Theory','Arrpoximation');

From: Quan on
sorry,I made some mistakes。The following is the modified,but the results are still not correct。
who can tell me the reason??thanks in advance!
close all;
clear all;
clc;

m=1000;
N =2*m; %sample points N=T/Ts=T/1/2W=2TW
Base= 0.01:0.02:1;
Pf =Base.^2; %False alarm
snr_avgdB =-20; %SNR=Ps/(N0*W)
snr_avg = power(10,snr_avgdB/10); %db to linear snr
for i=1:length(Pf)
Over_Num = 0;
Th(i) = chi2inv(1-Pf(i),2*m);
%Th(i) = gaminv(1-Pf(i),m,1)*2; %Threshold for given Pf(i)according to Digham 2003

for kk = 1:1000 %Monter-Carlo times
%If Transimmting Single tone
t = 1:N; %samples
x = sin(pi*t); %single tone, Fs=2F0
noise = randn(1,N); %Noise production
pn = mean(noise.^2); %Noise power average
% pn = (std(noise))^2; %Noise power average
amp = sqrt(noise.^2*snr_avg);
x=amp.*x./abs(x);
SNRdB_Sample=10*log10(x.^2./(noise.^2));

% %If trans. BPSK singnal
% Ns=10;%bits Number
% M=2;
% Fd=1;%code rate
% Fs=1;%sample freq.
% a=randint(1,Ns,M); %source
% x=dmodce(a,Fd,Fs,'psk',M); %BPSK baseband singal
% amp = sqrt(1/snr_avg); %amp. for noise,, avg power of signal=1
% noise = amp*randn(1,N); %Noise production
% pn = (std(noise)).^2; %Noise power average

signal = x(1:N);
ps = mean(abs(signal).^2); %signal average power

%awgn channels
Rev_sig = signal + noise; %received signal for detection
accum_power(i) = sum(abs(Rev_sig.^2))/pn; %accumulated received power(normalized to noise power)
if accum_power(i)>Th(i)
Over_Num = Over_Num +1; %decide 1 or 0, present or absent
end
end
Pd_sim(i) = Over_Num / kk; %Detection probability computation
Pd_theory(i) = marcumq(sqrt(snr_avg*2*m),sqrt(Th(i)),m);
%Theory results according Digham 2003 in awgn channel
Pd_appm(i) = 0.5*erfc((sqrt(2)*erfcinv(2*Pf(i))-snr_avg*sqrt(m))/sqrt(2+4*snr_avg));
%approximation when m>100, chi-squre appoximate gaussion
end

Pm_sim =1-Pd_sim;
Pm_theory = 1- Pd_theory;
Pm_appm =1-Pd_appm;

figure
loglog(Pf,Pm_sim,'*r',Pf,Pm_theory,'-.k',Pf,Pm_appm,'--b');%Complementary ROC Curve

title('Complementary ROC of ED under AWGN')
grid on
axis([0.0001,1,0.0001,1]);
xlabel('Pf');
ylabel('Pm');
legend('Simulation','Theory','Arrpoximation');
From: RAV R on
Hi.. could you get the solution to your problem?