version 1.1 (3.37 KB) by
Matt Tearle

Computes coefficients for the real Fourier series approximation to a data set.

Includes two functions: Fseries.m and Fseriesval.m

[a,b] = Fseries(X,Y,n) fits an nth-order Fourier expansion of the form

y = a_0/2 + Sum_k[ a_k cos(kx) + b_k sin(kx) ]

to the data in the vectors X & Y, using a least-squares fit.

Y = Fseriesval(a,b,X) evaluates the Fourier series defined by the coefficients a and b at the values in the vector X.

Extra arguments allow for rescaling of X data and sin-only or cosine-only expansions.

Example:

% Generate data

x = linspace(0,2,41)';

y = mod(2*x,1);

% Use FSERIES to fit

[a,b,yfit] = Fseries(x,y,10);

% Evaluate on finer grid

xfine = linspace(0,2)';

yfine = Fseriesval(a,b,xfine);

% Visualize results

plot(x,y,'x',x,yfit,'o',xfine,yfine)

This generates the attached image of a 10-term Fourier series approximation of a sawtooth wave.

Matt Tearle (2021). Simple real Fourier series approximation (https://www.mathworks.com/matlabcentral/fileexchange/31013-simple-real-fourier-series-approximation), MATLAB Central File Exchange. Retrieved .

Created with
R2010b

Compatible with any release

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Paul SafierThank you for making this.

oai daothanhDear everyone,

is it true that:

a=[a0 a1 a2 ... ak ...bn] and b=[b1 b2 ... bk ... bn]

and how can we can get frequency of 1-harmonic?

Thank you

Matt Tearle@Zhao Jin: See my reply to Amir further down. I suspect it has to do with the Curve Fitting app using the frequency as one of the fitted parameters. That means it's fitting a fundamentally different equation to a Fourier Series. From my experiments, if you call Fseries with the scaling option set to false, and run the Curve Fitting app with w forced to 1 (you can set bounds on the parameters with "Fit Options"), you get the same coefficient values.

Zhao JinI have a series of input named X1 and Y1, but after using the 'curve fitting',I can get some coefficient. However they are not the same as what i get from your Fseries.m ,in other words, they are different totally. I have no idea how the matter takes place?

Justin BrooksMatt Tearle@Alessandra Scarton: technically yes, but it's extremely messy unless x = 2*pi*(0:n-1)/n. FFTs don't care about the x values; they work only on frequency numbers. Also, I'm doing a least squares fit to a truncated expansion. So I don't know of any simple way to interpret the coefficients from my function in terms of the FFT values. (That's actually what motivated me in the first place - I got sick of trying to go via an FFT when all I really wanted was a simple sin/cos representation of a function.)

Alessandra ScartonIf I apply this to a signal, is there a connection between [a,b] and the result I get from a fft?

Matt Tearle@JaeHwang Jung: correct. a should have one more element than b. a(1) = a_0, a(2) = a_1, a(3) = a_2, etc, while b(1) = b_1, b(2) = b_2, etc.

JaeHwang Jungso a(1) is a_0 which is constant, and a(2) is a coefficient of sin(x), right? (in your example)

Matt TearleHi everyone. This has been updated to deal with the issue in 19b that broke the input parsing. It should now work for any version, but please let me know if you encounter further problems (I made a quick hack solution to get it back up and running...)

Habib AghasafariNot compatible with newer versions of Matlab like R2019b.

Matt Tearle@Fangzhou Wang: Thanks for letting me know. There was a change in R2019b to the precedence rules that happened to affect a bunch of my files (because of an old-school way I was doing some input parsing/checking). I didn't realize this one was busted, too. I'll work on updating it. In the meantime, it should work in any version up to 19a.

Fangzhou WangHello Matt,

thank you for sharing the code. I don't know why but the expamle in the description is not excutable. i got the error massage as following:

"Error:File: Fseries.m Line: 63 Column: 4

Identifier 'xrow' is not a function or a shared variable. To share 'xrow' with nested function, initialize it in the current scope. For more information,

see Sharing Variables Between Parent and Nested Functions."

Thank you very much if you could help me to solve this.

best regards

William Simmons@Matt - thank you for your quick response! This is extremely helpful information, and the function works perfectly.

Matt Tearle@William: I don't know of any deep, subtle mathematical reasons why you can't/shouldn't do that. If you treat it as a least-squares fit, where some of the x values happen to be repeated, it will work. One minor consideration, though, is that every data point is equally weighted, so you should really have the same number of y values for each x value (although that probably won't cause significant issues unless they're seriously imbalanced.) With this function, you'll need to "unroll" everything into one x vector (with repeated values) and one y vector - it won't allow a vector for x and a matrix for y.

William SimmonsThank you for this function!

I have a question - can you use a Fourier Series to fit a curve to a plot with multiple Y-values for each X-value? For example, if multiple temperature readings were taken per day throughout a region for many consecutive days, then plotted (X=day, Y=temperature), how would this function fit a Fourier curve to it? I have been able to fit a curve to analogous data using this function, and it looks reasonable, but I want to make sure it's statistically valid. Thanks!

Jie LIUMatt Tearle@Helena: Do you have a specific example function you were trying it on? When I try it on the sawtooth function in the example (x = linspace(0,2,41)'; y = mod(2*x,1);), I get the same values as the Curve Fitting app. Note that I use a different convention for the constant term than the CF app uses. My functions use y = a0/2 + ... whereas the CF app uses y = a0 + ... So there's a factor of 2 difference in that one term. The other values should be identical, though.

helena melzerThanks for this code. But I have a question too. I tried to get the same values with the curve fitting app. The coefficients for b are the same, but the values for a are different. Can you explain why?

I have set the Fseries with the scaling option to false, and run the Curve Fitting app with w forced to 1 (as you explained before)

Matt Tearle@Alexander: There's no weighting -- this is just pure Fourier series. It sounds like maybe you're encountering Gibb's phenomenon. If you fit a function f(x) on some interval [0 L] using Fourier series, there's an implicit assumption that f(x) is L-periodic -- ie f(x) repeats on [L 2L], [2L 3L], etc. Unless f(0) = f(L), that means you're introducing a discontinuity at the end points of the interval, which leads to Gibb's oscillations near the ends.

To your earlier question, I wonder if maybe an optimization/least-squares fitting approach might be better for your problem. The joins between the different periods could be treated as a constraint. If you can share some of the details, I'd recommend posting to MATLAB Answers (https://www.mathworks.com/matlabcentral/answers/). There you can show what your data looks like and the form of your model that you're trying to fit. Someone there will probably have useful suggestions about the best way to tackle the problem.

Alexander Graf@Matt Tearle: I saw sometimes that the fitting quality is very bad at the end of the data vector while the Fourier series fits very good with the points "in the middle"? Do you know why? Is there any kind of weightening included?

Alexander Graf@Matt Tearle: Thanks for your answer! Meanwhile everything is clear.

Your function works very well, nice Job!

Currently, i have a Problem that I don't how to solve. Maybe you can help.

My signal that i get from the experimental setup is sinusodial and has of course some errors. Not all of the periods of this Sinus function exhibit the same errors (different frequencies and amplitude). Since the errors are systematic and originate from the experimental Setup, they belong therefore always to the same period of the measurement signal. Therefore I thought about dividing the error function into periods (similar to the periods of the input signal) and fit each error period seperately with a Fourier Series. That works quite well, but a challenging problem is, that the remaining error curve isn't continuous any more. This could be a problem in the further proceeding. Maybe it will be not critical if the "jumps" between last point of the prevenient period and the first point of the consecutive period are not largely spaced.

Do you have an idea to solve this issue?

Matt Tearle@Alexander: there's no upper limit on n, but I wouldn't recommend this for your application, by the sounds of it. I think anything over a few thousand terms in the series is going to bring this function to a grinding halt. You may be better off using the FFT.

[It also depends on what the *range* of frequencies is. I'm assuming you mean a range from O(1) up to O(1e5). That's going to be difficult for this function. If it's just high frequencies (in Hz), then rescaling may help.]

Alexander GrafIs it possible to fit some data with an unknown and quite high frequency (around 10^5)?

Matt Tearle@Chaoyang Wu: Thanks for your feedback. I haven't thought through the statistical theory to know whether this is really valid, but naively you could calculate the coefficient of determination from the outputs, using the formula R^2 = 1 - (residual sum of squares)/(total sum of squares): >> R2 = 1 - sum((yfit - y).^2)/sum((y-mean(y)).^2)

Regarding your question about the "extra w", I assume you're referring to the fit that the Curve Fitting app is doing. See my reply to Amir Zakaria below for more details, but basically the Curve Fitting app allows the fundamental frequency to be a parameter. (This is then a nonlinear fitting problem.) Hope that helps.

CHAOYANG WUThank you.Since it uses least square method.Could you provide a value such as R square to evaluate the fitting?In matlab,there is extra w item.Could anyone tell me the difference, please?

goc3Inbar FrenkelPaul BarrattNeilTeresa LoThank you so much!!!!!!!!!!!!

Giorgio FagioliPerfect code and very useful since already built-in functions of matlab do not allow fitting beyond eight terms. Thank you!

Zuowei ZhuVery useful tool, thanks for your work!

Matt Tearle@Amir Zakaria: Thanks for the feedback. I just had a look at what the Curve Fitting app is doing at its "Fourier" option includes the fundamental frequency as one of the fit parameters. So it's fitting a_k cos(w*k*x), where the coefficients a_k *and* the frequency w are parameters. My function is intended for just plain Fourier series expansion (a_k cos(k*x)). If you call Fseries with the scaling option set to false, and run the Curve Fitting app with w forced to 1 (you can set bounds on the parameters with "Fit Options"), you get the same values. Hope that helps.

Amir ZakariaI found this very useful, but when i compare the same number of coefficients for example 5 using this function and using cftool, i have different values. can someone explain why?

Zachary Hugodingvery good tool

MoisésChihanKevin HammondsDoes anyone know if there is a way to extract a_o , a_n, and b_n from the command line or the script in the more general cos and sin terms?

TillmannReally great!

But, how can I visualize more than one period?

Miguel AngelI have a question: I've used your code to fit a Fourier series to a set of data (t,x). However, when I apply the same found coefficients to a different vector of time t', I do not obtain the same function (I expected that were the case). What can be occuring? Thank you so much in advance.

Miguel AngelVery fast and useful.

Thank you very much.

Miguel AngelChristiaan Tolvery helpful, thanks! I used it for my project

Matt TearleThe expansion is in terms of sin/cos(kx) for k = 0:n, so the frequencies are simply k/L (for k from 0 to n), where L = max(x) - min(x). The units would be in inverse units of x. In the example given in the description, you can see the dominant contribution from the 4th sine term, which corresponds to a frequency of 4/(2-0) = 2 Hz (if x represented time in seconds).

Matt TearleThanks for the feedback, too. Glad it could help.

Jong-Hwan KimCould you explain how to produce frequency from your code?

Jong-Hwan Kimthis is a very good tool. I really like it. It helps my project. Thanks