Spectral Analysis of Water Wave Gauge Data
49 views (last 30 days)
Show older comments
Wilson Meireles
on 18 Oct 2023
Commented: Star Strider
on 18 Oct 2023
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 Comments
Accepted Answer
Star Strider
on 18 Oct 2023
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 Comments
Star Strider
on 18 Oct 2023
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')
T1 = readtable(Uz{1}, 'VariableNamingRule','preserve')
VN = T1.Properties.VariableNames;
T1.Time = datetime(T1.Time, 'InputFormat','mm:ss.S', 'Format','mm:ss.S')
T1.DeltaTime = (0:size(T1,1)-1).'*0.25
% T1.DeltaTime.Format = 'mm:ss.S'
ix = find(abs(T1.Time - max(T1.Time)) > 0.9)
Dix = diff(ix)
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)
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.
.
More Answers (0)
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!