There are three posible methods to achieve this
1. Create a custom equation in the Curve Fitting Toolbox defining as many coefficients as you like.
You can refer to this documentation page for more info.
2. Compute the FFT of the signal and modify it as appropriate.
If you have N data points in your signal, then the FFT essentailly gives you the discrete Fourier coefficients for the first N/2 frequencies. Using all of these coefficients would allow you to reconstruct a curve that passes through all of your data exactly. If you wanted a filtered fit, then you could set some of the terms in the FFT corresponding to higher order frequencies equal to zero.
3. Use a submission on the File Exchange.
There is also a MathWorks File Exchange submission for computing Fourier series approximations in a very straightforward manner:
This function directly returns the coefficients of the series (to find them from the FFT you would have to do some manipulation of the complex-valued vector of FFT data using ABS and ANGLE), and it computes the fit based on least squares rather than simply setting higher frequencies equal to zero and using the lower order terms of the FFT.
Below is some sample code comparing the methods 2 and 3. To run the code first download the file Fseries.m from the File Exchange link above. Note that The MathWorks does not guarantee or warrant the use or content of these submissions. Any questions, issues, or complaints should be directed to the contributing author.
N = 128;
t = 1:N;
fftx = fft(x);
maxHar = 20;
if maxHar > ceil(N/2) + 1
maxHar = ceil(N/2) + 1;
fftxMod = fftx;
fftxMod(maxHar+2:end-maxHar) = 0;
[afit bfit yfit] = Fseries(t,x,maxHar);
legend('original data',['filtered FFT to ',num2str(maxHar),' frequencies'],['Fit with Fseries to ',num2str(maxHar),' frequencies'])