File Exchange

image thumbnail

Simple real Fourier series approximation

version 1.1 (3.37 KB) by Matt Tearle
Computes coefficients for the real Fourier series approximation to a data set.

27 Downloads

Updated 03 Mar 2020

View Version History

View License

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.

Cite As

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 .

Comments and Ratings (49)

Paul Safier

Thank you for making this.

oai daothanh

Dear 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 Jin

I 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 Brooks

Matt 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 Scarton

If 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 Jung

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

Matt Tearle

Hi 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 Aghasafari

Not 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 Wang

Hello 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 Simmons

Thank 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 LIU

Matt 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 melzer

Thanks 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 Graf

Is 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 WU

Thank 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?

goc3

Inbar Frenkel

Paul Barratt

Neil Fu

Teresa Lo

Thank you so much!!!!!!!!!!!!

Giorgio Fagioli

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

Zuowei Zhu

Very 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 Zakaria

I 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 Hugo

ding

very good tool

Moisés

Chihan

Kevin Hammonds

Does 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?

Tillmann

Really great!
But, how can I visualize more than one period?

Miguel Angel

I 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 Angel

Very fast and useful.
Thank you very much.

Miguel Angel

Christiaan Tol

very helpful, thanks! I used it for my project

Matt Tearle

The 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 Tearle

Thanks for the feedback, too. Glad it could help.

Jong-Hwan Kim

Could you explain how to produce frequency from your code?

Jong-Hwan Kim

this is a very good tool. I really like it. It helps my project. Thanks

MATLAB Release Compatibility
Created with R2010b
Compatible with any release
Platform Compatibility
Windows macOS Linux
Tags Add Tags

Community Treasure Hunt

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

Start Hunting!