How can I create isolated waves from data I have?
4 次查看（过去 30 天）
显示 更早的评论
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
采纳的回答
Star Strider
2022-10-9
编辑：Star Strider
2022-10-9
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 个评论
Star Strider
2022-10-10
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.
.
更多回答（1 个）
Image Analyst
2022-10-9
编辑：Image Analyst
2022-10-9
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 个评论
Image Analyst
2022-10-10
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);
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
您也可以从以下列表中选择网站：
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)