Finding a maximum signal variation in MATLAB

I have a signal which consist of one 1000000x1 array (data) and one 1000000x1 array (time). I'm having problem to detect the location where the signal have the largest variation with respect to the value of the array.

For example in this signal, I want to locate the variation of the signal from time 4 until 11. I found out that findpeaks command will detect all peaks but I want to detect only the highest peak(highest value from the array) and the data (eg: 5000 sampling) to the right and left from the highest peak.

Thank you in advance.

Answer 1

[Star Strider]

If you are only finding one peak each time, I would use the max and min functions, each with two outputs, to get the absolute maximum and minimum. The second output is the index location of the maximum or minimum.

You can also use the Signal Processing Toolbox findpeaks function on the original signal to get the value and location of the highest peak, the use findpeaks again on the negative of the original signal to get that highest value and location.

[Jimmy]

Thank you for the reply. Let say I'm using the max and min functions approach, how do I find the time of that particular max value

example:

data = [3 5 1 7 6 2 1]
time = [0 1 2 3 4 5 6]

since the max value of data is 7 which is at time = 3, how do I get the max value for that time?

[Star Strider]

My pleasure.

The code:

[data_max, max_idx] = max(data);
time_at_max = time(max_idx);

data_max =

       7

max_idx =

       4

time_at_max =

       3

The second output of max gives you the index of the maximum, so that will be the same index for the corresponding time vector as well. Here, ‘data_max’ is the maximum of ‘data’ and ‘time_at_max’ is the time the maximum occurs.

Use essentially the same syntax for the minimum, using the min function.

[Jimmy]

It worked. Thank you for your guidance.

One more inquiry if I may ask. Since now I found the max value of the signal, how to 'cut' the signal between those max value.

From the previous example, if the max_idx is 4, is the any method on how to take the two value from the left and right of the max value? which mean the left value will be 5 & 1 and the right value will be 6 & 2.

Thank you

[Star Strider]

My pleasure.

The two ‘left’ and ‘right’ values would be:

left2  = data(max_idx-[2:-1:1]);
right2 = data(max_idx+[1:2]);

Is that what you want?

[Jimmy]

I'm sorry but the output is not exactly what I want because the 'left'and 'right' values must also with respect to time in which 5 & 1 for time 1 & 2 and 6 & 2 for time 4 & 5.

Thank you

[Star Strider]

To get the times, do the same as for the maximum and minimum times.

left2_time = time(max_idx-[2:-1:1])
right2_time = time(max_idx+[1:2])

The same index references apply.

[Image Analyst]

But doesn't this get windows around each narrow peak in the burst, whereas he wants one single window around the whole burst?

[Jimmy]

Thank you. It worked. If I may ask how do I combine the left value and the right value starting at the max value because now I'm using the following coding:

plot(left2_time,left2);  %left value
plot(right2_time,right2); %right value

I want to construct those values in one graph. Am I using the wrong coding?

Thank you

[Star Strider]

Actually what you’re currently doing works, providing you use the hold function. I’m not sure what final result you want, because you haven’t described it. We’re doing this a bit at a time.

To combine all of them, do a simple horizontal concatenation:

time_vct = [left2_time  NaN  right2_time];
data_vct = [left2  NaN  right2];

The NaN is there to provide discontinuity, since I get the impression that you want to plot these separately.

If you want the entire sequence displayed at once, this works:

win_idx = (max_idx-2):(max_idx+2);
plot(time(win_idx), data(win_idx))

The entire code (in summary for continuity) now becomes:

data = [3 5 1 7 6 2 1];
time = [0 1 2 3 4 5 6];

[data_max, max_idx] = max(data);
time_at_max = time(max_idx);

left2  = data(max_idx-[2:-1:1]);
right2 = data(max_idx+[1:2]);

left2_time = time(max_idx-[2:-1:1]);
right2_time = time(max_idx+[1:2]);

time_vct = [left2_time  NaN  right2_time];
data_vct = [left2  NaN  right2];

figure(1)
plot(time_vct, data_vct, 'LineWidth',2)
grid

win_idx = (max_idx-2):(max_idx+2);

figure(2)
plot(time(win_idx), data(win_idx))
grid

[Jimmy]

I finally got the final result that I want. Thank you so much for your help.

The final coding are as follows:

load('SigData.mat');
[data_max,max_idx] = max(data);
time_at_max = time(max_idx);
left2  = data(max_idx-[200000:-1:1]);
right2 = data(max_idx+[1:200000]);
left2_time = time(max_idx-[200000:-1:1])
right2_time = time(max_idx+[1:200000])
win_idx = (max_idx-200000):(max_idx+200000);
plot(time(win_idx), data(win_idx))

The element data = [3 5 1 7 6 2 1] and time = [0 1 2 3 4 5 6] is just an example because my real data (SigData.mat) consist of 1000000x1 array. So with your help and guidance, I managed to get the exact result as I wanted by finding the max value and the 'left' and 'right' value. I really appreciated the additional tips and hints.

Again, thank you so much. Have a nice day

[Star Strider]

As always, my pleasure!

[Jimmy]

Greetings, Regarding my previous example, is it possible to perform Fast Fourier Transform on that signal? Let's assume the sampling frequency is 1000Hz. This is the coding that you have shown earlier to cut the signal. The attached file contains the original data. Thank you in advance

 load('TestSig.mat');
[data_max,max_idx] = max(data);     
% time_at_max = time(max_idx);      
% left2  = data(max_idx-[120000:-1:1]);
% right2 = data(max_idx+[1:200000]);   
% left2_time = time(max_idx-[120000:-1:1])
% right2_time = time(max_idx+[1:200000])  
win_idx = (max_idx-120000):(max_idx+200000);
plot(time(win_idx), data(win_idx))
xlabel('Time(s)');
ylabel('Data')

[Star Strider]

Sure. Just do the fft on the entire signal. Doing it on only part of the signal will introduce undesirable discontinuities that would likely show up as ‘ringing’ in the fft.

If you want to filter out noise, I would use a bandpass filter. You can get the approximate bandpass frequencies from the fft. (The code between the first two images in the fft documentation will provide you everything you need to get the frequency scaling correct.) Then choose the stopband frequencies heuristically to get the result you want. Begin by defining them as [0.5 1/0.5].*Wp where ‘Wp’ are the passband frequencies, and narrow them to [0.8 1/0.8].*Wp as necessary to get the behaviour you want. My filter design procedure is here.

[Jimmy]

Thank you for the reply. So I have to design a bandpass filter in order for me to design the fft is it? If so then how am I to do the fft on my data array and time array?

[Star Strider]

My pleasure.

Well, no. That’s backwards. Do the fft first, then use that information to design your bandpass filter. You do the fft on your data only. Your time vector is used to define the frequency vector for your fft. It otherwise doesn’t enter into the fft procedure.

[Jimmy]

I see. I've tried several attempts to perform the fft. My coding are as follows:

 load('TestSig.mat');
y = data;
Fs = 1000;              % Sampling frequency
L = length(y);
NFFT = 2^nextpow2(L);
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(Y(1:NFFT/2+1))) 
xlabel('Frequency (Hz)')
ylabel('data')

And the figure are as follows:

I'm not sure whether my coding is correct or not. Is my sampling frequency (Fs) too small in order to plot this fft? Thank you

[Star Strider]

Your code looks correct. You need to play with the axis limits of the figure so you can see the interesting parts of the spectrum. Add this line after the plot:

axis([0  10    ylim])

Experiment with the x-axis upper limit (here 10) to get the resolution you want.

[Jimmy]

Great. My figure looks much better after implementing the axis limit. Thank you so much for your help.

[Star Strider]

My pleasure.

If you have any problems with the filter design, I’ll do what I can to help.

[Jimmy]

Greetings, I would like to ask you one question.

I have one set of element such that:

    n = [9 7 4 8 13 2 6];

How do I find the first minimum value of this element? By right 2 is the minimum value but I want to get the value 4 for my first minimum value.

Thank you in advance.

[Image Analyst]

Here's something that works:

n = [9 7 4 8 13 2 6];
firstMinIndex = find(diff(n) > 0, 1, 'first')
firstMinValue = n(firstMinIndex)

If you have the Image Processing Toolbox, you can also do this:

allMinIndexes = find(imregionalmin(n))
firstMinValue = n(allMinIndexes(1))

[Star Strider]

To use findpeaks, negate your signal, then choose the first peak:

n = [9 7 4 8 13 2 6];
[pks,lcs] = findpeaks(-n);
FirstMin = n(lcs(1))

FirstMin =
     4

[Jimmy]

It worked. Thank you so much for your help.

[Star Strider]

My pleasure!

[Jimmy]

Greetings, I have a question that I would like to ask.

I have a set of element such that :

n = [-2.34 -3.87 -40.02 -13.91 -30.85 -19.49 -27.36 -17.60 -28.74 -17.94]
end

I want to find the min position between range 6 - 10 in which the value of -19.49 -27.36 -17.60 -25.74 -17.94 . By right, the min value should be -28.74 which is at position 9. My coding are as follows:

    [MinValue,MinPosition] = min(n(6:10));

However, when I check the MinPosition, it is at position 4 because I started at 6. I wish to find the min position of the element starting from position 1 so that the actual position will be at position 4. Please help me.

Thank you in advanced

[Star Strider]

Your Comment is a bit confusing. I believe you want the ‘MinPosition’ value to be 9, not 4. Since your selected vector subset begins at 6, add 5 to ‘MinPosition’ to get the correct position in the vector.

So add 1 less than the starting position of the subset vector index range to get the absolute index value.

[Jimmy]

Yes my mistake. It was supposed to be 9.

One more question, I would like to do high pass filter on a signal. Sampling frequency of 50000, cut-off frequency at 0.4. My coding are as follows:

n = signal;
Fc = 0.4;
Fs = 50000;
order = 3;
Wn = 2*Fc/Fs;
[c3,c2] = butter(order,Wn,'high');
NewNormalised = filter(c3,c2,n);

I would to confirm whether my coding for highpass filter is correct or not. Thank you in advanced

[Star Strider]

My pleasure.

That gives a passband frequency of 10 kHz. Be aware that high pass filters will include all the noise up to 25 kHz (the Nyquist frequency). I would use a bandpass filter instead to pass the frequencies of interest and limit noise as much as possible. Use the fft of your signal to decide on the passband.

My filter design procedure is here.

[Jimmy]

Figure below shows my signal and I want to filter out signal within the red area at the cut-off frequency approximately 0.1-0.5 Hz. The sampling frequency is 50000Hz. So do my highpass filter coding is correct or not? I would really appreciate your review and recommendation. Thank you in advanced.

[Star Strider]

You probably will not be able to get a stable filter with a Butterworth design, 0.1-0.5 Hz cutoff frequency and order = 3. You probably need to use a Chebyshev design. With a Nyquist frequency of 25 kHz, the normalised cutoff frequency Fc for the more realistic 0.5 Hz is 0.5/25000 = 2E-5, not 0.4.

Be certain to use the second-order section representation for your filter, and use freqz to be certain that it’s stable and gives you the characteristics you want.

Please see my filter design link that I posted earlier. That should tell you everything in detail.

[Jimmy]

Thank you for the tutorial. I managed to do my filter design. Basically, my signal processing methods for this signal consists of high pass filter, FFT, and power spectrum density. My next procedure is to perform the principal component analysis. I found this matlab function: [pc,score,latent,tsquare] = princomp(X);

What is the meaning of pc,score,latent and tsquare? Which output should I use to get the result?

Thank you in advance.

[Star Strider]

My pleasure.

With respect to the principal component analysis, the documentation is clear on what the outputs are for princomp. The interpretation in the context of your problem entirely depends on what ‘X’ is. I have used PCA in the past with the fft of my signals to determine the frequencies that best separated them so I could classify them more efficiently. (I was using them on the fft of electroencephalography data to determine — with statistically significant accuracy and repeatability — the task a person was performing. We published those results in the EEG journal in 1995.)

This is actually a new question. If you post it as such, describe in detail what ‘X’ is in your application, what you are studying, and the hypothesis you are testing. The more information you provide, the more information you will receive in return. If you include some or all of your data in your Question, attach it to your original post, rather than a subsequent comment.

[Jimmy]

I have 30 sets of data containing human motion signals (walking, running, etc). Each signal consists of time-domain signal of one 1000000x1 array (data) and one 1000000x1 array (time). I am required to extract information and determine whether the signal is either running or walking. From my signal processing procedures, I've done high pass filtering, FFT and power spectrum density(PSD) in which my current final result is a matrix of 1136364x1 (let's consider it as X). Now I am required to perform PCA to clustered the motion signals. All 30 sets of data will undergo the same procedures.

I'm trying to use the function princomp to get the PCA. So my coding to execute the PCA from my final PSD result (X) would be:

[pc,score,latent,tsquare] = princomp(X);

How do I use the [pc,score,latent and tsquare] to get the result? Thank you.

[Star Strider]

I would have done bandpass filtering to eliminate the d-c offset and any low-frequency baseline variations, as well as high-frequency noise. A Butterworth design would likely work for your data.

As I mentioned, I don’t have much recent experience with PCA, and without your data it would be difficult to give you a specific answer as to how PCA works with it.

You mentioned that your PSD records are column vectors (together comprising ‘X’), with your original data being time-domain data, and the data you’re presenting to the PCA being frequency-domain data, again as a column vector for each study subject. As I read the princomp documentation, you might have to transpose ‘X’ to use princomp correctly. The documentation states that the rows are observations (here, subjects either running or walking), and columns are variables (here, frequencies).

As I understand your study, your goal is to classify the FFT (or PSD) of each subject as either running or walking. The PCA would be appropriate for this, but you have to be certain you present your data correctly to the PCA, or your results will not be optimal.

It may be best for you to submit this as a new Question, and include the relevant parts of your code and at least a representative sample — if not all — of your data.

Save your data to a .mat file so we can simply load it and don’t have to write code to import and then parse a text or other file.

[Jimmy]

Thank you so much for the tutorial on performing fft. It worked out fine. I have a question regarding my previous signal. Is it possible to perform haar wavelet transform to the signal instead of using fft? I'm sorry as I dont have much knowledge on wavelet transform.

Is there any restriction in terms of the size of matrix array when performing haar wavelet because my signal consist of 1000000X1 array.

[Star Strider]

I have some experience with wavelets (and the Wavelet Toolbox), but I’ve not used them (or it) in a bit.

I finally had to go online to find Wavelet Families (it used to be easy to find in the Wavelet Toolbox Documentation, and it is less descriptive than previous versions). I would use a Daubechies wavelet, since they seem to more appropriately approximate the shape of your signal.

Note — In my absolutely not humble opinion, the documentation appears to be getting more difficult to follow rather than easier in the most recent releases. Look up the Wavelet Toolbox documentation for R2015a. It might be easier to follow. Then check with the current documentation for any functions you want to use, since the function behaviour may have changed.

Most signal processing functions don’t care how long your signal vectors are. Your computer hardware does.

[Jimmy]

Thank you for the reply. Regarding my signal processing technique, I am required to perform Haar wavelet transform. I'm having difficulties to use the scaling function and keep getting errors as I don't know what is the best scaling function to use for 1000000 data. Attached below is the coding for haar wavelet which I found from several references.

%assuming InpSig = data (which consists of 1000000x1 array)
Fs = 50000;
N = length(InpSig);                         %number of data points
freqRes = Fs/N;                             %frequency resolution
freqAxisLab = freqRes*[0:N-1];              %calculate the labels for the frequency axis
fc = centfrq('haar',1);                     %scaling function for haar is calculated based on central frequency and level of iteration
freqrange = [1 freqAxisLab];
scalerange = fc./(freqrange*(1/Fs)); 
scales = scalerange(end):0.05:scalerange(1);%scalerange is the initial value:steps of the scale:end value of the scale
inpSigWAV = cwt(InpSig,scales,'haar','plot');      %CWT calculation
cwt(InpSig, scales,'haar','plot');
figure
plot(freqAxisLab, inpSigWAV,'b')

I was wondering whether my Haar wavelet coding is correct or not? Thank you in advance.

[Star Strider]

I do not have enough recent experience with wavelets to provide the guidance you need.

It would be best if you post this as a new Question.

[Jimmy]

I see. Ok noted. Thank you for the tutorial and guidance that you have been given since I posted this question. It was a great help.

Answer 2

[Image Analyst]

Just threshold and call bwareaopen(). See my answer here: http://www.mathworks.com/matlabcentral/answers/231181#answer_187244

[Jimmy]

Thank you for the reply.

If I have a signal as shown below, how do I extract the signal only from 1 to 2 because this signal is from 1000000x1 array (x-axis) and 1000000x1 array (y-axis)

By the way, can you explain on how to use the bwareaopen() command?

Thank you

[Image Analyst]

You might have to take the mean of the last chunk of your signal and then subtract that from the whole signal so that the baseline is along the y=0, x axis. Then doing abs() flips the negative signal so it's all positive. Then the thresholding gets you all the elements above some level, like 0. But the problem is in that pulse/burst there are still elements that dip down to zero. However they are short, maybe only an element or two long. But if we want a window around the pulse, we need to get rid of those short dips in the main pulse. To do that we run bwareaopen() on the thresholded/binary/logical signal. We tell it to get rid of short segments of a few elements and just leave the longer ones. That way we'll clean up all the short noise and be left with just the overall main window of the whole burst. I would have tried it on your data but you didn't attach it.

[Image Analyst]

Attach your data if you want me to run the code on your data.

[Jimmy]

The attached file contains the 1000000x1 data array and 1000000x1 time array.

Figure below shows the output of the signal from the coding: plot(time,data)

I want to extract the signal only from area 1 to 2

Thank you

[Jimmy]

This is the attached file.

Thank you

Answer 3

[Greg Dionne]

You could also try envelope:

envelope(data-mean(data),20000,'rms')
rmsEnv = envelope(data-mean(data),20000,'rms');
envTime = time(rmsEnv>0.02);
envData = data(rmsEnv>0.02);
plot(envTime, envData)

[Image Analyst]

What toolbox is that in?

[Star Strider]

It’s a FEX contribution. There appear to be several like it.

[Image Analyst]

envelope() is now in the Signal Processing Toolbox. I haven't played around with it yet so I don't have a demo on this function.