How to obtain each each maximum value from each signal set?

Hi all

Because I am doing a special flow cytometry and I got 8 sets of wave, which means 1 set of wave is 1 cell pass through my sensor.

Now, I want to know how to pick up the maximum amplitude from each sets of wave?

Plus, are there any method to prevent double pick if one set of wave have same max value?

Thanks!

[Mehmed Saad]

what do you mean by cell? cell array?

can you attach that?

[Rik]

This becomes almost trivial once you have separated the cells. This should not be very difficult, since there are large portions with an (almost?) 0 voltage. What have you tried so far to separate them?

And what do you mean with your second question? Doesn't max(abs(wave_data)) do what you want?

[Peng Li]

Need some sort of segmentation first, and within each segment you could use max function to find the maximum.

are signals between each set of waves exactly zero or close to zero?

It is up to you which to keep if you do have multiple numbers that are the same. max function gives you the first one.

[Chen Kevin]

The cell means a typical biological cell.

The principle of detecting cell has been shown with my attachment, what I want to do is to pick up max value in each wave sets and count the cells

You can say that signals between each wave set is close to zero or noise.

I did consider about segmentation, but not sure how to do it and afraid the wrong segmentation will mistake signal.

[Peng Li]

To segment the signal, you can either do a thresholding process, or some analytical approach such as hilbert transform to get the signal envelope i think.

For the first one, you could try a small threshold close to zero, and try to detect a segment with values always below the threshold concecutively for how many point. This would be a candidate of a seperator between two waves.

For the second one, try this https://www.mathworks.com/help/signal/ref/envelope.html to get the envelope. After that, it might be easier to separate the waves i suppose.

[Chen Kevin]

I finally use Hilbert function to make these 8 sets become eight the envelopes.

The next step is how could I get the maximum amplitude from each sets?

Here is the code I am using:

%Generate signal
fs = 1e9;
tc = gauspuls('cutoff',10e3,0.5,[],-40); 
t = -tc:1/fs:tc; 
x = gauspuls(t,10e3,0.5); 
ts = 0:1/50e3:0.025;
d = [0:1/1e3:0.025;sin(2*pi*0.1*(0:25))]';
y = pulstran(ts,d,x,fs);
plot(ts,y)
xlim([0 0.01])
hold on
%Use Hilbert function to create envelope
hs= hilbert(y)
plot(ts,-abs(hs),'r-', ts,abs(hs),'-r', 'LineWidth',1.5)
xlabel('Time (s)')
ylabel('Volatge (a.u.)')
hold off

[Peng Li]

you can apply a threshold on the envelope, upper or lower, so that you can segment the waves. Within each segment, you could apply max function to the original signal to get the maximum. Does this sound clear to you?

[Peng Li]

if not clear simply upload an examplary file here, which would be easier for me to test.

[Chen Kevin]

Hi Peng Li,

May I hear about how do you use "threshold" way to find the max value from each wave sets? Maybe using "if" function?

Here is the signal code:

%Generate signal
fs = 1e9;
tc = gauspuls('cutoff',10e3,0.5,[],-40); 
t = -tc:1/fs:tc; 
x = gauspuls(t,10e3,0.5); 
ts = 0:1/50e3:0.025;
d = [0:1/1e3:0.025;sin(2*pi*0.1*(0:25))]';
y = pulstran(ts,d,x,fs);
plot(ts,y)
xlim([0 0.01])

[Peng Li]

Sorry didn't log on these days. quite stressful working from home...

see below an example using your synthetic data

%% Your Generate signal Copied here
fs = 1e9;
tc = gauspuls('cutoff',10e3,0.5,[],-40); 
t = -tc:1/fs:tc; 
x = gauspuls(t,10e3,0.5); 
ts = 0:1/50e3:0.025;
d = [0:1/1e3:0.025;sin(2*pi*0.1*(0:25))]';
y = pulstran(ts,d,x,fs);
plot(ts,y)
xlim([0 0.01])

%% segmentation
baseline = ismembertol(y, 0, 1e-3);

% convert 0-1 series baseline to segments
incl = baseline ~= 1;
incT = diff([0; incl(:); 0]);
stId = find(incT == 1);
edId = find(incT == -1) - 1;
seg  = [stId(:) edId(:)];

% you could possibly add some condition here to filter out segments that is
% too short
% alternatively (better), you find the between-waves segments, and add
% conditions to check the length of these segments. If too short, this
% might be a wrong detection and you delete it to merge the two waves
% let me know if this is clear for you

%% with each segment, detect maximum
% can use for, or arrayfun for better readability
pind   = nan(size(seg, 1), 1);
for iS = 1:size(seg,1)
    [~, ind] = max(y(seg(iS, 1):seg(iS, 2)));
    pind(iS) = ind + seg(iS, 1) - 1;
end

hold on;
plot(ts(pind), y(pind), 'ro');

Answer 1

[Star Strider]

You have the Signal Processing Toolbox, so use the findpeaks function to get the peak amplitudes. It has a number of options so you can be certain to get only the values you want. If you want to get the amplitudes of the ‘valleys’ use findpeaks on the negative value of the signal vector.

[Chen Kevin]

Hi Star Strider

I tried two ways to put threshold to seperate each wave sets, however, still don't know how to catch the only one maximum value in each wave sets.

The only way I can successfully get the max peaks from each wave sets is using findpeaks(y,fs,'MinPeakDistance', 2.5e-8)

but, I am not a fan for using this function because I can't garuntee that the flow speed of cells will be stable during the real experiment. Are there any alternative method to obtain the max amplitude from each wave sets?

Here are the two ways of code that I put threshold:

%Used findpeaks function

%Generate signal
clc,clear
close
fs = 1e9;
tc = gauspuls('cutoff',10e3,0.5,[],-40); 
t = -tc:1/fs:tc; 
x = gauspuls(t,10e3,0.5); 
ts = 0:1/50e3:0.025;
d = [0:1/1e3:0.025;sin(2*pi*0.1*(0:25))]';
y = pulstran(ts,d,x,fs);
plot(ts,y)
xlim([0 0.01])
%Set the threshold > 0.1
%minValue= min(y);
%y(y<0.1)=minValue;
findpeaks(y,fs,'Threshold',0.1)

%Generate signal
clc,clear
close
fs = 1e9;
tc = gauspuls('cutoff',10e3,0.5,[],-40); 
t = -tc:1/fs:tc; 
x = gauspuls(t,10e3,0.5); 
ts = 0:1/50e3:0.025;
d = [0:1/1e3:0.025;sin(2*pi*0.1*(0:25))]';
y = pulstran(ts,d,x,fs);
plot(ts,y)
xlim([0 0.01])
%Set the threshold > 0.1
minValue= min(y);
y(y<0.1)=minValue;
findpeaks(y)

[Star Strider]

You are not using findpeaks to its fuull advantage.

Try this:

fs = 1e9;
tc = gauspuls('cutoff',10e3,0.5,[],-40); 
t = -tc:1/fs:tc; 
x = gauspuls(t,10e3,0.5); 
ts = 0:1/50e3:0.025;
d = [0:1/1e3:0.025;sin(2*pi*0.1*(0:25))]';
y = pulstran(ts,d,x,fs);

[pks,locs] = findpeaks(y,'MinPeakProminence',0.5, 'MinPeakDistance',10);

figure
plot(ts,y)
hold on
plot(ts(locs), pks, '^r')
hold off
xlim([0 0.01])

producing:

That looks as though it does what you requested. You may need to experiment with the name-value pairs with your actual signal, however it should be possible to use findpeaks to get the cell pulse amplitudes.

If the minimum values (‘valleys’) are important, negate your signal to find them:

[vlys,vlocs] = findpeaks(-y,'MinPeakProminence',0.5, 'MinPeakDistance',10);

The actual values for the ‘valleys’ will be ‘-vlys’. The ‘vlocs’ value does not have to be adjusted at all. The same name-value pair arguments should work for both.

.

[Chen Kevin]

Thanks Star,

May I know the reason why you set 'MinPeakDistance' to 10?

[Star Strider]

My pleasure.

I experimented until I found a value that worked. It only required three iterations to get both name-value pairs reasonably optimised.

Otherwise, I would have requested all 4 outputs from findpeaks to search for the associated values for the peaks and their locations that were the ones I wanted, then used approximations of those for the values. I do that with more complicated problems. This one was relatively straightforward.

[Chen Kevin]

Hi Star,

I think I have the answer now, they were divided by the pulses(waves) in one set.

That's good, will put this function in account for future use.

One more extra question is, how could I set a code to prevent double count once there are more than 1 same max value in 1 set of wave?

[Star Strider]

‘One more extra question is, how could I set a code to prevent double count once there are more than 1 same max value in 1 set of wave?’

I do not know what this refers to, since ther is no illustration of it. I would assume the 'MinPeakDistance' value would need to be adjusted to account for that. (I am familiar with flow cytometry, however I have no direct experience with it or processing signals from it.)

[Chen Kevin]

Thanks, I think I will figure out that.

Sorry for annoying, but here is one extra question: How to measure width of each wave sets?

[Star Strider]

As always, my pleasure!

That is likely best done after using the envelope function. If you use findpeaks on the signal envelopes, it will output the full-width-half-maximum values for the peaks it detects, if you request that output. If you want the widths of the peaks at the some other value, that is more complicated and depends on the characteristics of the signals themsleves. There is no general way to do that.

.

[Chen Kevin]

What I have considered method by used for and if function is :

1. Set a threshold for switch on and switch off record the data
2. Once the amplitude (y>threshold) which means the cell just enter the sensing zone, start the record
3. When the amplitude (y<threshold) which means the cell almost leave the zone, stop recording
4. Find the max value from this recording time interval
5. The width can be get when the time interval calculated

However, I am quite fresh in Matlab so don't know how to wirte down them...

[Star Strider]

I do not understand. If you want to do this in real time, I doubt MATLAB is capable of tha. That might require additional hardware. If so, I will not be able to help you with it.

If you are only interested in recorded signals, there are several functions that might be able to help you with that, including findchangepts, and ischange.