- 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=

# Peaks and Valleys of Noisy, Increasing Sinewave

49 views (last 30 days)

Show older comments

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;

##### 0 Comments

### Answers (2)

Shubham
on 15 Oct 2024 at 6:17

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:

I hope this helps!

##### 0 Comments

William Rose
on 15 Oct 2024 at 6:55

Edited: William Rose
on 15 Oct 2024 at 7:12

@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.

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!