FFT of tidal data - problems identifying Mf component in MATLAB

I am running an FFT in MATLAB of some tidal data. I successfully pick up most of the major tidal species, however I can't seem to be able to pick up the fortnightly component (Mf), even though it is apparent in the time domain. Please see the link to the image below. Does anyone have any suggestions as to why this is?

[David Young]

Is it possible to attach the actual data, or a link to it?

[Donald John]

Hi, thank you for your response. Attached is the raw data.

[Donald John]

The output from the FFT is also attached.

Answer 1

[David Young]

The apparent fortnightly frequency is in fact a beat pattern - there is little energy at this frequency and the low-frequency noise masks any peak there in the FFT. I think this is the point that Youssef KHMOU is making. [ Edit: I initially thought there was no energy there. In fact, there is a component in the gravitational potential with a fortnightly frequency, but its amplitude is relatively small.]

In fact, your data do show good agreement with the expected peaks. The code below demonstrates it - if you run it you will see that the red tick marks at the twice-daily frequencies line up very nicely with the peaks in your data. Only the largest diurnal component is visible above the noise, but you also have strong thrice-daily and four-times-daily peaks.

Note that it's useful to plot only the low-frequency end of the spectrum, and taking logs makes the lower peaks much easier to see. Also removing the mean helps visibility and accuracy.

Incidentally, the FFT isn't a good way to estimate the amplitude and phase of tidal components, if that's what you aim to do. It's much better to do a direct least-squares fit of waves at the known frequencies.

%%Get data
d = load('Tide_data');
tide = d.NH_tide_data;
t = d.NH_tide_data_time;

%%Take FFT and plot with correct frequencies on x-axis
deltat = mean(diff(t));    % have checked that t has uniform increments
N = length(tide);
deltaf = 1 / (N * deltat);   % frequency step in day^(-1)
f = (0:N-1)*deltaf;          % frequency vector, correct origin and scale
ftide = fft(tide - mean(tide));    % fft

nplot = floor(N/20);    % choose how much to plot
plot(f(2:nplot), log(abs(ftide(2:nplot))));

%%Compare with astronomical tidal periods
% periods in hours of M2 K1 S2 O1 M1 N2 Q1 L2 2N2 J1 OO1 
% (main diurnal and twice-daily components)
tideperiodsHours = [12.421 23.934 12.0 25.819 24.833 12.658 26.868 ...
    12.192 12.905 23.098 22.306];
tidefreqs = 24./tideperiodsHours;    % main frequencies in day^(-1)
hold on
for k = 1:length(tidefreqs)
    freq = tidefreqs(k);           
    ind = round(freq/deltaf) + 1;
    tideval = log(abs(ftide(ind)));  % position tick mark vertically
    plot([freq freq], tideval + [0.1 0.5], 'r-');
end
hold off

[ Edit: additional information about tidal component estimation added.]

The FFT is fundamentally a bad method for estimating tidal component amplitudes, because it fits sinusoids whose periods divide exactly into the time span of the series. These periods won't coincide with the tidal periods, so errors are introduced.

The simplest alternative is to fit a model made of sinusoids at the known tidal frequencies, and this works very well. The least-squares approach is known as "classical harmonic analysis", for example in this PowerPoint presentation by Russ Herman. It has the additional advantage that missing data are not a problem. A basic function that carries out harmonic analysis is attached.

You can check out the results by estimating the amplitudes, then seeing how a resynthesised time series, plotted below in red, compares with the original, plotted in blue:

%%Compare fitted model with original data
% estimate tidal amplitudes
tidenomean = tide - mean(tide);
[a, b] = harmonicFit(t, tidenomean, tidefreqs);
disp(hypot(a, b));    % display amplitudes
% synthesise time series using only these frequencies
synth = harmonicSynth(a, b, tidefreqs, t);
% plot initial section only for clarity
nplot = 500;
plot(t(1:nplot), tidenomean(1:nplot), 'b-', ...
    t(1:nplot), synth(1:nplot), 'r-');

There are visible mismatches. I suspect this is because my array of frequencies does not include the 3-cycles-per-day and higher frequencies which are prominent in the Newhaven data. You can also compare in the frequency domain, by looking at the power spectrum of the residuals, like this:

% Look at spectrum of residuals
fresid = fft(tidenomean.' - synth);
nplot = floor(N/20);
plot(f(2:nplot), log(abs(fresid(2:nplot))));

Adding more frequencies would clearly improve the fit.

However, harmonic analysis is still relatively unsophisticated, and there's a big literature in the area, going back to George Darwin. One way to improve on the simple function attached is to get error estimates, and a good paper dealing with this is here. Pawlowicz has MATLAB code relating to this here. It may also be well worth checking out some of the other toolboxes around, such as this one on the FEX.

Finally, I think it's worth having a look at the nineteenth-century version of my attached file harmonicSynth.m - it's much nicer to look at and is here.

[Donald John]

Thank you for this very comprehensive reply, it is incredibly useful. Would you be able to point me in the direction of an explanation as to the beat pattern effect, I was under the impression that the fortnightly frequency is in fact a significant energy component of the tide. Finally do you have any pointers for the least square method you outline?

[David Young]

I don't have an obvious reference, but textbooks on acoustics, signal processing and the Fourier transform should be helpful. The Wikipedia article on beats is quite helpful. The first equation shows how beats arise from the sum of two sine waves of close frequencies. Since the FT picks apart signals into a sum of sine waves, it will analyse the signal into those two original frequencies (f1 and f2 in the Wikipedia article), not into the beat frequency.

I'll try and get back shortly on fitting sinusoids to time series.

[David Young]

A correction to my answer. You're right that there is a real peak with a fortnightly period in the spectrum of the gravitational potential; as you say, it's Mf. I had forgotten about it. It's quite small - amplitude about 1/6 that of M2.

However, you don't see it in your data because it's lost in the low-frequency noise. The effect that you see in the graph of the time series - the change in the envelope between spring and neap tides - is beating between nearby semidiurnal and diurnal frequencies, mainly between S2 and M2, and is not due to the Mf component.

[Donald John]

Okay, I think that makes sense. Is there other methods that could be applied to pick up the Mf component? For example would the least squares method you mentioned prove useful here?

[David Young]

I think the Mf component is below the level of the low-frequency noise, and you are not likely to get a reliable estimate of its amplitude unless you can get hold of a much longer time series. I'm about to edit the answer to include information about tidal component estimation.

[Donald John]

Thanks for the incredibly useful edit and extension to your answer. I am now trying to use the Lomb Scargle method, which is a least squares method. It is very useful because I have other data that is unevenly sampled. However I am having an issue with a different data set, that has a much shorter time-series. Is there anyway to improve the accuracy of peak estimates using short-time series with these least squares based methods?

[David Young]

The least-squares methods will work correctly with short time series (and the advantage over FT/periodogram-based methods is even greater). To improve the accuracy you need to introduce some constraints on the model - the classical way is to assume that the response varies smoothly with frequency so that for nearby frequencies the amplitudes are in proportion to the amplitudes of the driving potential. This seminal paper for this is Munk and Cartwright (1966); there may be more recent work that is less complex.

[Adi Purwandana]

Excuse me, Since I am a beginner in matlab data files; could you help me to modify "%% Get data" into xlsx file instead of using .mat file (first column is time, second column is tide)? sadly, I use matlab r2013a, so i have to use the syntax "X=dir('blabla.xlsx')" then to call the file.

Many thanks, Adi

[Adi Purwandana]

solved

Answer 2

[Youssef Khmou]

From the result in the second figure; you can remark that there are two frequency components. which means that there are two wave forms with different magnitudes and infinitesimal difference of frequencies, without using the data, the following simulation produces the same pattern ;

   Fs=20;t=0:1/Fs:40-1/Fs;
   y=sin(2*pi*4*t)+0.2*sin(2*pi*4.2*t);
   plot(t,y)
   figure; plot((abs(fft(y))))

basic solution consists of making the region of interest smaller, in the given example it is :

axis([100 200 0 450]);

[Donald John]

Thanks, I agree. These components are the 12.42 hour M2 and the 12 hour, S2, component of the tide. However there is clearly a fortnightly frequency within the data from the top figure, the FFT analysis does not seem to pick this up. Otherwise there would be another frequency spike at about 10^1 days. This is not seen in the FFT results in the lower figure. Any ideas why?

[Youssef Khmou]

how do you recognize the fortnightly frequency in the first graph? try :

spectrogram(y);

which is short term transform,

Answer 3

[Erik S.]

Hi,

Have you tried to plot the power spectral density instead? If you provide the sampling frequency you can get better insight to which frequencies you have.

You get the following graph by pwelch( NH_tide_data-mean(NH_tide_data) );

[David Young]

Yes, but it makes more sense if you give the x-axis the right units - cycles/day in this case. Also, it's not clear that the smoothing implicit in Welch's method is appropriate in the case where there are well-defined sinusoids in the data - the method trades off the frequency resolution against noise, and this may not be desirable here,

Answer 4

[qj fan]

Dear David Young,

 i have a question.    There are some daily recorded series for temperature, and show evidently seasonal trend. how can i get the main frequency and get rid of the periodic trend by fft. one way is the following algorithm:
Smoothing Filter :
1. Let the power-law noise be represented by x and the sinusoidal trend by t. Therefore the noise with the trend is given by y = x + t; with Fourier transform F (f ).
2. Let the frequency in the Fourier transform corresponding to the periodic trend occur at f = fk.
3. Replace the power at the frequency fk by a smoothing filter
|F∗(fk)| = 0.5(|F(fk−1)| + |F(fk+1)|) |F∗(f)| = |F(f)| ,f ̸= fk
Assign random phases so as to satisfy the conjugacy constraints.
4. Determine the inverse Fourier transform to obtain the filtered data yf .

because fft need sample frequency, i don't know how to solve it.

would you like to give me some suggestion.