SHARE
TWEET

4_Energy_Detection_Elyes (plage_N_Pd_vs_Pfa)

famansour May 23rd, 2019 75 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. clc;
  2. clear all;
  3. close all;
  4. %% Signal generation
  5. t=[0:1/(2*pi*400):1]; % intervalle de temps pour la construction du signal.
  6. x=sin(50*2*pi*t); %signal message.
  7. y = ssbmod(x,12e+7,48e+7,pi/4); % modulation AM du message x.
  8.  
  9. SNRv= -5; % plage de valeur de SNR pour tracer les courbes ROC et comparer les perfs suivant le SNR.
  10. [y1, noise] = add_awgn_noise(y,SNRv); % ajout du bruit au signal modulé.
  11. mesur_sig=[noise noise y1 y1]; % construction du signal de transmission.
  12. for i=1:length(mesur_sig(:,1))
  13.     mesur_sig(i,:)=mesur_sig(i,:)-mean(mesur_sig(i,:));
  14. end
  15. siglabel=zeros(length(mesur_sig),1); % construction du siglabel (0 quand noise et 1 quand signal utile)
  16. siglabel(end-2*length(y1)+1:end)=1;
  17.  
  18. Pf = 0.0001:0.1:0.9001; % de meme que pour le SNR
  19.  
  20.  
  21. [Tn,PU,N]=Energy(mesur_sig,Pf,noise,0); % appel de la fonction qui calcule l'energie du signal recu.
  22.  
  23.  
  24. ROC(siglabel, PU, SNRv,Pf, Tn, N) % tracage des courbes de perfs.
  25.  
  26. %%%%%%%%%%%%%%%%%%%%%%      Start Pasting      %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  27. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  28. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  29. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  30. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  31. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  32. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  33.  
  34. function [y, noise] = add_awgn_noise(x,SNR_dB)
  35.  
  36. %Function to add AWGN to a given signal
  37.  
  38. %Authored by Mathuranathan Viswanathan
  39.  
  40. %How to generate AWGN noise in Matlab/Octave by Mathuranathan Viswanathan
  41.  
  42. %is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0  International License.
  43.  
  44. %You must credit the author in your work if you remix, tweak, and build upon the work below
  45.  
  46.     %y=awgn_noise(x,SNR) adds AWGN noise vector to signal 'x' to generate a
  47.     %resulting signal vector y of specified SNR in dB
  48.     rng('default');%set the random generator seed to default (for comparison only)
  49.     L=length(x);
  50.     SNR = 10.^(SNR_dB/10); %SNR to linear scale
  51.     Esym=sum(abs(x).^2)/(L); %Calculate actual symbol energy
  52.     N0=Esym./SNR; %Find the noise spectral density
  53.     if(isreal(x)),
  54.         noiseSigma = sqrt(N0);%Standard deviation for AWGN Noise when x is real
  55.         noise = noiseSigma.'.*randn(1,L);%computed noise
  56.     else
  57.         noiseSigma=sqrt(N0/2);%Standard deviation for AWGN Noise when x is complex
  58.         noise = noiseSigma.'.*(randn(1,L)+1i*randn(1,L));%computed noise
  59.     end    
  60.     y = x + noise; %received signal    
  61. end
  62. function [Tn,PU,N]=Energy(sig,Pf,noise,gr)
  63.  
  64. %% Filtre passe bande Butterworth
  65. [A,B,C,D]=butter(2,[(12e+7-50) (12e+7+50)]/(48e+7/2));
  66. sos = ss2sos(A,B,C,D);
  67. fltpb = sosfilt(sos,sig);
  68.  
  69.  
  70. %% Energie du signal recu
  71. N = [50 500 1000 1500 2000 2500 3000 3500];
  72. for m=1:length(N)
  73.     tic
  74.     for i= N(m)+1:length(fltpb)
  75.         Tn(m,i) = (1/N(m))*sum(fltpb(i-N(m):i).^2); %calcul de l'energie du signal recu apres filtrage.
  76.     end
  77.     Tn(m,1:N(m)) = Tn(m,N(m)+1:2*N(m)); % on fixe les N premieres valeurs d'energie a cause du decalage introduit par le calcul de l'energie sur une fenetre.
  78.    
  79.     %TN(m,:) = Tn(m,:)./max(Tn(m,:));
  80.     timenrg(m)=toc;
  81. end
  82.  
  83. %% determination du seuil(valeur theorique pour un environnement gaussien)
  84.  
  85. for m=1:length(N)
  86.     n=length(Pf);
  87.     %determination seuil en fct du Pf.
  88.     for i=1:length(Tn)
  89.     thresh(m*n-n+1:m*n,i) = (sqrt(-2*var(noise).*log(Pf'))./sqrt(N(m)))+1; % reference, Detection en environnement non gaussien, Emmanuelle Jay.
  90.     end
  91.     %(qfuncinv(Pf(m))./sqrt(N))+ 1; % Valeur theorique du seuil, reference, Sensing Throughput Tradeoff in Cognitive Radio, Y. C. Liang
  92. end
  93.  
  94. for k = 1:length(N)
  95.     tic
  96.     for m=1:length(Pf)        
  97.         for i=1:length(thresh)
  98.             if(Tn(k,i) >= thresh(m+length(Pf)*(k-1),i)) %comparaison par rapport au seuil
  99.                 % PU present
  100.                 PU(m+length(Pf)*(k-1),i) = 1;
  101.             else
  102.                 % PU absent
  103.                 PU(m+length(Pf)*(k-1),i) = 0;
  104.             end
  105.         end
  106.     end
  107.     time(k)=toc;
  108. end
  109.  
  110.  
  111. %% Courbes
  112. if gr
  113. figure('Name','seuil de detection');
  114. plot(1:length(Tn),Tn/max(Tn),1:length(Tn),thresh)
  115. grid on;
  116. figure,
  117. plot(1:length(Tn),Tn/max(Tn),1:length(PU),PU)
  118. grid on;
  119. end
  120. figure(1),
  121. plot(N,time)
  122. xlabel('N (nombre echantillon)'),
  123. ylabel('T (temps)'),
  124. % figure(2),
  125. % plot(N,timenrg)
  126.  
  127. end
  128. function ROC(siglabel, PU, SNRv,Pf, Tn, N)
  129. %% ROC SNR vs Pd
  130.        
  131. % for k =1:length(PU(:,1))
  132. %     bd=0;
  133. %     md=0;
  134. %     for kk = 1:length(PU)
  135. %
  136. %         if ( siglabel(kk)==PU(k,kk)) % calcul de la proba de bonne detection (Pd) et de fausse detection (Pfa).
  137. %             bd=bd+1;
  138. %         else
  139. %             md=md+1;
  140. %         end
  141. %     end
  142. %
  143. %     Pd(k,:)=bd/length(PU);
  144. %     Pfa(k,:)=md/length(PU);
  145. %
  146. % end
  147. %
  148. % for m = 1:length(SNRv)
  149. %     n=length(Pf);
  150. %     Pdd(m,:)=Pd(m*n-n+1:m*n,:)';        
  151. % end
  152. % %Pdd=sort(Pdd);
  153. % Pdd=Pdd';
  154. %
  155. % [Pd_th,SNR_dB_th] = rocpfa(Pf); % calcul de Pd et SNR theorique suivant le Pfa que j'ai fixé.
  156. %
  157. % figure('Name', 'Pd vs SNR');
  158. % for k=1:length(Pf)
  159. %     Color_val=rand(length(Pf),3);
  160. %     hold on,
  161. %     plot(SNR_dB_th,Pd_th(:,k))
  162. %     macourbe1(k)=plot(SNRv,'O');
  163. %     set( macourbe1(k),'color',  Color_val(k,:), 'LineStyle', '-');
  164. %     text(SNRv,Pdd(k), ['Pfa =',num2str(Pf(k)),' '], 'color', Color_val(k,:), 'fontsize', 8)
  165. %     xlabel('SNR');ylabel('P_d');grid on;
  166. % end
  167.  
  168.  
  169. %% ROC Pd vs Pfa
  170.  
  171. figure('Name', 'Pd vs Pfa'),
  172. plot([0 1], [0 1], '--k');
  173.  
  174. for ki=1:length(N)
  175.  
  176.     Color_val=rand(length(N),3); % variable de couleur pour tracer chaque courbe d'une couleur differente.
  177.  
  178.     [tpri,fpri,thresholdsi] = roc(siglabel',Tn(ki,:)); %calcul de Pd et Pfa theorique.
  179.  
  180. hold on;
  181.  
  182. macourbe1(ki)=plot(fpri,tpri,'LineWidth',1.3);
  183. set( macourbe1(ki),'color',  Color_val(ki,:), 'LineStyle', '-');
  184. text(fpri(fix(length(fpri)/2)),tpri(fix(length(tpri)/2)), ['N =',num2str(N(ki)),' '], 'color', Color_val(ki,:), 'fontsize', 8)
  185. legend('limit theorical','N = 50','N = 500','N = 1000','N = 1500','N = 2000', 'N = 2500','N = 3000','N = 3500')
  186. title('ROC curve', 'fontsize', 9);
  187. xlabel('Probability of false alarm (Pfa)', 'fontsize', 11); ylabel('Probability of detection (Pd)', 'fontsize', 11);
  188. grid on;
  189. hold off;
  190.  
  191.  
  192.  
  193.  
  194.  
  195. end
  196.  
  197. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top