Making something like 'findpeaks' function without using threshold

18 次查看(过去 30 天)
Hi! i want to ask about the code i am making. i am making myself a for loop like a finpeaks (but i don't want to use that and any threshold for a more automatic approach) on my data and i am having a problem with showing the valid peaks/valleys because some ot the shown peaks/valleys were not the one i want to be shown. I have my code here and was cracking my brains out for a solution. If you look at the photo, the plot cannot show some valid peaks/valleys (colored greed), and some non valid ones are shown in some. Can you help me?
I attached the data and the photo of the plot.
% Design for bandpass filter (Butterworth filter)
fs = 1000; % Sampling frequency
f_low = 0.1; % Low cutoff frequency
f_high = 10; % High cutoff frequency
order = 2; % Filter order
% Bandpass filtering
[b, a] = butter(order, [f_low, f_high]/(fs/2));
mag_filtered = filtfilt(b, a, acc);
% Set up plot with aspect ratio adjustment
figure;
% Plot for the main data
h_main = plot(NaN, NaN, 'k'); % Main plot in black
hold on;
h_peaks = plot(NaN, NaN, 'or'); % Red markers for peaks
h_valleys = plot(NaN, NaN, 'ob'); % Blue markers for valleys
xlabel('Time (s)');
ylabel('Filtered Magnitude');
title('Animated Peaks and Valleys');
% Set up additional plot for the entire data
h_total = plot(NaN, NaN);
% Set up additional plot for the animated entire data
h_total_animated = plot(NaN, NaN, 'k');
% Set up a new figure for temporary data
figure;
h_temp = plot(NaN, NaN);
xlabel('Time (s)');
ylabel('Filtered Magnitude');
title('Temporary Data');
% Animation parameters
step = 300; % Number of points/data to display in each step
delay = 0.01; % Delay between steps in seconds
% Initialize arrays to store peaks and valleys
peaks_x = [];
peaks_y = [];
valleys_x = [];
valleys_y = [];
% Initialize additional variables
time_total = [];
mag_total = [];
time_data_new = [];
mag_data_new = [];
% Initialize variables to track highest peak and lowest valley
highest_peak = -inf;
lowest_valley = inf;
% Initialize variables to track all peaks and valleys
all_peaks_x = [];
all_peaks_y = [];
all_valleys_x = [];
all_valleys_y = [];
% Initialize variables to track the number of steps since the last peak or valley
steps_since_peak = 0;
steps_since_valley = 0;
% Main loop for animation
for i = 1:step:numel(time)
% Extract a portion of data for the current step
time_data_new = time(i:min(i+step-1, numel(time)));
mag_data_new = mag_filtered(i:min(i+step-1, numel(time)));
% Find the highest peak and lowest valley within the current step
max_value = max(mag_data_new);
min_value = min(mag_data_new);
% Update the highest peak and lowest valley if found
if max_value > highest_peak
highest_peak = max_value;
peaks_x = time_data_new(mag_data_new == highest_peak);
peaks_y = highest_peak;
steps_since_peak = 0;
else
steps_since_peak = steps_since_peak + 1;
end
if min_value < lowest_valley
lowest_valley = min_value;
valleys_x = time_data_new(mag_data_new == lowest_valley);
valleys_y = lowest_valley;
steps_since_valley = 0;
else
steps_since_valley = steps_since_valley + 1;
end
% Update all peaks and valleys arrays
all_peaks_x = [all_peaks_x peaks_x];
all_peaks_y = [all_peaks_y peaks_y];
all_valleys_x = [all_valleys_x valleys_x];
all_valleys_y = [all_valleys_y valleys_y];
% Update plot for the entire data
set(h_total, 'XData', time(1:i), 'YData', mag_filtered(1:i));
% Update plot for animated entire data
set(h_total_animated, 'XData', time(1:i), 'YData', mag_filtered(1:i));
% Update main plot
set(h_main, 'XData', time_data_new, 'YData', mag_data_new);
% Update plot
drawnow; % Force the plot to refresh
pause(delay);
% Plot all stored peaks and valleys
set(h_peaks, 'XData', all_peaks_x, 'YData', all_peaks_y);
set(h_valleys, 'XData', all_valleys_x, 'YData', all_valleys_y);
% Update plot for time_temp and mag_temp
set(h_temp, 'XData', time_data_new(:), 'YData', mag_data_new(:));
drawnow;
pause(delay);
% Reset the highest peak and lowest valley if no new peak or valley is found within 10 steps
if steps_since_peak >= 10
highest_peak = -inf;
end
if steps_since_valley >= 10
lowest_valley = inf;
end
end
  2 个评论
Steven Lord
Steven Lord 2024-3-14
If one of the reasons you don't want to use findpeaks is because you don't have Signal Processing Toolbox, take a look at the islocalmax and islocalmin functions in MATLAB.
Stella
Stella 2024-3-15
Hi! I have the signal processing toolbox, but my PI doesn't want to use it because we'll use different data on it and and there were instances that when I used it, he saw that I used different MinPeakDistance values on the data. I also used islocalmin/max, but still not good enough for him 😅

请先登录,再进行评论。

采纳的回答

Kevin Holly
Kevin Holly 2024-3-14
Is there a particular reason that you don't want to use findpeaks function?
To be clear on what you would like as an output, would the peaks and valleys detected below work for you?
%Load data
Data = importfile('TV_500ml_Hong.csv');
figure('Color','w')
plot(Data.TimePlot1,Data.YgPlot1,'Color','k')
title('Animated Peaks and Valleys')
xlabel('Time (s)')
ylabel('Filtered Magnitude')
prominencethreshold = 0.01;
[peaks, locsPeaks] = findpeaks(Data.YgPlot1, 'MinPeakProminence', prominencethreshold);
hold on; % Keep the plot for adding peaks and valleys
plot(Data.TimePlot1(locsPeaks), peaks, 'v', 'MarkerFaceColor', 'g','MarkerEdgeColor','none'); % Plot peaks
% Find and plot valleys by inverting the data
[invertedValleys, locsValleys] = findpeaks(-Data.YgPlot1, 'MinPeakProminence', prominencethreshold);
valleys = -invertedValleys; % Convert back to original scale
plot(Data.TimePlot1(locsValleys), valleys, '^', 'MarkerFaceColor', 'r','MarkerEdgeColor','none'); % Plot valleys
hold off;
legend('Data', 'Peaks', 'Valleys','Location','northeastoutside');
In the next section, I applied the peak and valley detection to the other 3 plots found in the spreadsheet. A helper function called plotpeaksandvalleys was created to reduce repetitive code.
figure
tiledlayout(4,1)
nexttile
plot(Data.TimePlot0,Data.XgPlot0,'Color','k')
title('Animated Peaks and Valleys')
xlabel('Time (s)')
ylabel('Magnitude')
hold on
plotpeaksandvalleys(Data.TimePlot0,Data.XgPlot0,prominencethreshold)
nexttile
plot(Data.TimePlot1,Data.YgPlot1,'Color','k')
title('Animated Peaks and Valleys')
xlabel('Time (s)')
ylabel('Magnitude')
hold on
plotpeaksandvalleys(Data.TimePlot1,Data.YgPlot1,prominencethreshold)
nexttile
plot(Data.TimePlot2,Data.ZgPlot2,'Color','k')
title('Animated Peaks and Valleys')
xlabel('Time (s)')
ylabel('Magnitude')
hold on
plotpeaksandvalleys(Data.TimePlot2,Data.ZgPlot2,prominencethreshold)
nexttile
plot(Data.TimePlot3,Data.TotalgPlot3,'Color','k')
title('Animated Peaks and Valleys')
xlabel('Time (s)')
ylabel('Magnitude')
hold on
plotpeaksandvalleys(Data.TimePlot3,Data.TotalgPlot3,prominencethreshold)
function plotpeaksandvalleys(x,y,prominencethreshold)
[peaks, locsPeaks] = findpeaks(y, 'MinPeakProminence', prominencethreshold);
hold on; % Keep the plot for adding peaks and valleys
plot(x(locsPeaks), peaks, 'v', 'MarkerFaceColor', 'g','MarkerEdgeColor','none'); % Plot peaks
% Find and plot valleys by inverting the data
[invertedValleys, locsValleys] = findpeaks(-y, 'MinPeakProminence', prominencethreshold);
valleys = -invertedValleys; % Convert back to original scale
plot(x(locsValleys), valleys, '^', 'MarkerFaceColor', 'r','MarkerEdgeColor','none'); % Plot valleys
hold off;
legend('Data', 'Peaks', 'Valleys','Location','northeastoutside');
end
  1 个评论
Stella
Stella 2024-3-15
Hi! Thank you for your answer! Before i was using the findpeaks function in my code, but my PI don't want to use it because he wants to make an automated algorithm that for example every 50% that the data was being updated per step (300 data) on the plot, it will look for the valid peaks/valleys and will continue to do so for the rest of the data. i did try it in the code above, but the peaks/valleys are not really good. I also did the custom function you did, but i don't know why he doesn't like to have that too so i sticked with the if-loop in the for-loop.

请先登录,再进行评论。

更多回答(0 个)

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by