Fourier decomposition in Matlab

Question

I am doing a Fourier decomposition of sea level series ( to find out one particular contribution within the sea level height). it should be applied subsequently on shifted (by one hour) windows of 96 hours length. I have already the $a$, $b$ as the following:

$$a_{i,n}=\frac{\displaystyle\sum_{j=0}^{95} x_{i+j}\cos^2\left(2\pi jT^{-1}_n\right)}{\displaystyle\sum_{j=0}^{95} \cos^2\left(2\pi jT^{-1}_n\right)}$$

$$b_{i,n}=\frac{\displaystyle\sum_{j=0}^{95} x_{i+j}\sin^2\left(2\pi jT^{-1}_n\right)}{\displaystyle\sum_{j=0}^{95} \sin^2\left(2\pi jT^{-1}_n\right)}$$

where $x$ is the sea level series (hourly dataset)

As a start I tried to code only the first time window as the following code to simplify the coding because I am just a beginner in Matlab but unfortunately I am getting wrong results as the amplitude is higher than the actual sea level here is my code so far

hour_window = 96 ;
T= 19,8 ;
f = 1 / T ;
% sea_level is a vector of the first 96 values ( as I am doing only the first interval)
for j = 0: hour_window -1 
   a = sealevel .* cos(2pifj) ./ cos(2pifj).^2 ;
   b = sealevel .* sin(2pifj) ./ sin(2pifj).^2 ;
end
A = sqrt(a.^2 + b^.2) ;
phase= atan2(a,b) ;