- added finding peak width

git-svn-id: http://moon:8086/svn/software/trunk/projects/GnuRadio@445 b431acfa-c32f-4a4a-93f1-934dc6c82436

matlab/peak_detect.m

@@ -38,10 +38,12 @@ pa = 0.3;
 d_max = 30;
 
 am = 0.04;
-as = 0.9;
+as = 0.2;
-Ad = 15;
+Ad = 20;
 Cm = 0;
-pbh = 5;
+
+pb_min = 5
+pb_max = 37
 
 s = zeros(1, Nh);
 a = zeros(1, Nh);
@@ -83,11 +85,12 @@ for k=1:numFrames
 	% Init 
 	if k == 1
 		Xm = Xm + mean(X);
+		Xs = Xs + mean(X);
 	endif
 	
 	% Peak detection
 	Xs = (1-as)*Xs + as*X;
-	gate_peak = Xs > (Ad+Xm);
+	gate_peak = (Xs-Xm) >= Ad;
 	peak_count = gate_peak.*peak_count + gate_peak;
 
 	% Output
@@ -133,11 +136,50 @@ for k=1:numFrames
 	peak_boxes_rel = zeros(1, Nh);
 	peak_boxes_abs = zeros(1, Nh);
 	for x = peak_cand_sorted,
+		
+		% Find peak widths
+		left = x;
+		right = x;
+		nom = X(x) - 10;
+		left_found = 0;
+		right_found = 0;
+		while(~(left_found & right_found))
+			% Determine peak width left from center (x)
+			if X(left) > nom
+				if left > 1
+					left = left - 1;
+				else
+					left_found = 1;
+				end
+			else
+				left_found = 1;
+			endif
+			% Determine peak width right from center (x)
+			if X(right) > nom
+				if right < Nh
+					right = right + 1;
+				else
+					right_found = 1;
+				end
+			else
+				right_found = 1;
+			endif
+			pbh_left = x - left;
+			pbh_right = right - x;
+			
+			% Take larger peak width
+			pbh = max(pbh_left, pbh_right);
+			
+			% Ensure peak box is at least pb_min
+			pbh = max(pbh, pb_min);
+		end
+		
+		% Use peak width for peak box
 		pbox_start = max(1, x-pbh);
 		pbox_end = min(Nh, x+pbh);
 		
 		% Find intersection
-			box_height_rel = Xs(x) - Xm(x);
+		box_height_rel = X(x) - Xm(x);
 		box_height_abs = X(x);
 		does_intersect = 0;
 		for n=pbox_start:pbox_end
@@ -154,14 +196,13 @@ for k=1:numFrames
 	endfor
 
 	% Update Xm at no-peak bins
-	gate_upd = peak_boxes_rel > 0;
 	Xm = (1-am)*Xm + am*X;
 
 	% Fill peak gaps with guessed values
 	last = Xm(1);
 	for n=2:Nh,
-		if gate_upd(n) == 0
+		if (Xm(n) < (last+2)) or (peak_boxes_rel(n) == 0);
-			last = 0.5*last + 0.5*Xm(n);
+			last = (1-am)*last + am*Xm(n);
 		else
 			Xm(n) = last;	
 		end
@@ -175,15 +216,8 @@ for k=1:numFrames
 	subplot(2, 1, 1)
 	plot(1:Nh, X, peak_final, X(peak_final), 'ro', 1:Nh, peak_box_values, 'm-'); grid; axis([0, Nh, -60, 40]);
 	subplot(2, 1, 2)
-	plot(1:Nh, Xs, 1:Nh, Xm); grid; axis([0, Nh, -60, 40]);
+	plot(1:Nh, Xs-Xm, 1:Nh, Xm); grid; axis([0, Nh, -60, 40]);
 	pause(1/30);
 end
 
-function s = matchedFilter(N, B)
-n2 = (N-1)/2
-n = (0:N-1)-n2
-s = sinc(n/B)
-
-endfunction
-	
 endfunction