Spectral Analysis of Water Wave Gauge Data

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!

Answer 1

[Star Strider]

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.

[Wilson Meireles]

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]

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.

.