The source of confusion might be that with FFT we need to reorder the positive and negative frequency bins with fftshift().
Please check the code below where complex coefficients are used.
using FFTW, LaTeXStrings, Plots; gr()
function FourierSeries(x, fx, T)
ck = 1/N * fftshift(fft(fx))
n = 1 + length(x)÷2
yr = zero(fx) .+ 0.0im
for (i, ci) in pairs(ck)
yr .+= ci * exp.(im*2π*(i-n)/T * x)
end
return real(fftshift(yr)), ck
end
N = 100 # use N even
x = LinRange(-π, π, N+1)[1:N]
T = 2π
fx = @. sin(2π*x/T) + cos(2π*4*x/T)
fsx, ck = FourierSeries(x, fx, T)
p1 = plot(xlabel=L"x \in [-\pi,\pi]", legend=:topleft)
plot!(x, fx, lw=3, label="f(x)")
plot!(x, fsx, lw=1, lc=:black, ls=:dash, label="Fourier Series")
xp = 2π*(-N÷2:N÷2-1)/T
p2 = plot(xp, real.(ck), st=:stem, lc=:blue, label="Real(ck)")
plot!(xp, imag.(ck), st=:stem, lc=:red, xlabel="k", label="Imag(ck)")
plot(p1,p2)
