Making something like 'findpeaks' function without using threshold

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

[Steven Lord]

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]

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 😅

Answer 1

[Kevin Holly]

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

[Stella]

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.