How can I create isolated waves from data I have?
2 views (last 30 days)
Show older comments
Andrew Mosqueda
on 9 Oct 2022
Commented: Image Analyst
on 10 Oct 2022
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
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
Star Strider
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 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);
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!