Consider real signal of the form:
$$
\begin{align*} x(n) = &\sum_{m=1}^{M} \left( c_m \cos(2\pi f_m n) + s_m \sin(2\pi f_m n) \right) =\\
&\sum_{m=1}^{M} a_m \cos(2\pi f_m n + \phi_m)
\end{align*}
$$
where signal values of $x(n)$ are known for $n = 1,\dots, N >> M$ and some frequencies $f_m$ are known real numbers in $(0.0,0.5)$.

What is the best way to estimate $c_m$ and $s_m$ real parameters?

Estimation of the form:
$$c_m + j s_m = 2 \frac{\sum_{n=1}^{N} x(n) \exp(2 \pi j f_m n)}{N}$$
is not very accurate ($j^2=-1$), but windowing (e.g. cosine window) improves estimation accuracy significantly:
$$c_m + j s_m = 2 \frac{\sum_{n=1}^{N} w(n) x(n) \exp(2 \pi j f_m n)}{\sum_{n=1}^{N} w(n)}$$

Why window is so effective here? Is it some kind of well-known formula?

The least-squares method also gives good results if frequencies are known for all terms.

Are there any other good methods for the case when not all frequencies are known? Perhaps based on SVD or some kind of orthogonalization.

Test code in Mathematica:
```Mathematica
F1 = 0.123000 ;
F2 = 0.312874 ;

LENGTH = 1024 ;
RANGE = Range[LENGTH] ;

SIGNAL = 0.0 ;
SIGNAL = SIGNAL + 1.0*Cos[2*Pi*F1*RANGE] + 0.5*Sin[2*Pi*F1*RANGE] ;
SIGNAL = SIGNAL + 0.05*Cos[2*Pi*F2*RANGE] + 0.01*Sin[2*Pi*F2*RANGE] ;

WINDOW = HannWindow[N[Rescale[Range[LENGTH],{1,LENGTH},{-1,1}]]]^2 ;

(* NO WINDOW *)
2*(Total[SIGNAL*Exp[2*Pi*I*F1*RANGE]]/LENGTH)               
2*(Total[SIGNAL*Exp[2*Pi*I*F2*RANGE]]/LENGTH)               
(* OUTPUT: 0.9994929536863852`\[VeryThinSpace]+0.49994639046951617` I *)
(* OUTPUT: 0.04819835064062154`\[VeryThinSpace]+0.012127390241918076` I *)

(* WITH WINDOW *)
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F1*RANGE]]/Total[WINDOW]) 
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F2*RANGE]]/Total[WINDOW])
(* OUTPUT: 0.9999999999775656`\[VeryThinSpace]+0.5000000000025541` I *)
(* OUTPUT: 0.05000000004559643`\[VeryThinSpace]+0.009999999958900167` I *)

(* LEAST SQUARES *)
(* ONE KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE]}//Transpose,SIGNAL]
(* ALL KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE],Cos[2*Pi*F2*RANGE],Sin[2*Pi*F2*RANGE]}//Transpose,SIGNAL]
(* OUTPUT: {0.999957445314386`,0.4999499273288035`} *)
(* OUTPUT: {1.`,0.49999999999999983`,0.05000000000000022`,0.010000000000000056`} *)
```