Efficiently compute the Fourier series coefficients of a piecewise trigonometric function

Question

I am searching for the most efficient algorithm to compute the Fourier series coefficients of a number of given functions $f(\theta)$ $$a_k = \int_0^{2\pi} f(\theta)\cos(k\theta) \\ b_k = \int_0^{2\pi} f(\theta)\sin(k\theta)$$ for $k \in {1,...,N}$. The functions $f(\theta)$ are continuous piecewise trigonometric functions, i.e. in each interval $l$ they can be described as follows: $$f(\theta) = c_{1,l}\cos(\omega\theta) + c_{2,l}\sin(\omega\theta)$$ The functions $f(\theta)$ are also known to be continuous and periodic at $2\pi$. Note that the value of $\omega$ is the same for all sub-intervals.

The goal is to compute $a_k,b_k$ with fastest speed having known the values of $c_{1,l},c_{2,l}$ and $\omega$. Note that there are many number of $f(\theta)$ functions which require great computational efficiency. I know that there is closed form solution for this, but it going to be a trigonometric function itself which at high number of evaluations it becomes costly. One approach that came to mind was using FFT, but then it will require evaluation of $f(\theta)$ at a high number of discrete points which itself requires a lot of trigonometric evaluations. For now, I use the closed form formula. For a 100 number of different $f(\theta)$s each having 100 number of intervals for $N=100$, it takes approximately 1 second. I need to reduce the time to an order of milliseconds if possible.