How can I create isolated waves from data I have?
Show older comments
I have data from a buoy that measured the wave height every 0.25 seconds. Applying a downward crossing or upward crossing method (everytime the slope goes down or up and crosses 0 the wave starts/begins), how would I isolate the waves to be able to measure their height (crest minus trough) and period?
EDIT: Attached Data File
1 Comment
Andrew Mosqueda
on 9 Oct 2022
Accepted Answer
Star Strider
on 9 Oct 2022
Edited: Star Strider
on 9 Oct 2022
The most robust way I am aware of to detect zero-crossings is:
zci = find(diff(sign(y)))
where ‘y’ is the signal youwant to analyse, the interpolate to get the exact times —
t = linspace(0, 10, 250);
y = sin(2*pi*t);
zci = find(diff(sign(y)));
for k = 1:numel(zci)
idxrng = max(1,zci(k)-1) : min(numel(t),zci(k)+1);
tx(k) = interp1(y(idxrng), t(idxrng), 0);
end
figure
plot(t, y)
hold on
plot(tx, zeros(size(tx)), 'rs')
hold off
grid
That will give you the time of each zero-crossing, and from those you can extract each period.
Use the findpeaks on the positive and negative of the data to get the peaks and troughs respectively, and their times and amplitudes. Other options are islocalmax and islocalmin.
Make appropriate changes to work with your data.
EDIT — Corrected typograhpical errors.
EDIT — (10 Oct 2022 at 18:16)
Analysing provided data —
T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1150580/waves_n_130.dat.txt');
DT = datetime(T1{:,2}, 'ConvertFrom','datenum');
figure
plot(DT, T1{:,4} )
grid
xlabel('Time')
ylabel('Amplitude')
title('Summary Data')
t = T1{:,2};
y = T1{:,4};
zci = find(diff(sign(y)));
for k = 1:numel(zci)
idxrng = max(1,zci(k)-1) : min(numel(t),zci(k)+1);
tx(k) = interp1(y(idxrng), t(idxrng), 0);
end
tx = tx(:);
tx = datetime(tx, 'ConvertFrom','datenum');
t = DT;
[pks, plocs] = findpeaks(y, 'MinPeakProminence',0.01, 'MinPeakDistance',45);
[trs, tlocs] = findpeaks(-y, 'MinPeakProminence',0.01, 'MinPeakDistance',45);
figure
plot(t, y, 'DisplayName','Data')
hold on
plot(tx, zeros(size(tx)), 'rs', 'DisplayName','Zero-Crossings')
plot(t(plocs), pks, '^r', 'DisplayName','Peaks')
plot(t(tlocs), -trs, 'vr', 'DisplayName','Troughs')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
title('Detailed Data')
legend('Location','best')
xlim([min(DT) min(DT)+minutes(2.5)])
% ZeroCrossingTimes = tx
% Periods = diff(tx)
% PeakTimes = t(plocs)
% PeakAmplitudes = pks
% TroughTimes = t(tlocs)
% TroughAmplitudes = -trs
I did not do any filtering of these data, although that could be an appropriate option.
EDIT — (9 Oct 2022 at 20:02)
Adding a filter —
L = numel(t);
Fs = 1/0.25; % Hz
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTy = fft(y.*hanning(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTy(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform')
xlim([0 0.5])
xline(0.09,'--r')
tv = T1{:,2};
y = T1{:,4};
y = lowpass(y, 0.09, Fs, 'ImpulseResponse','iir');
zci = find(diff(sign(y)));
for k = 1:numel(zci)
idxrng = max(1,zci(k)-1) : min(numel(tv),zci(k)+1);
tvx(k) = interp1(y(idxrng), tv(idxrng), 0);
end
tx = tvx(:);
tx = datetime(tx, 'ConvertFrom','datenum');
t = DT;
[pks, plocs] = findpeaks(y, 'MinPeakProminence',0.01, 'MinPeakDistance',45);
[trs, tlocs] = findpeaks(-y, 'MinPeakProminence',0.01, 'MinPeakDistance',45);
figure
plot(t, y, 'DisplayName','Data')
hold on
plot(tx, zeros(size(tx)), 'rs', 'DisplayName','Zero-Crossings')
plot(t(plocs), pks, '^r', 'DisplayName','Peaks')
plot(t(tlocs), -trs, 'vr', 'DisplayName','Troughs')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
title('Detailed Data')
legend('Location','best')
xlim([min(DT) min(DT)+minutes(2.5)])
ZeroCrossingTimes = tx
Periods = diff(tx)
PeakTimes = t(plocs)
PeakAmplitudes = pks
TroughTimes = t(tlocs)
TroughAmplitudes = -trs
.
2 Comments
Andrew Mosqueda
on 10 Oct 2022
Edited: Andrew Mosqueda
on 10 Oct 2022
Sure! Those are quite similar to findpeaks.
Using islocalmin and isolcalmax (linked to earlier and introduced in R2017b) and changing the lowpass filter to a smoothdata call —
T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1150580/waves_n_130.dat.txt');
DT = datetime(T1{:,2}, 'ConvertFrom','datenum', 'Format','dd-MMM-yyyy HH:mm:ss.SSS');
figure
plot(DT, T1{:,4} )
grid
xlabel('Time')
ylabel('Amplitude')
title('Summary Data')
t = T1{:,2};
y = T1{:,4};
L = numel(t);
Fs = 1/0.25; % Hz
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTy = fft(y.*hanning(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTy(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform')
xlim([0 0.5])
xline(0.09,'--r')
tv = T1{:,2};
y = T1{:,4};
% % y = lowpass(y, 0.09, Fs, 'ImpulseResponse','iir');
% y = smoothdata(y, 'sgolay', 55, 'omitnan', 'Degree',3)
y = smoothdata(y, 'rlowess', 55);
zci = find(diff(sign(y)));
for k = 1:numel(zci)
idxrng = max(1,zci(k)-1) : min(numel(tv),zci(k)+1);
tvx(k) = interp1(y(idxrng), tv(idxrng), 0);
end
tx = tvx(:);
tx = datetime(tx, 'ConvertFrom','datenum', 'Format','dd-MMM-yyyy HH:mm:ss.SSS');
t = DT;
% [pks, plocs] = findpeaks(y, 'MinPeakProminence',0.01, 'MinPeakDistance',45);
% [trs, tlocs] = findpeaks(-y, 'MinPeakProminence',0.01, 'MinPeakDistance',45);
LvPks = islocalmax(y, 'MinProminence',0.01, 'MinSeparation',45); % 'islocalmax'
plocs = find(LvPks);
pks = y(plocs);
LvTrf = islocalmin(y, 'MinProminence',0.01, 'MinSeparation',45); % 'islocalmin'
tlocs = find(LvTrf);
trs = y(tlocs);
figure
plot(t, y, 'DisplayName','Data')
hold on
plot(tx, zeros(size(tx)), 'rs', 'DisplayName','Zero-Crossings')
plot(t(plocs), pks, '^r', 'DisplayName','Peaks')
plot(t(tlocs), trs, 'vr', 'DisplayName','Troughs')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
title('Detailed Data')
legend('Location','best')
xlim([min(DT) min(DT)+minutes(2.5)])
ZeroCrossingTimes = tx
Periods = diff(tx) % 'duration' Array — Does Not Allow Milliseconds
PeakTimes = t(plocs)
PeakAmplitudes = pks
TroughTimes = t(tlocs)
TroughAmplitudes = -trs
There are two smoothdata calls, one using 'sgolay' and the other 'rlowess' because the 'sgolay' option might require the Signal Processing Toolbox. The 'rlowess' option does not, and the both give comparable results to the lowpass filter. Change the window (55 here) to get different results.
.
More Answers (1)
Image Analyst
on 9 Oct 2022
Edited: Image Analyst
on 9 Oct 2022
You could use image analysis. Do you have the Image Processing Toolbox? If so, something like this (untested of course because you forgot to attach your data):
meanHeight = mean(signal)
crests = signal >= meanHeight;
troughs = signal < meanHeight;
% Measure max values in the crests segments
propsCrests = regionprops(crests, signal, 'MaxIntensity')
maxCrestLevels = [propsCrests.MaxIntensity]
% Measure min values in the troughs segments
propsTroughs = regionprops(troughs, signal, 'MinIntensity')
minTroughLevels = [propsTroughs.MinIntensity]
% subtract them but make sure you align them first, like maybe you want to
% get the crest height as the crest level minus the average or min of the
% two troughs on either side.
If you have any more questions, then attach your data and code to read it in with the paperclip icon after you read this:
8 Comments
Image Analyst
on 9 Oct 2022
Just put it inside a .zip file or if it's text, change the extension to .txt or .csv.
Andrew Mosqueda
on 9 Oct 2022
Edited: Andrew Mosqueda
on 9 Oct 2022
Image Analyst
on 9 Oct 2022
Edited: Image Analyst
on 9 Oct 2022
Type
>> ver
to see if you have it.
Here's the code with your file:
baseFileName = 'waves_n_130.dat.txt';
fullFileName = fullfile(pwd, baseFileName);
data = readmatrix(fullFileName);
x = data(:, 1);
signal = data(:, end);
subplot(2, 1, 1);
plot(x, signal, 'b-');
grid on;
title('Wave Height');
xlabel('Reading Number (Time Point)')
meanHeight = mean(signal)
crests = signal >= meanHeight;
troughs = signal < meanHeight;
% Measure max values in the crests segments
propsCrests = regionprops(crests, signal, 'MaxIntensity')
maxCrestLevels = [propsCrests.MaxIntensity]
% Measure min values in the troughs segments
propsTroughs = regionprops(troughs, signal, 'MinIntensity')
minTroughLevels = [propsTroughs.MinIntensity]
% subtract them but make sure you align them first, like maybe you want to
% get the crest height as the crest level minus the average or min of the
% two troughs on either side.
crestToTroughDistance = maxCrestLevels - minTroughLevels;
subplot(2, 1, 2);
plot(crestToTroughDistance, 'b-');
grid on;
caption = sprintf('%d Crest-to-Trough Distances', numel(crestToTroughDistance));
title(caption);
xlabel('Crest Number')
% Take a histogram
subplot(3, 1, 3);
histogram(crestToTroughDistance)
grid on;
title('Histogram of Crest-to-Trough Distances')
xlabel('Crest-to-Trough Distance')
ylabel('Count')
Image Analyst
on 9 Oct 2022
What if you've got small crests and troughs on a large crest? How do you define crest? For example here is a wave with small crests and small troughs superimposed on much larger crests and troughs. What vertical distances would you want to measure in this case? Is a crest the whole thing if the height is above the long-term mean? Or do you want to measure the small peaks and valleys even if they're on something that is a large trough?
x = linspace(0, 5*pi, 2000);
% Create slowly varying, large amplitude wave.
period1 = 5;
amplitude1 = 7;
y1 = amplitude1 * sin(2 * pi * x / period1);
plot(x, y1, 'b-')
% Create rapidly varying, small amplitude wave.
period2 = 0.6;
amplitude2 = 2;
y2 = amplitude2 * sin(2 * pi * x / period2);
% Add them together and plot the sum.
yBoth = y1 + y2;
plot(x, yBoth, 'r-', 'LineWidth', 2)
grid on;
Andrew Mosqueda
on 10 Oct 2022
Image Analyst
on 10 Oct 2022
Not sure I understand. So in the plot above are there 4 peaks, 1 peak (highest of the high), or 27 (local peaks)?
Andrew Mosqueda
on 10 Oct 2022
Image Analyst
on 10 Oct 2022
OK, so the code I gave should work except for the case where you have a tiny peak briefly crossing the meanHeight line. In that case you will have tiny, short regions for crests and troughs. You can filter those out with bwareaopen or bwareafilt
meanHeight = mean(signal)
crests = signal >= meanHeight;
% Throw out small crests
minAllowableWidth = 10; % elements, or whatever length you want.
crests = bwareaopen(crests, minAllowableWidth);
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
