there are some issues in this code. First of all, you should remove the mean from your signal. Second, using the toeplitz assumption, you also assume that your signal is stationary. Third, I would use SVD instead of eig for a lot of reason (if you don't want to use the toepliz  assumption anymore, your covariance matrix could be not a full rank matrix, and svd gives you the closest subspace approximation for your covariance matrix representation; furthermore, svd return the absolute values of the eigenvalues).

Good luck with your research,

Francesco Onorati

@ Francesco

Thanks for answering my concern, I think I still have some issue in understanding (L) which is the consecutive samples; Could you please highlight if it is implemented in the right way based on the references [1] and [2]? 

After your corrections the code becomes as follows:

L=4; % consecutive samples

y=randn(1,L);

y=y-mean(y);

r = xcorr(y);

N=length(r);

R=toeplitz(r);

[EG,S,V] = svd(R);

MX=max(EG(:));

MN=min(EG(:));

TH=((sqrt(N)+sqrt(L))/(sqrt(N)-sqrt(L)))^2*(1+((sqrt(N)+sqrt(L))^-2/3/(N*L)^1/6)*0.45);

Best regards 

Francesca is right 

 I have tried to simulate MME algorithm by following the procedure as you guys discussed in the above section. But it gives Pd always 1. Additionally, How to plot ROC curve for Pd Vs SNR. Please help me in this regard, Your help will highly be appreciated. Here is code for MME algorithm:

clc

close all

clear all

N = 60000; %# of samples

L=4; %Length of recieved signal vector

snr_dB = -10; % SNR in decibels

snr = 10.^(snr_dB./10); % Linear Value of SNR

Pf = 0.1:0.1:1; % Pf = Probability of False Alarm

%% Simulation to plot Probability of Detection (Pd) vs. Probability of False Alarm (Pf)

for m = 1:length(Pf)

m

i = 0;

for kk=1:100 % Number of Monte Carlo Simulations

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

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

y = s+n; % Received signal at SU

r = xcorr(y); % Auto-correlation of y

R=abs(toeplitz(r(1:L)))/N; % sample covariance matrix

[EG,S,V] = svd(R); %S computes eigenvalues of R

MX=max(S(:)); %Maximum eigenvalue

MN=min(diag(S)); %Minimum eigenvalue

TH(m)=((sqrt(N)+sqrt(L))/(sqrt(N)-sqrt(L)))^2*(1+((sqrt(N)+sqrt(L))^-2/3/(N*L)^1/6)*0.45);

test(m)=MN*TH(m);

if (MX>test(m))

display('Signal is present (As max_eigenvalue>Threshold*min_eigenvalue)....H1')

% Checks 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(Pf, Pd,'linewidth',2)

hold on

grid on

Francesco Onorati Please check the code given above, I 'm not geeting correct  answer

I am also getting the  same problem as discussed in my above post. Can anyone help me out?

Thanks!

@Muhammad Karam Shehzad 

I think there might be problem with your threshold. It is constant and not changing as per your Pf. The MX is always greater than test(m). Therefore you are getting Pd as always 1.

I have implemented MAXIMUM EIGENVALUE based detection for white noise signal... this might help. The plot is for comparison of energy detection with MED. reference used is "Maximum Eigen Value Detection: Theory and Application: Theory and Application" by Yonghong Zeng. This code will plot Fig. 4 in this reference. I am facing problem with generating wireless microphone signal which they have used as primary signal. if anyone having any idea please reply.

If you find any error please reply. 

clc;

clear;

close all;

tic;

%% SNR vs pd using maximum eigenvalue detection

Ns=100000;

SNR=-25:-16;

snr=10.^(SNR./10);

L=8;

pf=0.1;

F1_inv=0.45;

num_iter=2000;

a=((sqrt(Ns)+sqrt(L))^2)/(Ns);

b=1+((sqrt(Ns)+sqrt(L))^(-2/3)*F1_inv)/((Ns*L)^(1/6));

threshold=a*b;

%

w=waitbar(0,'calculating probability of false alarm');

l=1;

for i=1:length(SNR)

count=0;

for h=1:num_iter

signal=randn(1,Ns);

noise=randn(1,Ns);

noise_power=norm(noise)^2;

signal_power=norm(signal)^2;

mult=sqrt(snr(i)*noise_power/signal_power);

signal=mult*signal;

signal=signal+noise;

%noise=noise/std(noise);

k=0;

Cx=zeros(1,L);

for n=0:L-1

for j=1:1:Ns-L-1

Cx(n+1)=Cx(n+1)+signal(j)*signal(j+k);

end

k=k+1;

end

Cx=Cx/Ns;

Cx_mtx=toeplitz(Cx);

eig_value=eig(Cx_mtx);

max_eig=max(max(eig_value));

if max_eig>threshold

count=count+1;

end

end

pd(i)=count/2000;

waitbar(i/length(SNR));

end

close(w);

plot(SNR,pd,'-^');

toc;

%% SNR vs pd using energy detection

M=Ns/2;

Pfa=0.1;

threshold=gammaincinv(1-Pfa,M)*2;

w=waitbar(0,'SNR vs Pd using energy detection');

for i=1:length(SNR)

count=0;

for h=1:num_iter

signal=randn(1,Ns);

noise=randn(1,Ns);

signal_power=norm(signal)^2;

noise_power=norm(noise)^2;

mult=sqrt(snr(i)*noise_power/signal_power);

signal=mult*signal;

signal=signal+noise;

energy=norm(signal)^2;

if energy>threshold

count=count+1;

end

end

Pd1(i)=count/num_iter;

waitbar(i/length(SNR));

end

close(w);

hold on;

plot(SNR,Pd1,'-s');

legend('MED','ED');

xlabel('SNR');

ylabel('P_d');

@Kamal I am facing same problem for wirless microphone signals. Have you implemented that successfully?