CODEWITHSUNDEEP
matlabsignal-processing
dpdp
dpdp

Reputation: 322

finding peaks that have some structure in 1-d

I have signal with peaks that has some structure on top of some background, and I'm trying to find a robust way to locate their positions and amplitudes.

For example, assume the peak has this form:

t=linspace(0,10,1e3);
w=0.25;
rw=@(t0) 2/(sqrt(3*w)*pi^0.25)*(1-((t-t0)/w).^2).*exp(-(t-t0).^2/(2*w.^2));

and I have several peaks a bit too close so their structure starts to interfere with the others, and one separated, so you can see it's structure:

pos=[ 2 2.5 3 8]; % positions
total_signal=0.*t;
for i=1:length(pos) 
   total_signal=total_signal+rw(pos(i));
end
plot(t,total_signal);

screenshot

So the goal is to find all peaks found in total_signal and check that their positions agree with the original positions that were used to generate them in pos.

Upvotes: 1

Views: 165

Answers (1)

Cris Luengo
Cris Luengo

Reputation: 60514

Your signal can be thought of as the convolution of a pulse train with the peak shape. Deconvolution can be used to retrieve the pulse train, whose peaks are the locations you are looking for. This assumes that the peak shape is known and constant.

I will start by rewriting your code as follows:

t = linspace(0,10,1e3);
w = 0.25;
rw = @(t) 2/(sqrt(3*w)*pi^0.25)*(1-((t)/w).^2).*exp(-(t).^2/(2*w.^2));

pos = [2, 2.5, 3, 8];
total_signal = zeros(size(t));
for i=1:length(pos) 
   total_signal = total_signal + rw(t-pos(i));
end
total_signal = total_signal + randn(size(t))*1e-2;
plot(t,total_signal);

I've changed rw to take t-t0, rather than just t0. This will allow me later to create a clean signal with just one peak in the middle. I have also added noise to the signal, for a more realistic problem. Without the noise, the problem is a lot easier to solve.

Wiener deconvolution is the simplest approach to solve this problem. In short, we assume that

total_signal = conv(pulse_train, shape)

In the frequency domain, this is written as

G = F .* H

(with .* the element-wise multiplication, using MATLAB syntax here). The Wiener deconvolution is:

F = (conj(H) .* G) ./ (abs(H).^2 + k)

with k some constant that we can tune to regularize the solution.

We implement this as follows:

shape = rw(linspace(-5,5,1e3));
G = fft(total_signal);
H = fft(ifftshift(shape)); % ifftshift moves the origin to sample #0, as expected by FFT.
k = 1;
F = (conj(H) .* G) ./ (abs(H).^2 + k);
pulse_train = ifft(F);

Now, findpeaks (requires the Signal Processing Toolbox) can be used to find the prominent peaks:

findpeaks(pulse_train, 'MinPeakProminence', 0.02)

output, showing 4 detected peaks at close to the same height

Note how the four peaks are approximately the same height. It's not exact, because we're regularizing to deal with the noise. An exact solution is only possible in the noise-free case. Without noise, k=0, and the expression simplifies to F = G ./ H.

Also, the x-axis is off in the plot produced by findpeaks, but this shouldn't affect the results. The locations returned by it are indices into the array, those same indices can be used to index into t and find the actual locations of the peaks.

Upvotes: 1

Related Questions