Yes, to calculate by two series.

I am sorry, but I didn't get your point. Do you mean to say that I should split f(x) into two functions f1(x)  and f2(x), such that, 

      f(x) = f1(x) + f2(x)

No, it needs to define "the coordinates" a_n and b_n of function f in common basis u_1, v_1, u_2, v_2,.... They would be the corresponding integrals: a_n=\int{f,u_n}, b_n=\int{f,v_n}. After that it needs to consider the sum \Sigma_n=1^n=N{a_nu_n+b_nv_n} and N->\infty.

The expansion that you have added in your original posting show that your basis functions are actually sine or cosine functions, and \omega and \xi are the frequencies.  In any case, assuming that the basis functions are, say {\phi_n} and {\zeta_n}, then it is possible to expand f(x) as f(x)= \sum_{n=1}^{n=\infty} a_n \phi_n(x) + \sum_{n=1}^{n=\infty} b_n \zeta_n(x) as long as \phi_n and \zeta_m are orthogonal.

There are many different sets of ortogonal functions known as special functions. You could expand the function in terms of any number of these.

However in another one it would no longer be called the Fourier series.

Actually the bases are not required to be orthogonal. One can even use nonorthogonal bases in the series or transform. If you can accept an iterative computation there is also the possibility to have polar coordinates in the series. Essentially one can use almost any function as a basis and even change the basis at each frequency. Of course the orthogonal bases provide the fastest convergence. See:

https://arxiv.org/abs/0804.4347