clear all
close all;
load mydata noise1 SIGMA   
echo on
to=0.2; 
[R,r] = wireless_microphone(to); 
r=r(1:1000);
%S=40;
tal=3;
N=length(r);
q=20;
for snr=-10:1:-1
    Td(snr+11)=0
R2=0
for time=1:1:10
     R2=R2+covariance(noise1(snr+11,:),q,tal)
end
R2=R2/10
R2=inv(R2)

for m=1:1:1000
    noise_result(m)=0
    R1=estimate2(noise1(snr+11,:),q,tal)
    noise_result(m)=R1*R2*R1'
end
noise_result=sort(noise_result)
door=noise_result(990)


    for m=1:1:100   
          signal(snr+11,:)=sqrt(2*10^(snr/10))*r+noise1(snr+11,:)
          R1=estimate2(signal(snr+11,:),q,tal)
%           R2=covariance(signal,20,3)
%           R2=inv(R2)
          signal_result=R1*R2*R1'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          
       Pfa=0.1;                                                     % 虛警概率可改：0.3,0.2,0.1，0.05
       aa=sqrt(2*log10(1./Pfa)); 
       th(snr+11)=SIGMA(snr+11)*aa;  
       threshold=chi2inv(0.9,q);            
               % snr_avg(snr)=(sum(rr(snr,:).^2))/(sum(noise1(snr,:).^2));
                snr_avg(snr+11)=((signal(snr+11,:).^2))/((noise1(snr+11,:).^2));
                FT=0.5;                                             % 加權指數因子FT可改為0.05,0.6,0.8,1.2,3
                th4(snr+11)=threshold*(sqrt(snr_avg(snr+11))/ signal_result).^FT;
                
                K=0.5;                                              % 乘性因子K可改為 0.02,0.03,0.05
                th6(snr+11)=K*SIGMA(snr+11);     
          
          
          T(m)=0
          if signal_result>door
              T(m)=1
          end
           D(m)=0
          if signal_result>threshold
              D(m)=1
          end
    end
     Td(snr+11)=sum(T)
     Dd(snr+11)=sum(D)
end
figure(1)
 k=1:1:10;
 th4=sort(th4);
 plot(k,threshold,'-ko',k,th4(k),'--r*',k,th6(k),'b+-',k,th(k),'-cd',k,door,'-bo')
legend('(1)固定閾值1','(2)指數加權控制閾值','(3)自適應恒虛警檢測閾值','(4)自適應門限加權控制閾值','(5)固定閾值2')
 grid on ;  
xlabel('SNR/db')                
ylabel('幅度/ V')  
title('門限閾值比較圖');
 figure(2)
 plot(k,door,'-bo')
 figure(3)
 plot( noise_result,'-ro')
 figure(4)
 plot( signal_result,'-k+')

figure(5)
plot(Td/100,'b+-')
xlabel('SNR/db')                
ylabel('Pd')
figure(6)
plot(Dd/100,'k+-')
xlabel('SNR/db')                
ylabel('Pd')