Hi Ziyi,
It is my understanding that you have tried using the "findpeaks()" function, but it detects peaks in the flat parts as well, which is not desired.
Additionally, when you tried using a threshold, it ignored some peaks in the oscillation part.
Here is an example where you can use the "findpeaks()" function to focus on just the oscillatory portions of the waveforms for computing the phase order:
%% Generate example data (for complete demonstration)
t = linspace(0, 10, 1000);
initialPeak = 1; % Amplitude of the initial peak
frequency = 3; % Frequency of the oscillations (in Hz)
amplitude = 0.2; % Amplitude of the regular oscillations
irregularAmplitude1 = sin(2*pi*0.005*(1:100) + pi/4); % Amplitude of the irregular oscillations during the initial peak for waveform 1
irregularAmplitude2 = sin(2*pi*0.005*(1:100)); % Amplitude of the irregular oscillations during the initial peak for waveform 2
% Generate waveform 1 with an initial peak followed by oscillations
waveform1 = [initialPeak * irregularAmplitude1 .* ones(1, 100), amplitude * sin(2*pi*frequency*t(101:end))];
% Generate waveform 2 with an initial peak followed by oscillations with a different phase
phaseShift = pi/4; % Phase shift between waveform 1 and waveform 2
waveform2 = [initialPeak * irregularAmplitude2 .* ones(1, 100), amplitude * sin(2*pi*frequency*t(101:end) + phaseShift)];
% Introduce a decaying nature in the oscillation portions
decayFactor = 0.99; % Decay factor for the oscillation portions
decay1 = decayFactor .^ (1:length(waveform1)-100);
decay2 = decayFactor .^ (1:length(waveform2)-100);
waveform1(101:end) = waveform1(101:end) .* decay1;
waveform2(101:end) = waveform2(101:end) .* decay2;
%% Considering only the oscillation portions
oscillation1 = waveform1(101:end);
oscillation2 = waveform2(101:end);
[~, locs1] = findpeaks(oscillation1, 'MinPeakProminence', 0.05); % Adjust the threshold as per your data
[~, locs2] = findpeaks(oscillation2, 'MinPeakProminence', 0.05); % Adjust the threshold as per your data
% Plot the waveforms
subplot(2, 1, 1);
plot(t, waveform1);
hold on;
scatter(t(locs1+100), waveform1(locs1+100), 'r');
hold off;
xlabel('Time');
ylabel('Amplitude');
title('Waveform 1');
grid on;
subplot(2, 1, 2);
plot(t, waveform2);
hold on;
scatter(t(locs2+100), waveform2(locs2+100), 'r');
hold off;
xlabel('Time');
ylabel('Amplitude');
title('Waveform 2');
grid on;
As you can see, using the "MinPeakProminence" argument inside the function "findpeaks()", we can focus on to the portion of the waveforms where the oscillations are prominant. Finally, with this knowledge, you can apply your algorithm to determine the phase order.
Additionally, you can refer to the MathWorks documentation of the function "findpeaks()" for better understanding of different use-cases in the link:
I hope this helps!
