determine location of certain turning points of curve
9 次查看(过去 30 天)
显示 更早的评论
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! :)
6 个评论
Star Strider
2014-11-2
编辑:Star Strider
2014-11-2
Thank you. That will help if I need to design a filter.
回答(1 个)
Star Strider
2014-11-2
编辑:Star Strider
2014-11-2
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)')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/176467/image.png)
6 个评论
Star Strider
2014-11-4
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Digital Filter Design 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!