Peaks and Valleys of Noisy, Increasing Sinewave

Hello!

I am analyzing data from a cyclic test in which axial loading is applied from an Instron actuator to a screwed in material. The data recorded is the amount of loosening of that material. Attached is an image of the raw data and the data file. I only have the loosening data as time was not recorded. I am trying to find the peaks and valleys of the data and take a difference of the peaks and valleys. The idea is that when the difference is plotted, I can fit a trendline and get the rate of loosening from the data. However, I am having difficulties finding the peaks/valleys due to the nature of the data.

I am currently trying to use the findpeaks function for both peaks/valleys and adjusting the MinPeakProminence, MinPeakWidth, and MinPeakHeight. I am also trying to segment the sinewaves to help find the peaks/valleys better. However, it is not working very well. Some of the peaks are recorded while others aren't I have gotten it realtively close but then I have a different number of peaks vs. valleys which causes and error when trying to find the difference. I would do it my hand but I have over 50 tests to anaylze and would like to batch them eventually.

Any and all advice is appreicated! Thank you :)

%% Load Data file
% Load the data (replace with your file path if needed)
orignal_decoup_data = load('dataData.txt');

%% Remove Any Outliers from Official Data - Based on Z-scores

% Compute z-scores for each data point
z_scores = zscore(orignal_decoup_data);

% Define a z-score threshold (e.g., 4 standard deviations from the mean)
threshold = 4;

%Remove outliers based on z-score
decoup_data = orignal_decoup_data(abs(z_scores) <= threshold);
cleaned_decoup_angles = decoup_data; % done for preserving decoup_data

%% Data Filtering

% Filter parameters
order = 3; % Filter order

% Design a low-pass Butterworth filter
[b, a] = butter(order, 0.2, 'low'); 

% Apply the filter
full_signal = filtfilt(b, a, cleaned_decoup_angles);

%% Plot Original data, removed outlier data, and filtered data

subplot(3,1,1)
plot(orignal_decoup_data);
title(['Original Data']);
xlabel('Sample Number');
ylabel('Amplitude');

subplot(3,1,2)
plot(cleaned_decoup_angles);
title('Data After Removing Outliers (Z-Score Method)');
xlabel('Sample Number');
ylabel('Amplitude');
hold on;

subplot(3,1,3)
plot(full_signal);
title('Data with Butterworth Filtered');
xlabel('Sample Number');
ylabel('Amplitude');
hold off;

%% Find initital threshold value based on first few cycles

% find a decent threshold to use the MinPeakHeight option
minPeakHeight = max(findpeaks(full_signal(1:500)))*0.5;

% Makes sure the start is a good start
[int_pk, int_loc] = findpeaks(full_signal);
above_threshold_idx = find(int_pk > minPeakHeight, 1, 'first'); % Find the first peak above the threshold
int_pk_min_loc = abs(int_loc(above_threshold_idx) - 50);

% initial_full_signal = full_signal(1:int_pk_min_loc,1);
initial_full_signal = full_signal(1:500);

thresh_per = 60; % Will take 60% of the max peak foudn int he first few cycles. Can be adjusted if needed.

[int_pk, int_loc] = findpeaks(initial_full_signal);
[int_pk_max, int_max_idx] = max(int_pk);
int_pk_max_loc = int_loc(int_max_idx);

threshold = int_pk_max*(thresh_per/100);  % Amplitude threshold to distinguish between flat and sine waves

% Plot found Threshold Value
figure(2)
plot(full_signal,'color','blue')
hold on;
plot(initial_full_signal,'color','red', LineStyle='--')
hold on;
plot(int_pk_max_loc,int_pk_max,'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'g')
hold off;
legend('Full Signal', 'Initial Signal',['Threshold = ', num2str(round(int_pk_max,3)),' * ',num2str(thresh_per),'%'],'location','best')
title('Find threshold')
xlabel('Sample Number');
ylabel('Rotation Degree')

%% Identify and Segment Sine Waves

% Detecting the start and end points of sine waves
sine_wave_indices = find(abs(full_signal) > threshold);  % Indices where the signal exceeds the threshold

% Identify sections
sections = [];
start_idx = sine_wave_indices(1); % Start from the first detected sine wave index
segement_number = 1;

for i = 2:length(sine_wave_indices)
    % Check if current index is not continuous with the previous index
    if sine_wave_indices(i) > sine_wave_indices(i-1) + 1
        end_idx = sine_wave_indices(i-1);  % End index of the previous segment
        sections = [sections; start_idx, end_idx, segement_number];  % Save the start and end indices
        start_idx = sine_wave_indices(i);  % Update start index
        segement_number = segement_number + 1;
    end
end

% Capture the last section
if ~isempty(start_idx)
    sections = [sections; start_idx, sine_wave_indices(end),segement_number];
end

%% Find Peaks and Valleys for Each Detected Section

% Initialize a cell array to store peaks and valleys for each section
peaks_valleys = cell(size(sections, 1), 1);
peaks_store = [];
valleys_store = [];
peaks_location_store = [];
valleys_location_store = [];

for i = 1:size(sections, 1)-1
    section_signal = full_signal(sections(i, 1):sections(i+1, 2));  % Extract section

    [pks, loc_pks] = findpeaks(section_signal,'MinPeakHeight', threshold);  % Find peaks
    warning('off')
    loc_pks = loc_pks + sections(i, 1)-1;
    [valls, loc_val] = findpeaks(-section_signal,'MinPeakHeight', threshold);  % Find peaks
    warning('off')
    loc_val = loc_val + sections(i, 1)-1;
    warning('off')

    % Find the maximum peak and its location
    [max_peak, max_idx] = max(pks);  % Get the maximum peak value and its index
    max_peak_loc = loc_pks(max_idx);          % Get the location of the maximum peak

    % Find the minimum peak and its location
    [min_valley, min_idx] = max(valls);  % Get the maximum peak value and its index
    min_valley = -min_valley;
    min_valley_loc = loc_val(min_idx);          % Get the location of the maximum peak

    % Live plot the peaks/valleys found
figure(3)
    plot((sections(i, 1):sections(i+1, 2)),full_signal(sections(i, 1):sections(i+1, 2)),'linewidth',1)
    hold on;
    plot(min_valley_loc,min_valley,'ko','MarkerSize', 8,'MarkerFaceColor','auto')
    hold on;
    plot(max_peak_loc,max_peak,'ko','MarkerSize', 8,'MarkerFaceColor','auto')

    peaks_store = [peaks_store;max_peak];
    valleys_store = [valleys_store;min_valley];
    peaks_location_store = [peaks_location_store; max_peak_loc];
    valleys_location_store = [valleys_location_store; min_valley_loc];

    j = i +2;
    i=j;
end

% For some reason the above loop is recording it twice. So this just
% takes the first instance (which happens to be every odd index in the
% data)

peaks_store = peaks_store(1:2:end);
valleys_store = valleys_store(1:2:end);
peaks_location_store = peaks_location_store(1:2:end);
valleys_location_store = valleys_location_store(1:2:end);
% difference_peaks_valleys = peaks_store - valleys_store;

% Plot the Peaks and Valleys with Area
hold on;
plot(full_signal, 'LineWidth', 1.5, 'Color', 'b', 'DisplayName', 'Original Signal');  % Original signal

for i = 1:size(sections, 1)
    % Highlight the sine wave sections in red
    area(sections(i, 1):sections(i, 2), full_signal(sections(i, 1):sections(i, 2)), 'FaceColor', [1 0.5 0.5], 'EdgeColor', 'none', 'FaceAlpha', 0.5, 'DisplayName', 'Sine Wave Section');
end

title('Sine Wave Sections with Peaks and Valleys');
xlabel('Sample Number');
ylabel('Amplitude');
hold off;

Answer 1

[Shubham]

Hey @kjb,

I went through your script and I could not find anything wrong with it. However I have found a query similar to yours on MATLAB Answers forum. I would suggest you to check out the following answer: https://www.mathworks.com/matlabcentral/answers/1946024-find-peaks-and-valley-of-sinusoidal-curve#comment_2781914

I believe you should give a try to the following File Exchange submissions and check if they help you in resolving the issue:

• PeakSeek- generally faster than findpeaks: https://www.mathworks.com/matlabcentral/fileexchange/26581-peakseek?s_tid=srchtitle
• peakfinder- finds local peaks/valleys in a noisy signal: https://www.mathworks.com/matlabcentral/fileexchange/25500-peakfinder-x0-sel-thresh-extrema-includeendpoints-interpolate?requestedDomain=

I hope this helps!

Answer 2

[William Rose]

@kjb,

[Edit: fixed spelling error and incorrect high frequency cutoff value in my description. Code is not changed.]

orig_data = load('dataData.txt');
z_scores = zscore(orig_data);       % z-scores for each data point
threshold = 4;
data = orig_data(abs(z_scores) <= threshold); % remove outliers
N=length(data);

The period of the sinusoidal stimuli is about 20 samples, i.e. frequency ~= 0.05 samples^-1, i.e. w=0.10, where 1=Nyquist freq. Therefore we choose a high pass filter with w=0.07.

[b,a]=butter(2,0.07,'high');
y=filtfilt(b,a,data);

Plot the unfiltered and high-pass-filtered signals, for the full duration, and for short segments from the start, middle, and end.

subplot(211), plot(1:N,data,'-r',1:N,y,'-b'); grid on;
subplot(234), plot(1:N,data,'-r',1:N,y,'-b'); grid on; xlim([250 350])
subplot(235), plot(1:N,data,'-r',1:N,y,'-b'); grid on; xlim([1400 1500])
subplot(236), plot(1:N,data,'-r',1:N,y,'-b'); grid on; xlim([2720 2820])

Find the maxima

[pks1,locs1]=findpeaks(y,'MinPeakProminence',0.1,'MinPeakHeight',0.1,'MinPeakDistance',10);
Np=length(pks1);
figure
subplot(211), plot(1:N,y,'-b',locs1,pks1,'g^'); grid on
subplot(234), plot(1:N,y,'-b',locs1,pks1,'g^'); grid on; xlim([250 350]);
subplot(235), plot(1:N,y,'-b',locs1,pks1,'g^'); grid on; xlim([1400 1500]);
subplot(236), plot(1:N,y,'-b',locs1,pks1,'g^'); grid on; xlim([2720 2820]);

The minima are trickier to find, at least early in the recording, because they are not much lower than the noisy baseline. It will suffice to find the minimum between each of the maxima already found, for maxima separated by less than 30 samples.

pks2=[]; locs2=[]; delta=[];    % arrays for minima, deltas, locations
j=0;                            % counter for minima and deltas
for i=2:Np
    if locs1(i)-locs1(i-1)<30   % if distance between peaks is <30 samples
        j=j+1;                  % increment counter
        [pks2(j),idx]=min(y(locs1(i-1):locs1(i)));  % jth minimum
        locs2(j)=idx+locs1(i-1)-1;                  % location of minimum
        % delta=difference between the mean of the 2 maxes and the minimum in between
        delta(j)=(pks1(i-1)+pks1(i))/2-pks2(j);     % delta
    end
end
figure;
subplot(211), plot(1:N,y,'-b',locs1,pks1,'g^',locs2,pks2,'vr'); grid on
subplot(234), plot(1:N,y,'-b',locs1,pks1,'g^',locs2,pks2,'vr'); grid on; xlim([250 350]);
subplot(235), plot(1:N,y,'-b',locs1,pks1,'g^',locs2,pks2,'vr'); grid on; xlim([1400 1500]);
subplot(236), plot(1:N,y,'-b',locs1,pks1,'g^',locs2,pks2,'vr'); grid on; xlim([2720 2820]);

Plot delta versus sample number. If the amplitude of the applied force is constant, and delta is rising with time, then the system is getting looser or more compliant.

figure
plot(locs2,delta,'k*')
grid on; xlabel('Time (samples)'); ylabel('Delta')

The plot above shows the difference between the minimum and maximum during sinusoidal stresses. If the sinusoidal stress amplitude is constant, then the gradual increase in delta with time indicates the system is becoming looser or more compliant.