determine location of certain turning points of curve

I'm working on a set of running gait kinematics data and I'm supposed to find the maximum and minimum points of the curve to get the maximum and minimum degree of motion at each part of the gait cycle.

I've managed to do so with the help of the function peakdex (<http://www.billauer.co.il/peakdet.html>) This is my code:

[maxpts, minpts] = peakdet(hipangle, 0.2);
hold on; plot(minpts(:,1), minpts(:,2), 'g*');
plot(maxpts(:,1), maxpts(:,2), 'r*');
plot(hipangle);

However, now i'm not very sure how should i go about separating the data such that i can separate them into different gait cycles i.e. made up of 2 min points and one peaks (one is larger than the other). Can i get some suggestions on how I can do so?

Currently, what I have in mind is to (1) obtain all the minimum points (those nearer to the baseline) and storing them into an array minpts_low, and (2)use a for loop to obtain the values and positions of maxpts within 2 minpts_low and put them into another matrix gaitcycles. So, gait cycles will be made up of 2 cols (first maxpt and 2nd maxpt) and each row represent one gait cycle.

This is the code i've come up with based on my idea:

minpts_low = [];
  % minpts_low = Inf;
  for i= 1:length(minpts)
      if (minpts(i+1,1) - minpts(i,1)) <= 0.03 
          minpts_low(j) = i;
      end
  end

gaitcycles_hip = [Inf,Inf];
for k = 1:length(maxpts)
    for j = 1:length(minpts_low)
        if (maxpts(k,1) > minpts_low(j,1)) && (maxpts(k,1) < minpts_low(j+1,1))
            gaitcycles_hip = [maxpts]; % i put maxpts only as maxpts is made up of 2 cols
(values position)
        end
    end
end

However, I have some problems translating my ideas into code. I'm not sure how to start with part (2) of my idea. For part (2), what I get from my code is Inf Inf, which means that my for loop didnt work. :/

As for part 1, I get the error "Attempted to access minpts(30,1); index out of bounds because size(minpts)=[29,2]". I cant define the size of minpts beforehand as I'm not sure what the size will be.

really appreciate if anyone can offer some suggestions regarding this problem! :)

[Star Strider]

It would help if you attached the time-hip angle data you plotted. It’s difficult to explore various possibilities without it.

[Mich]

Hi, this is the data that I had used.

[Star Strider]

What is the sampling time or sampling frequency?

[Mich]

100Hz

[Star Strider]

Thank you. That will help if I need to design a filter.

[Mich]

Can I check the rough approach I should take if i were to design a filter ? I've actually considered this option, based on the method that you had helped me with in my previous qn (<http://www.mathworks.com/matlabcentral/answers/159212-finding-value-that-satisfies-2-conditions)>. However, i realised that I will not know which specific threshold or specific quantitative pattern to consider when designing the filter due to the high variability in the data.

Answer 1

[Star Strider]

The filter is an excellent idea — and the approach I have implemented. I kept everything I used in developing my code including the fft call to present my approach to the problem and to illustrate the frequency content of your signal. The filter design is relatively straightforward. (I used freqz for quality assurance — the filter is stable, well-behaved, and does what I designed it to do. I didn’t include that call here.)

The signal processing filters out high-frequency noise and the baseline drift, giving a stable signal without losing any of the signal characteristics you probably want to preserve. I used the findpeaks function to detect the peaks (part of the Signal Processing Toolbox), and negated your signal to flip it upside-down to detect the troughs as peaks to detect what appear to be the onsets of each gait cycle, then flipped it back to plot it. I set an arbitrary threshold and plotted what seem to me to be ‘valid’ strides with blue stars. Note that my ‘gaitlocs’ variable contains the indices of the ‘peaks’ (troughs) that meet the criterion.

ha = dlmread('Mich_hipangle.txt');
Fs = 100;                                           % Sampling Frequency (Hz)
Fn = Fs/2;                                          % Nyquist Frequency
Ts = 1/Fs;                                          % Sampling Interval (s)
Tv = linspace(0, (length(ha)-1)*Ts, length(ha));    % Time Vector

ftha = fft(ha)/length(ha);
Fv = linspace(0, 1, length(ha)/2+1)*Fn;
Iv = 1:length(Fv);

figure(1)
plot(Fv, abs(ftha(Iv)))
grid
axis([0  5    ylim])
xlabel('Frequency (Hz)')
ylabel('Amplitude')

% Design & Implement Bandpass Filter
[n,Wn] = buttord([0.75 3.5]/Fn, [0.5 4.0]/Fn, 1, 10);   
[z,p,k] = butter(n,Wn);
[sos,g] = zp2sos(z,p,k);
haf = filtfilt(sos,g,ha);                           % Filtered Data

[pks,locs] = findpeaks(-haf);                       % Find Minima Of Inverted Data
gaitlocs = locs(pks >= 10);                         % Detect Stride End-Points

figure(2)
plot(Tv, haf)
hold on
plot(Tv(locs), -pks, '+r')
plot(Tv(gaitlocs), haf(gaitlocs), 'pb', 'MarkerFaceColor','b')
hold off
grid
title('Filtered Signal')
xlabel('Time (s)')

[Mich]

hi, thanks! the blue points are what I would like to obtain.sorry for the late reply, i was trying to understand the filters used.

Can I check how did you derive at the stopband and bandpass parameters for the butterworth filter? I'm also not sure what is the purpose of using linspace. As f the filtfilt(zero phase filter), can I say that the filter serves to make the baseline as zero instead of a positive value in the original value?

[Star Strider]

My pleasure!

The indices of the blue stars are given by my ‘gaitlocs’ variable. The third plot call in figure(2) illustrates their use to define the time and the hip angle data of the filtered signal. To get the hip angles of your original data at the positions of the blue stars, use ha(gaitlocs).

Real-world filter design is heuristic. I designed the filters using the results of the fft call, the only reliable way to design filters. One other thing I did in the process was to subtract the mean of the signal before taking the fft in order to discover the frequency of the baseline drift. (I omitted that step as well as the freqz call in the code I posted.) Then I designed the lower stopband to eliminate it and the constant (mean) offset, resulting in a mean of zero. There was nothing above about 4 Hz, so that defined the upper stopband. The other steps in the filter design and implementation are routine.

I use linspace to be sure the length of the resulting vector is what I want it to be. Personal preference.

The filtering step simply makes it easier to detect the positions of the gait limits. It does not alter their positions w.r.t. time, so you can use those indices with your original data.

[Mich]

sorry actually i am still pretty confused over the filters used. actually can i check why is a bandpass filter used instead of a high pass filter only? could i use the high pass filter only to filter out the low freq signal that is causing the drift in the signal? can i also check if you mean this: fft(ha- mean(ha)). Im not sure how I should go about with the subsequent steps of designing a filter after that. could i trouble you to explain more?

[Star Strider]

I use bandpass filters almost exclusively because they filter out noise on the high end and filter out the D-C offset and baseline drift on the low end. Using only a highpass filter would filter out the baseline drift and D-C offset, but would also pass all the noise (up to the Nyquist frequency).

If you want to see the low frequency signals clearly, subtracting the mean from the signal or from the fft of the signal (either will produce the same result after you calculate the fft) will show them. So fft(ha-mean(ha)) would do this.

When I design a filter, I start from the fft of the signal I want to process. That tells me the frequencies that are present in the signal and where to start designing its passbands and stopbands. The choice of the filter type (Butterworth, etc.) depends on what I want the filter to do and the tradeoffs I’m willing to accept. More often than not, I try a few filter designs before I settle on one that gives me the results I want. (The exception is a filter for EKG signals, where the accepted standard passband limits are 0.05 Hz to 100 Hz.) I use the second-order-section (SOS) implementation if possible, and I always check the stability with freqz to be sure the filter has the characteristics I want, before I use it. I use filtfilt for the actual filtering.

The filter I designed is about optimal for your signal. You can probably use it (and the rest of my code) for all your data, providing your sampling frequency doesn’t change.

Filter design and digital signal processing (DSP) in general are more detailed than I can possibly explain here. I encourage you to take a course in DSP when you have the opportunity. MATLAB also has a number of online resources that go into this in detail.

[Mich]

thank you for your explanation! I've read up more about high freq noises which are prevalent in electric devices.

can i check if this is the correct process of deriving the parameters?

qq = fft(ha-mean(ha))/length(ha)
plot(Fv, abs(qq(Iv)) %this shows the range of frequencies in the signal, Fv just changes the scale of the x axis from time points to secs

this gives me the graph which i've attached. i'm still unsure of how this graph provides an indication of the passband and stopband limits that i can use. As for the function freqz, im still unsure of how it can test for stability. i read from the matlab page that it is usually used as freqz(sos,no of evaluation points). So for this case, I would assume that the my input arguments are freqz(sos,1500).

However, this code doesnt work and there is an error that says that the numerator and demominator of the function must be stored in vectors. Im confused regarding this part as based on my understanding, the numerators and denominators are stored in sos.

As for filtfilt, i understand that it is used to correct for phase distortion for the filter. however, i read that one fault of the zero phase filtering is that causality is lost and this has result some small peaks being detected during the calibration period (the first 5 seconds). im still unable to appreciate the use of a zero phase filtering, and tried to plot the graph prior to the application of filtfilt but was unable to do so.

With regards to the usefulness of the filter, i agree that it is very applicable for me to help me process my data, e.g. for my ankle angle data. However, for hip angle in particular, i would like to keep the original baseline (without the drift) because based on physiological data, the baseline should be positive with the majority of the data in the positive range.

[Star Strider]

That looks correct. However your interpretation of my ‘Fv’ is not quite correct.

‘Fv’ is the frequency vector in Hz. (The other vector, ‘Iv’, is an index vector that makes it easier to plot the fft.) The passband cutoffs for my filter are 0.75 Hz and 3.5 Hz (in the filter design they’re then normalised to the Nyquist frequency as the filter design functions require), so you can see from your plot and those frequencies that it would include your frequencies of interest, while eliminating the low-frequency baseline drift signal and the high-frequency noise. It keeps as much of the data as possible that actually characterise your hip angle information, and makes your data easier to analyse.

You’re correct about using the freqz function. To have freqz correspond with your signal, you would want it to be the length of ‘Fv’ (the length of the fft), not the length of your original signal. I’m at a loss as to your freqz call throwing an error. I called it the same way and it executed without error:

Fn = 50;
[n,Wn] = buttord([0.75 3.5]/Fn, [0.5 4.0]/Fn, 1, 10);   
[z,p,k] = butter(n,Wn);
[sos,g] = zp2sos(z,p,k);
freqz(sos, 1500)

Zero-phase filtering prevents phase distortion in the filtering process. (In the analogue filters I designed for biomedical signal processing, I always used Gaussian filters for the same reason.) You can use the filter function instead of filtfilt, but you will likely see more transients and phase distortion than with filtfilt.

The original baseline is the mean of the original signal, so add that to the filtered signal if you want it to have your signal close to its original baseline. It will include a small contribution from the baseline drift signal, but that is unavoidable. The idea behind the filtering is to make the identification of the times (actually the indices) of the minimum hip angles easier to identify. I would do all the analyses on the filtered signals, and report the baselines separately.

In my experience, physiological data are never negative unless they are derivatives of the actual signal or are in reference to something else. That is the reason the lognormal distribution is generally used to describe them, rather than the normal or other distributions that can take on negative values.

While I’m thinking about it, note that for consistency, you should filter all your data with the same filter, and present the filtered data the same way with the same thresholding and other analysis criteria. Also, report the filter passbands in your ‘Materials and Methods’ section along with the type of filter (Butterworth using second-order-sections implementation), and your sampling frequency.