Spectral Analysis of Water Wave Gauge Data

Question

Wilson Meireles 2023-10-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2035369-spectral-analysis-of-water-wave-gauge-data

评论： Star Strider 2023-10-18

Hello all! I have a 2 week long data set of wave data collected in the field and I am looking to conduct a spectral analysis. I am struggling a bit with the initial plotting of a spectrum of this data and I am also looking to run a filter to remove tide data. I have tried running a high pass butterworth filter and using the fft function but have not had success. Any help would be appreciated!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-10-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2035369-spectral-analysis-of-water-wave-gauge-data#answer_1336009

It would help to have the data.

Be certain that the independent variable data are regularly sampled, so that the sampling frequency is constant. (The nufft function can take non-uniformly sampled data, however it is the exception. All other signal processing functions assume a uniform sampling frequency.) Determine this by calculating the mean and standard deviation (std) of diff(x) where ‘x’ is the independent variable vector. The standard deviation should be several (5 or more) orders of magnitude less than the mean if the data are regularly sampled. If it is not regularly sampled, use the resample function to resample it to a uniform sampling frequency.

Design the filter by first using fft (or pspectrum) to see the frequency characteristics of the signal. Choose the passband frequency on that basis. In calculating the fft, subtract the mean of the signal from the signal in order to see the peaks clearly. Use the second-order-section (zp2sos) filter implementation for best results and the most stable filter. You can design your own filter, however it might be easier to use the highpass function with the 'ImpulseResponse','iir' name-value pair for best results. (That’s just easier.)

After that, experiment to get the result you want.

2 个评论
显示无隐藏无

Wilson Meireles 2023-10-18

LH data test 3.zip

Thank you very much for the help! I have attached a portion of the data here as well, it is just raw pressure data that I have no problem converting to wave heights and doing a wave by wave analysis of, just struggling with the spectral portion.

Star Strider 2023-10-18

在 MATLAB Online 中打开

LH data test 3.zip

The ‘Time’ variable made this a bit of a challenge. I have no idea what the ‘Time’ format is, or the units (it also ‘wraps’ at one hour), so I created something that looks reasonably representative as ‘DeltaTime’ and am using that.

Try this —

Uz = unzip('LH data test 3.zip')

Uz = 1×1 cell array

{'LH data test 3.csv'}

T1 = readtable(Uz{1}, 'VariableNamingRule','preserve')

T1 = 328823×4 table

Time Pressure Sea pressure Depth ___________ ________ ____________ _______ {'02:34.0'} 10.275 0.1424 0.14124 {'02:34.0'} 10.275 0.14238 0.14121 {'02:34.5'} 10.276 0.14327 0.1421 {'02:34.7'} 10.276 0.14301 0.14185 {'02:35.0'} 10.275 0.14226 0.1411 {'02:35.2'} 10.274 0.1419 0.14074 {'02:35.5'} 10.275 0.14277 0.1416 {'02:35.8'} 10.279 0.14606 0.14486 {'02:36.0'} 10.272 0.13987 0.13873 {'02:36.2'} 10.271 0.13884 0.13771 {'02:36.5'} 10.273 0.14079 0.13964 {'02:36.7'} 10.271 0.13859 0.13746 {'02:37.0'} 10.271 0.13823 0.13711 {'02:37.2'} 10.271 0.13815 0.13702 {'02:37.5'} 10.271 0.13801 0.13688 {'02:37.7'} 10.271 0.13898 0.13785

VN = T1.Properties.VariableNames;

T1.Time = datetime(T1.Time, 'InputFormat','mm:ss.S', 'Format','mm:ss.S')

T1 = 328823×4 table

Time Pressure Sea pressure Depth _______ ________ ____________ _______ 02:34.0 10.275 0.1424 0.14124 02:34.0 10.275 0.14238 0.14121 02:34.5 10.276 0.14327 0.1421 02:34.7 10.276 0.14301 0.14185 02:35.0 10.275 0.14226 0.1411 02:35.2 10.274 0.1419 0.14074 02:35.5 10.275 0.14277 0.1416 02:35.8 10.279 0.14606 0.14486 02:36.0 10.272 0.13987 0.13873 02:36.2 10.271 0.13884 0.13771 02:36.5 10.273 0.14079 0.13964 02:36.7 10.271 0.13859 0.13746 02:37.0 10.271 0.13823 0.13711 02:37.2 10.271 0.13815 0.13702 02:37.5 10.271 0.13801 0.13688 02:37.7 10.271 0.13898 0.13785

T1.DeltaTime = (0:size(T1,1)-1).'*0.25

T1 = 328823×5 table

Time Pressure Sea pressure Depth DeltaTime _______ ________ ____________ _______ _________ 02:34.0 10.275 0.1424 0.14124 0 02:34.0 10.275 0.14238 0.14121 0.25 02:34.5 10.276 0.14327 0.1421 0.5 02:34.7 10.276 0.14301 0.14185 0.75 02:35.0 10.275 0.14226 0.1411 1 02:35.2 10.274 0.1419 0.14074 1.25 02:35.5 10.275 0.14277 0.1416 1.5 02:35.8 10.279 0.14606 0.14486 1.75 02:36.0 10.272 0.13987 0.13873 2 02:36.2 10.271 0.13884 0.13771 2.25 02:36.5 10.273 0.14079 0.13964 2.5 02:36.7 10.271 0.13859 0.13746 2.75 02:37.0 10.271 0.13823 0.13711 3 02:37.2 10.271 0.13815 0.13702 3.25 02:37.5 10.271 0.13801 0.13688 3.5 02:37.7 10.271 0.13898 0.13785 3.75

% T1.DeltaTime.Format = 'mm:ss.S'

ix = find(abs(T1.Time - max(T1.Time)) > 0.9)

ix = 0×1 empty double column vector

Dix = diff(ix)

Dix = 0×1 empty double column vector

figure

plot(T1.Time,'.-')

% figure

% plot(T1.DeltaTime, T1.('Sea pressure'))

% xlabel('\Delta Time (s?)')

% ylabel('Amplitude (Sea pressure Units)')

% title('Highpass-Filtered Sea pressure')

[pks,locs] = findpeaks(T1.('Sea pressure'), 'NPeaks',1, 'MinPeakHeight',7)

pks = 8.4417

locs = 23575

figure

plot(T1.DeltaTime, T1.('Sea pressure'))

hold on

plot(T1.DeltaTime(locs), T1.('Sea pressure')(locs),'^r')

hold off

grid

xlabel('\Delta Time (s?)')

ylabel('Amplitude (Sea pressure Units)')

title('Original Sea pressure')

t = T1.DeltaTime(locs:end);

Sp = T1.('Sea pressure')(locs:end);

Fs = 1/(t(2)-t(1));

Fn = Fs/2;

figure

plot(t, Sp)

grid

xlabel('\Delta Time (s?)')

ylabel('Amplitude (Sea pressure Units)')

title('Unfiltered Sea pressure')

L = numel(t);

NFFT = 2^nextpow2(L);

FTSp = fft((Sp-mean(Sp)).*hann(L), NFFT)/L;

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

figure

plot(Fv, abs(FTSp(Iv))*2)

grid

xlim([0 0.0005])

xline(0.00005, '--', 'Cutoff Frequency')

xlabel('Frequency (Hz?)')

ylabel('Magnitude (Sea pressure Units)')

SpHPF = highpass(Sp, 0.0005, Fs, 'ImpulseResponse','iir');

figure

plot(t, SpHPF)

grid

xlabel('\Delta Time (s?)')

ylabel('Amplitude (Sea pressure Units)')

title('Highpass-Filtered Sea pressure')

It would help to have the time variable explained, however that would only require some changes in the units, and likely not the rest of the code. A bandpass filter could minimise some of the high-frequency noise if you want to experiment with that. Use the same frequency as I did here for the lower passband, and experiment with the upper passband.

.

请先登录，再进行评论。