How can I create isolated waves from data I have?
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Andrew Mosqueda
am 9 Okt. 2022
Kommentiert: Image Analyst
am 10 Okt. 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
Akzeptierte Antwort
Star Strider
am 9 Okt. 2022
Bearbeitet: Star Strider
am 9 Okt. 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 Kommentare
Star Strider
am 10 Okt. 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.
.
Weitere Antworten (1)
Image Analyst
am 9 Okt. 2022
Bearbeitet: Image Analyst
am 9 Okt. 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 Kommentare
Image Analyst
am 10 Okt. 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);
Siehe auch
Kategorien
Mehr zu Spectral Measurements finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!