Spectral Analysis of Water Wave Gauge Data
49 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Wilson Meireles
el 18 de Oct. de 2023
Comentada: Star Strider
el 18 de Oct. de 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 comentarios
Respuesta aceptada
Star Strider
el 18 de Oct. de 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 comentarios
Star Strider
el 18 de Oct. de 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.
.
Más respuestas (0)
Ver también
Categorías
Más información sobre Spectral Measurements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!