How can I create isolated waves from data I have?

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

[Andrew Mosqueda]

The total rows for number of times measured is 8191 (2047.5 seconds in total) if it helps

Answer 1

[Star Strider]

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

ZeroCrossingTimes = 258×1 datetime array
   23-Dec-2021 03:00:00
   23-Dec-2021 03:00:09
   23-Dec-2021 03:00:17
   23-Dec-2021 03:00:25
   23-Dec-2021 03:00:33
   23-Dec-2021 03:00:41
   23-Dec-2021 03:00:49
   23-Dec-2021 03:00:57
   23-Dec-2021 03:01:05
   23-Dec-2021 03:01:12
   23-Dec-2021 03:01:22
   23-Dec-2021 03:01:30
   23-Dec-2021 03:01:41
   23-Dec-2021 03:01:49
   23-Dec-2021 03:01:57
   23-Dec-2021 03:02:04
   23-Dec-2021 03:02:13
   23-Dec-2021 03:02:22
   23-Dec-2021 03:02:30
   23-Dec-2021 03:02:38
   23-Dec-2021 03:02:47
   23-Dec-2021 03:02:55
   23-Dec-2021 03:03:03
   23-Dec-2021 03:03:13
   23-Dec-2021 03:03:20
   23-Dec-2021 03:03:33
   23-Dec-2021 03:03:43
   23-Dec-2021 03:03:51
   23-Dec-2021 03:03:59
   23-Dec-2021 03:04:07

Periods = diff(tx)

Periods = 257×1 duration array
   00:00:08
   00:00:08
   00:00:08
   00:00:07
   00:00:08
   00:00:08
   00:00:07
   00:00:08
   00:00:07
   00:00:09
   00:00:08
   00:00:10
   00:00:08
   00:00:07
   00:00:07
   00:00:08
   00:00:08
   00:00:08
   00:00:07
   00:00:09
   00:00:08
   00:00:07
   00:00:10
   00:00:06
   00:00:12
   00:00:10
   00:00:08
   00:00:07
   00:00:07
   00:00:08

PeakTimes = t(plocs)

PeakTimes = 123×1 datetime array
   23-Dec-2021 03:00:05
   23-Dec-2021 03:00:21
   23-Dec-2021 03:00:37
   23-Dec-2021 03:00:53
   23-Dec-2021 03:01:08
   23-Dec-2021 03:01:26
   23-Dec-2021 03:01:45
   23-Dec-2021 03:02:00
   23-Dec-2021 03:02:18
   23-Dec-2021 03:02:33
   23-Dec-2021 03:02:51
   23-Dec-2021 03:03:06
   23-Dec-2021 03:03:30
   23-Dec-2021 03:03:48
   23-Dec-2021 03:04:03
   23-Dec-2021 03:04:20
   23-Dec-2021 03:04:37
   23-Dec-2021 03:04:56
   23-Dec-2021 03:05:11
   23-Dec-2021 03:05:26
   23-Dec-2021 03:05:44
   23-Dec-2021 03:05:59
   23-Dec-2021 03:06:16
   23-Dec-2021 03:06:35
   23-Dec-2021 03:06:51
   23-Dec-2021 03:07:09
   23-Dec-2021 03:07:25
   23-Dec-2021 03:07:42
   23-Dec-2021 03:07:59
   23-Dec-2021 03:08:22

PeakAmplitudes = pks

PeakAmplitudes = 123×1

0.2370 0.2493 0.2663 0.3552 0.1826 0.1672 0.2007 0.3074 0.3644 0.3105

TroughTimes = t(tlocs)

TroughTimes = 122×1 datetime array
   23-Dec-2021 03:00:13
   23-Dec-2021 03:00:29
   23-Dec-2021 03:00:46
   23-Dec-2021 03:01:00
   23-Dec-2021 03:01:18
   23-Dec-2021 03:01:37
   23-Dec-2021 03:01:53
   23-Dec-2021 03:02:09
   23-Dec-2021 03:02:26
   23-Dec-2021 03:02:42
   23-Dec-2021 03:02:59
   23-Dec-2021 03:03:17
   23-Dec-2021 03:03:36
   23-Dec-2021 03:03:55
   23-Dec-2021 03:04:10
   23-Dec-2021 03:04:29
   23-Dec-2021 03:04:43
   23-Dec-2021 03:05:04
   23-Dec-2021 03:05:19
   23-Dec-2021 03:05:35
   23-Dec-2021 03:05:52
   23-Dec-2021 03:06:08
   23-Dec-2021 03:06:25
   23-Dec-2021 03:06:43
   23-Dec-2021 03:06:59
   23-Dec-2021 03:07:17
   23-Dec-2021 03:07:32
   23-Dec-2021 03:07:50
   23-Dec-2021 03:08:06
   23-Dec-2021 03:08:29

TroughAmplitudes = -trs 

TroughAmplitudes = 122×1

-0.2801 -0.2842 -0.2995 -0.2769 -0.1457 -0.1330 -0.2806 -0.3341 -0.4409 -0.2659

.

[Andrew Mosqueda]

Your stuff has helped me get the periods down, but I cannot use find peaks because I don't have the signal processing toolbox. Can you show me how I'd do this with islocalmax/min?

[Star Strider]

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

ZeroCrossingTimes = 246×1 datetime array
   23-Dec-2021 03:00:09.117
   23-Dec-2021 03:00:17.480
   23-Dec-2021 03:00:25.324
   23-Dec-2021 03:00:33.402
   23-Dec-2021 03:00:41.503
   23-Dec-2021 03:00:50.991
   23-Dec-2021 03:00:56.940
   23-Dec-2021 03:01:05.595
   23-Dec-2021 03:01:13.325
   23-Dec-2021 03:01:22.513
   23-Dec-2021 03:01:31.514
   23-Dec-2021 03:01:40.688
   23-Dec-2021 03:01:49.402
   23-Dec-2021 03:01:57.140
   23-Dec-2021 03:02:05.343
   23-Dec-2021 03:02:13.693
   23-Dec-2021 03:02:22.766
   23-Dec-2021 03:02:30.619
   23-Dec-2021 03:02:38.249
   23-Dec-2021 03:02:47.364
   23-Dec-2021 03:02:55.285
   23-Dec-2021 03:03:03.809
   23-Dec-2021 03:03:15.272
   23-Dec-2021 03:03:21.545
   23-Dec-2021 03:03:33.534
   23-Dec-2021 03:03:41.705
   23-Dec-2021 03:03:51.250
   23-Dec-2021 03:03:59.525
   23-Dec-2021 03:04:07.183
   23-Dec-2021 03:04:15.787

Periods = diff(tx)                      % 'duration' Array — Does Not Allow Milliseconds

Periods = 245×1 duration array
   00:00:08
   00:00:07
   00:00:08
   00:00:08
   00:00:09
   00:00:05
   00:00:08
   00:00:07
   00:00:09
   00:00:09
   00:00:09
   00:00:08
   00:00:07
   00:00:08
   00:00:08
   00:00:09
   00:00:07
   00:00:07
   00:00:09
   00:00:07
   00:00:08
   00:00:11
   00:00:06
   00:00:11
   00:00:08
   00:00:09
   00:00:08
   00:00:07
   00:00:08
   00:00:08

PeakTimes = t(plocs)

PeakTimes = 123×1 datetime array
   23-Dec-2021 03:00:05.097
   23-Dec-2021 03:00:21.600
   23-Dec-2021 03:00:37.843
   23-Dec-2021 03:00:55.123
   23-Dec-2021 03:01:08.342
   23-Dec-2021 03:01:26.918
   23-Dec-2021 03:01:45.408
   23-Dec-2021 03:02:01.651
   23-Dec-2021 03:02:18.585
   23-Dec-2021 03:02:34.656
   23-Dec-2021 03:02:51.590
   23-Dec-2021 03:03:08.352
   23-Dec-2021 03:03:28.396
   23-Dec-2021 03:03:47.836
   23-Dec-2021 03:04:03.388
   23-Dec-2021 03:04:20.150
   23-Dec-2021 03:04:37.344
   23-Dec-2021 03:04:55.401
   23-Dec-2021 03:05:12.336
   23-Dec-2021 03:05:26.851
   23-Dec-2021 03:05:43.872
   23-Dec-2021 03:05:59.596
   23-Dec-2021 03:06:16.617
   23-Dec-2021 03:06:35.107
   23-Dec-2021 03:06:51.868
   23-Dec-2021 03:07:08.889
   23-Dec-2021 03:07:25.910
   23-Dec-2021 03:07:42.153
   23-Dec-2021 03:07:59.347
   23-Dec-2021 03:08:21.897

PeakAmplitudes = pks

PeakAmplitudes = 123×1

0.1118 0.1416 0.1384 0.0610 0.0924 0.0909 0.1265 0.1545 0.2139 0.1654

TroughTimes = t(tlocs)

TroughTimes = 122×1 datetime array
   23-Dec-2021 03:00:13.651
   23-Dec-2021 03:00:29.894
   23-Dec-2021 03:00:46.396
   23-Dec-2021 03:01:00.912
   23-Dec-2021 03:01:17.414
   23-Dec-2021 03:01:36.336
   23-Dec-2021 03:01:53.875
   23-Dec-2021 03:02:10.377
   23-Dec-2021 03:02:27.657
   23-Dec-2021 03:02:43.123
   23-Dec-2021 03:02:59.366
   23-Dec-2021 03:03:17.596
   23-Dec-2021 03:03:36.604
   23-Dec-2021 03:03:55.872
   23-Dec-2021 03:04:11.856
   23-Dec-2021 03:04:28.617
   23-Dec-2021 03:04:45.897
   23-Dec-2021 03:05:04.387
   23-Dec-2021 03:05:20.112
   23-Dec-2021 03:05:35.145
   23-Dec-2021 03:05:50.611
   23-Dec-2021 03:06:08.409
   23-Dec-2021 03:06:25.344
   23-Dec-2021 03:06:43.660
   23-Dec-2021 03:06:59.904
   23-Dec-2021 03:07:17.875
   23-Dec-2021 03:07:34.118
   23-Dec-2021 03:07:50.880
   23-Dec-2021 03:08:07.123
   23-Dec-2021 03:08:29.846

TroughAmplitudes = -trs 

TroughAmplitudes = 122×1

0.1552 0.1614 0.1823 0.1449 0.0989 0.0832 0.1183 0.1066 0.1844 0.1800

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.

.

Answer 2

[Image Analyst]

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:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[Image Analyst]

Just put it inside a .zip file or if it's text, change the extension to .txt or .csv.

[Andrew Mosqueda]

I don't have image processing toolbox, at least I dont think so. I'm a novice still at MATLab and don't know where it'd be.

[Image Analyst]

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]

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]

I was able to get back to my code and I don't have image procesing tools. The crests and troughs are the local maximum and minimums of the wave; so only lowest of the low and highest of the high. Also the data has already been "normalized" to where the amplitudes are the distance from the mean sea level (y=0) where the data was recorded.

[Image Analyst]

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]

Yeah, I'm sorry, reading back it's kind of contradictory, but it's 4 peaks in that plot.

[Image Analyst]

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);