jens/GnuRadio
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

- improved peak detector

Browse Source 
git-svn-id: http://moon:8086/svn/software/trunk/projects/GnuRadio@434 b431acfa-c32f-4a4a-93f1-934dc6c82436
This commit is contained in:
jens
2019-05-21 18:28:54 +00:00
parent 7cfba591a9
commit e523b779d3
1 changed files with 76 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files
matlab/peak_detect.m
+76 -28
View File
@@ -25,61 +25,109 @@
function retval = peak_detect ()
N = 1024
frames = load ('-ascii', '/home/jens/frames.dat');
[numFrames, N] = size(frames)
Nh = N/2;
kn = 0.1
ks = 0.1;
ks = 0.5;
pa = 0.3;
d_max = 30;
am = 0.1;
as = 0.9;
Ad = 10;
Cm = 1;
Ad = 20;
Cm = 0;
s = zeros(1, Nh);
a = zeros(1, Nh);
Xm = zeros(1, Nh);
Xm = -40*ones(1, Nh);
Xs = zeros(1, Nh);
peak_count = zeros(1, Nh);
close all;
figure
for k=1:128
for k=1:numFrames
% Peak generation
k = fix(Nh*rand() + 1);
if (rand() < pa)
if (s(k) == 0)
a(k) = ks;
s(k) = fix(d_max*rand());
if 0
z = fix(Nh*rand() + 1);
if (rand() < pa)
if (s(z) == 0)
a(z) = ks;
s(z) = fix(d_max*rand());
end
end
end
x = kn*randn(1, N);
for n = 1:Nh
if s(n) > 0
s(n) = s(n) - 1;
f = 0.5-rand() + n;
x = x + a(n)*cos(2*pi*f*(0:N-1)/N);
x = kn*randn(1, N);
for n = 1:Nh
if s(n) > 0
s(n) = s(n) - 1;
f = 0.5-rand() + n;
x = x + a(n)*cos(2*pi*f*(0:N-1)/N);
end
end
X = abs(fft(x))/sqrt(N);
X = X(1:Nh);
X = 10*log10(X.*X);
else
X = frames(k,:);
X = X(1:Nh);
end
X = abs(fft(x))/sqrt(N);
X = X(1:Nh);
X = 10*log10(X.*X);
% Peak detection
Xs = (1-as)*Xs + as*X;
Xm_gate = Xs > (Ad+Xm);
peak_count = Xm_gate.*peak_count + Xm_gate;
peak = peak_count > Cm;
gate_peak = Xs > (Ad+Xm);
gate_upd = Xs > (10+Xm);
peak_count = gate_peak.*peak_count + gate_peak;
peak_cand = peak_count > 0;
% Update Xm at no-peak bins
Xm = (1-peak) .* ((1-am)*Xm + am*X) + peak.*Xm;
Xm = (1-gate_upd) .* ((1-am)*Xm + am*X) + gate_upd.*Xm;
% Filtering
if 0
Nhf = 21;
Ndf = (Nhf-1)/2;
hf = [ones(1, Ndf), 0, -ones(1, Ndf)];
Xf = abs(filter(hf, 1, X)/Nhf);
end
% Output
marker_X = find(peak ~= 0);
numPeaksLast = 0;
while(1)
numPeaks = 0;
for x = find(peak_cand > 0),
_x = x;
while(1)
y0 = X(_x);
xp = min(Nh, _x+1);
xn = max(1, _x-1);
if X(xp) > y0
peak_count(_x) = 0;
peak_count(xn) = 0;
_x = xp;
elseif X(xn) > y0
peak_count(_x) = 0;
peak_count(xp) = 0;
_x = xn;
else
peak_count(xn) = 0;
peak_count(xp) = 0;
numPeaks = numPeaks + 1;
break;
end
end
endfor
if numPeaks == numPeaksLast
break;
end
numPeaksLast = numPeaks;
end
peak = peak_count > Cm;
marker_Y = X.*peak;
marker_X = find(peak > 0);
subplot(2, 1, 1)
plot(1:Nh, X, marker_X, marker_Y(marker_X), '+'); grid; axis([0, Nh, -60, 40]);
subplot(2, 1, 2)