Hi Bartosz,
The functions you used for computation of the fourier coefficient is erroneous. The formula for discrete time fourier coefficients is as follows :
Furthermore, the functions defined in the code for the coefficients of the fourier series are erroneous. f(1) corresponds to the coefficient a0, not coefficient a(1) and therefore, the index of either the function or the sinusoidals needs to be changed appropriately as shown in the code snippet below:
san = san + (f(j)*cos((2*pi/T)*(i-1)*(j-1)));
sbn = sbn - (f(j)*sin((2*pi/T)*(i-1)*(j-1)));
sfn = sfn + (fan(jfn)*cos((jfn-1)*(2*pi/T)*(ifn-1))-fbn(jfn)*sin((jfn-1)*(2*pi/T)*(ifn-1)));
Presented below is the output that was obtained following the aforementioned changes:
Hope this helps