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.
d = load('Tide_data');
tide = d.NH_tide_data;
t = d.NH_tide_data_time;
deltat = mean(diff(t));
N = length(tide);
deltaf = 1 / (N * deltat);
f = (0:N-1)*deltaf;
ftide = fft(tide - mean(tide));
nplot = floor(N/20);
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;
for k = 1:length(tidefreqs)
freq = tidefreqs(k);
ind = round(freq/deltaf) + 1;
tideval = log(abs(ftide(ind)));
plot([freq freq], tideval + [0.1 0.5], 'r-');
[ 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:
tidenomean = tide - mean(tide);
[a, b] = harmonicFit(t, tidenomean, tidefreqs);
synth = harmonicSynth(a, b, tidefreqs, t);
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:
fresid = fft(tidenomean.' - synth);
nplot = floor(N/20);
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.