一维ca-cfar(均值类恒虚警检测)与二维ca-cfar算法matlab实现

东方耀

一维ca-cfar(均值类恒虚警检测)与二维ca-cfar算法matlab实现

1. % Implement 1D CFAR using lagging cells on the given noise and target scenario.
   
2. 
   
3. % Close and delete all currently open figures
   
4. close all;
   
5. 
   
6. % Data_points
   
7. Ns = 1000;
   
8. 
   
9. % Generate random noise
   
10. s=abs(randn(Ns,1));
    
11. 
    
12. %Targets location. Assigning bin 100, 200, 300 and 700 as Targets with the amplitudes of 8, 9, 4, 11.
    
13. s([100 ,200, 350, 700])=[8 15 7 13];
    
14. 
    
15. %plot the output
    
16. figure ('Name','original_signal'),plot(s);
    
17. 
    
18. % TODO: Apply CFAR to detect the targets by filtering the noise.
    
19. 
    
20. % 1. Define the following:
    
21. % 1a. Training Cells
    
22. T = 12;
    
23. % 1b. Guard Cells
    
24. G = 4;
    
25. 
    
26. % Offset : Adding room above noise threshold for desired SNR
    
27. offset=5;
    
28. 
    
29. % Vector to hold threshold values
    
30. threshold_cfar = [];  % 存储阈值的数组
    
31. 
    
32. %Vector to hold final signal after thresholding
    
33. signal_cfar = [];     % 是否有目标？
    
34. 
    
35. % 2. Slide window across the signal length
    
36. for i = 1:(Ns-(G+T+1))     
    
37. 
    
38.     % 2. - 5. Determine the noise threshold by measuring it within the training cells
    
39.    
    
40.     noise_level =sum(s(i:i+T-1));
    
41.    
    
42.     % 6. Measuring the signal within the CUT
    
43.    
    
44.     threshold = (noise_level/T)*offset;   % 这里是乘以offset
    
45.     threshold_cfar=[threshold_cfar,{threshold}];
    
46.    
    
47.     signal=s(i+T+G);
    
48.    
    
49.     % 8. Filter the signal above the threshold
    
50.     if (signal<threshold)
    
51.         signal=0;   % 不是目标
    
52.     end
    
53.     signal_cfar = [signal_cfar, {signal}];
    
54. end
    
55. 
    
56. 
    
57. 
    
58. 
    
59. % plot the filtered signal
    
60. % signal_cfar2 = cell2mat(signal_cfar);
    
61. figure,plot (cell2mat(signal_cfar),'g--');
    
62. 
    
63. % plot original sig, threshold and filtered signal within the same figure.
    
64. figure,plot(s);
    
65. hold on,plot(cell2mat(circshift(threshold_cfar,G)),'r--','LineWidth',2)
    
66. hold on, plot (cell2mat(circshift(signal_cfar,(T+G))),'g--','LineWidth',4);
    
67. legend('Ori_Signal','CFAR Threshold','detection')

复制代码

1. clear all
   
2. close all;
   
3. clc;
   
4. 
   
5. %% Radar Specifications
   
6. %%%%%%%%%%%%%%%%%%%%%%%%%%%
   
7. % Frequency of operation = 77GHz
   
8. % Max Range = 200m
   
9. % Range Resolution = 1 m
   
10. % Max Velocity = 100 m/s
    
11. %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
12. 
    
13. max_range=200;
    
14. c = 3e8;
    
15. range_resolution= 1;
    
16. 
    
17. %Operating carrier frequency of Radar
    
18. fc= 77e9;             %carrier freq
    
19. wavelength = c/fc;
    
20. 
    
21. %% User Defined Range and Velocity of target
    
22. % *%TODO* :
    
23. % define the target's initial position and velocity. Note : Velocity
    
24. % remains contant
    
25. target_pos=100;
    
26. target_speed=30;
    
27. 
    
28. 
    
29. %% FMCW Waveform Generation
    
30. 
    
31. % *%TODO* :
    
32. %Design the FMCW waveform by giving the specs of each of its parameters.
    
33. % Calculate the Bandwidth (B), Chirp Time (Tchirp) and Slope (slope) of the FMCW
    
34. % chirp using the requirements above.
    
35. 
    
36. B_sweep = c/(2*range_resolution); %Calculate the Bandwidth (B)
    
37. 
    
38. T_chirp = 5.5*2*max_range/c;   % 为何啊？
    
39. 
    
40. slope=B_sweep/T_chirp;
    
41. 
    
42. 
    
43. 
    
44.                                                          
    
45. %The number of chirps in one sequence. Its ideal to have 2^ value for the ease of running the FFT
    
46. %for Doppler Estimation.
    
47. Nd=128;                   % #of doppler cells OR #of sent periods % number of chirps
    
48. 
    
49. %The number of samples on each chirp.
    
50. Nr=1024;                  %for length of time OR # of range cells
    
51. 
    
52. % Timestamp for running the displacement scenario for every sample on each
    
53. % chirp
    
54. t=linspace(0,Nd*T_chirp,Nr*Nd); %total time for samples
    
55. 
    
56. 
    
57. %Creating the vectors for Tx, Rx and Mix based on the total samples input.
    
58. Tx=zeros(1,length(t)); %transmitted signal
    
59. Rx=zeros(1,length(t)); %received signal
    
60. Mix = zeros(1,length(t)); %beat signal
    
61. 
    
62. %Similar vectors for range_covered and time delay.
    
63. r_t=zeros(1,length(t));  % 每时每刻的距离
    
64. td=zeros(1,length(t));   % 距离算出来的时延
    
65. 
    
66. 
    
67. %% Signal generation and Moving Target simulation
    
68. % Running the radar scenario over the time.
    
69. 
    
70. doppler_resolution = 1 / (Nd * T_chirp);
    
71. 
    
72. v_resolution = (wavelength*doppler_resolution)/2;  % 速度分辨率为2.07 为何？
    
73. 
    
74. for i=1:length(t)         
    
75.    
    
76.    
    
77.     % *%TODO* :
    
78.     %For each time stamp update the Range of the Target for constant velocity.
    
79.    
    
80.     r_t(i) = target_pos+(target_speed*t(i));
    
81.     td(i) = 2*r_t(i)/c; % Time delay
    
82.    
    
83.     % *%TODO* :
    
84.     %For each time sample we need update the transmitted and
    
85.     %received signal.
    
86.     Tx(i) = cos(2 * pi * (fc * t(i) + slope * (t(i)^2)/2));
    
87.     Rx(i) = cos(2 * pi * (fc * (t(i) - td(i)) + slope * ((t(i)-td(i))^2)/2));
    
88.    
    
89.     % *%TODO* :
    
90.     %Now by mixing the Transmit and Receive generate the beat signal
    
91.     %This is done by element wise matrix multiplication of Transmit and
    
92.     %Receiver Signal
    
93.     Mix(i) = Tx(i) .* Rx(i);    % 不是发射信号的复共轭 不是复信号啊
    
94.    
    
95. end
    
96. 
    
97. %% RANGE MEASUREMENT
    
98. 
    
99. 
    
100. % *%TODO* :
     
101. %reshape the vector into Nr*Nd array. Nr and Nd here would also define the size of
     
102. %Range and Doppler FFT respectively.
     
103. 
     
104. Mix_reshape = reshape(Mix,[Nr,Nd]);
     
105. 
     
106. % *%TODO* :
     
107. %run the FFT on the beat signal along the range bins dimension (Nr) and
     
108. %normalize.
     
109. 
     
110. sig_fft1_0=fft(Mix_reshape,Nr)./Nr;
     
111. 
     
112. % *%TODO* :
     
113. % Take the absolute value of FFT output
     
114. 
     
115. sig_fft1_1=abs(sig_fft1_0);
     
116. 
     
117. % *%TODO* :
     
118. % Output of FFT is double sided signal, but we are interested in only one side of the spectrum.
     
119. % Hence we throw out half of the samples.
     
120. sig_fft1=sig_fft1_1(1:Nr/2, 1);   % 这个切片从矩阵(1024*128)变向量(1*512)了  只取一列？
     
121. 
     
122. 
     
123. %plotting the range
     
124. figure ('Name','Range from First FFT')
     
125. subplot(1,1,1)
     
126. 
     
127. % *%TODO* :
     
128. % plot FFT output
     
129. 
     
130. plot(sig_fft1);
     
131. %plot(range_resolution*Nr/2, sig_fft1);
     
132. axis ([0 200 0 1]);
     
133. 
     
134. 
     
135. 
     
136. %% RANGE DOPPLER RESPONSE
     
137. % The 2D FFT implementation is already provided here. This will run a 2DFFT
     
138. % on the mixed signal (beat signal) output and generate a range doppler
     
139. % map.You will implement CFAR on the generated RDM
     
140. 
     
141. 
     
142. % Range Doppler Map Generation.
     
143. 
     
144. % The output of the 2D FFT is an image that has reponse in the range and
     
145. % doppler FFT bins. So, it is important to convert the axis from bin sizes
     
146. % to range and doppler based on their Max values.
     
147. 
     
148. Mix=reshape(Mix,[Nr,Nd]);
     
149. 
     
150. % 2D FFT using the FFT size for both dimensions.
     
151. sig_fft2 = fft2(Mix,Nr,Nd);
     
152. 
     
153. % Taking just one side of signal from Range dimension.
     
154. sig_fft2 = sig_fft2(1:Nr/2,1:Nd);
     
155. sig_fft2 = fftshift (sig_fft2);
     
156. RDM = abs(sig_fft2);
     
157. RDM = 10*log10(RDM) ;
     
158. 
     
159. %use the surf function to plot the output of 2DFFT and to show axis in both
     
160. %dimensions
     
161. doppler_axis = linspace(-100,100,Nd);
     
162. range_axis = linspace(-200,200,Nr/2)*((Nr/2)/400);
     
163. figure ('Name','2D FFT'),surf(doppler_axis,range_axis,RDM);
     
164. xlabel('Speed');
     
165. ylabel('Range');
     
166. zlabel('db of rdm');
     
167. 
     
168. %% CFAR implementation
     
169. 
     
170. %Slide Window through the complete Range Doppler Map
     
171. 
     
172. % *%TODO* :
     
173. %Select the number of Training Cells in both the dimensions.
     
174. 
     
175. Tr=10;
     
176. Td=8;
     
177. 
     
178. % *%TODO* :
     
179. %Select the number of Guard Cells in both dimensions around the Cell under
     
180. %test (CUT) for accurate estimation
     
181. 
     
182. Gr=4;
     
183. Gd=4;
     
184. 
     
185. % *%TODO* :
     
186. % offset the threshold by SNR value in dB
     
187. 
     
188. off_set=1.4;
     
189. 
     
190. 
     
191. % *%TODO* :
     
192. %design a loop such that it slides the CUT across range doppler map by
     
193. %giving margins at the edges for Training and Guard Cells.
     
194. %For every iteration sum the signal level within all the training
     
195. %cells. To sum convert the value from logarithmic to linear using db2pow
     
196. %function. Average the summed values for all of the training
     
197. %cells used. After averaging convert it back to logarithimic using pow2db.
     
198. %Further add the offset to it to determine the threshold. Next, compare the
     
199. %signal under CUT with this threshold. If the CUT level > threshold assign
     
200. %it a value of 1, else equate it to 0.
     
201. 
     
202. 
     
203. % Use RDM[x,y] as the matrix from the output of 2D FFT for implementing
     
204. % CFAR
     
205. 
     
206. RDM = RDM/max(max(RDM)); % Normalizing
     
207. 
     
208. % *%TODO* :
     
209. % The process above will generate a thresholded block, which is smaller
     
210. %than the Range Doppler Map as the CUT cannot be located at the edges of
     
211. %matrix. Hence,few cells will not be thresholded. To keep the map size same
     
212. % set those values to 0.
     
213. 
     
214. %Slide the cell under test across the complete martix,to note: start point
     
215. %Tr+Td+1 and Td+Gd+1
     
216. for i = Tr+Gr+1:(Nr/2)-(Tr+Gr)   % 这里只需要一半
     
217.     for j = Td+Gd+1:(Nd)-(Td+Gd)   % trip个数 慢时间维 多普勒维
     
218.         %Create a vector to store noise_level for each iteration on training cells
     
219.         noise_level = zeros(1,1);
     
220.         %Step through each of bins and the surroundings of the CUT
     
221.         for p = i-(Tr+Gr) : i+(Tr+Gr)
     
222.             for q = j-(Td+Gd) : j+(Td+Gd)
     
223.                 %Exclude the Guard cells and CUT cells
     
224.                 if (abs(i-p) > Gr || abs(j-q) > Gd)
     
225.                     %Convert db to power   噪声水平求和
     
226.                     noise_level = noise_level + db2pow(RDM(p,q));
     
227.                 end
     
228.             end
     
229.         end
     
230.         
     
231.         %Calculate threshould from noise average then add the offset
     
232.         threshold = pow2db(noise_level/(2*(Td+Gd+1)*2*(Tr+Gr+1)-(Gr*Gd)-1));
     
233.         %Add the SNR to the threshold   db值就是加offset
     
234.         threshold = threshold + off_set;
     
235.         %Measure the signal in Cell Under Test(CUT) and compare against
     
236.         CUT = RDM(i,j);
     
237.         
     
238.         if (CUT < threshold)
     
239.             RDM(i,j) = 0;
     
240.         else
     
241.             RDM(i,j) = 1;
     
242.         end
     
243.         
     
244.     end
     
245. end
     
246. 
     
247. RDM(RDM~=0 & RDM~=1) = 0;
     
248. 
     
249. 
     
250. % *%TODO* :
     
251. %display the CFAR output using the Surf function like we did for Range
     
252. %Doppler Response output.
     
253. % figure,surf(doppler_axis,range_axis,'replace this with output');
     
254. % colorbar;
     
255. figure('Name', 'CA-CFAR Filtered RDM')
     
256. surf(doppler_axis,range_axis,RDM);
     
257. colorbar;
     
258. title( 'CA-CFAR Filtered RDM surface plot');
     
259. xlabel('Speed');
     
260. ylabel('Range');
     
261. zlabel('Normalized Amplitude');
     
262. 
     
263. 
     
264. 
     
265. 
     
266. 
     

复制代码

zouqiqi

让天下人人学会人工智能！人工智能的前景一片大好！