• API
• FAQ
• Tools
• Archive
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.
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.

Top