Fourier Series Curve Fitting and giving its coefficient with respect to eqaution

% Given x and y coordinates
x = [0.015404238	0.027389034	0.03937383	0.051358625	0.063343421	0.075328217	0.087313012	0.099297808	0.111282604	0.1232674	0.135252195	0.147236991	0.159221787	0.171206582	0.183191378	0.195176174	0.20716097	0.219145765	0.231130561	0.243115357	0.255100152	0.267084948	0.278320694	0.28880739	0.298919562	0.308657208	0.318394855	0.329256076	0.341240872	0.353225667	0.365210463	0.377195259	0.389180054	0.40116485	0.413149646	0.425134442	0.437119237	0.449104033	0.461088829	0.473073624	0.48505842	0.497043216	0.509028012	0.521012807	0.532997603	0.544982399	0.556967194	0.56895199	0.580936786	0.592921582	0.6];  % x coordinates
y = [0.271524849	0.264952735	0.272505755	0.283825478	0.294533112	0.304063651	0.311663755	0.317992597	0.323473931	0.328484427	0.332647415	0.337846246	0.342809658	0.348526411	0.357021107	0.365421635	0.36892545	0.366355457	0.358088326	0.351139541	0.353325011	0.366245582	0.388184835	0.412362705	0.43601428	0.459267874	0.483100962	0.506751175	0.524474293	0.533957748	0.535295709	0.530277359	0.519938542	0.503149247	0.487537046	0.480682429	0.477877017	0.475872029	0.472266193	0.466965341	0.462229494	0.459235747	0.453746559	0.44288982	0.429349306	0.415385037	0.405564142	0.400828295	0.398634972	0.396629985	0.393772585];  % y coordinates

% Define time range
t = 0:0.001:1;

% Number of Fourier coefficients to compute
numCoefficients = 5;

% Compute Fourier coefficients
coefficients = zeros(1, numCoefficients);
for k = 1:numCoefficients
    coefficients(k) = sum(y .* exp(-1i * 2 * pi * (k-1) * x));
    coefficients(k) = coefficients(k) / length(x);
end

% Generate Fourier series curve
curve = zeros(size(t));
for k = 1:numCoefficients
    curve = curve + coefficients(k) * exp(1i * 2 * pi * (k-1) * t);
end

% Plot original data and fitted curve
plot(x, y, 'ro', 'DisplayName', 'Original Data');
hold on;
plot(t, real(curve), 'b-', 'DisplayName', 'Fitted Curve');
hold off;
xlabel('x');
ylabel('y');
title('Fourier Series Curve Fitting');
legend('Location', 'best');
grid on;

% Display Fourier coefficients
disp('Fourier Coefficients:');

Fourier Coefficients:

disp(coefficients);

   0.4034 + 0.0000i  -0.1227 - 0.1916i   0.0415 + 0.0041i  -0.0248 - 0.0159i  -0.0041 - 0.0481i

Good evening;

I am trying to get fourier series curve fitting and when it runs it should give me the coefficient of the series and show me the plot.

The equation is vel=(a0/2)+(sum(a.*cos(w*n.*t')+b.*sin(w*n.*t'),2)) this but i couldnt manage to do it.

May someone please help me?

numCoefficients and time should be changable.

[Torsten]

@OLuxanna

Where did you get this formula from ?

coefficients(k) = sum(y .* exp(-1i * 2 * pi * (k-1) * x));

[OLuxanna]

Hi;

A friend of mine did try to help me to build up the code. I'm afraid I don't know.

But the one that I used before is in the previous comment.

Answer 1

[Torsten]

Two ways:

x = [0.015404238	0.027389034	0.03937383	0.051358625	0.063343421	0.075328217	0.087313012	0.099297808	0.111282604	0.1232674	0.135252195	0.147236991	0.159221787	0.171206582	0.183191378	0.195176174	0.20716097	0.219145765	0.231130561	0.243115357	0.255100152	0.267084948	0.278320694	0.28880739	0.298919562	0.308657208	0.318394855	0.329256076	0.341240872	0.353225667	0.365210463	0.377195259	0.389180054	0.40116485	0.413149646	0.425134442	0.437119237	0.449104033	0.461088829	0.473073624	0.48505842	0.497043216	0.509028012	0.521012807	0.532997603	0.544982399	0.556967194	0.56895199	0.580936786	0.592921582	0.6];  % x coordinates
y = [0.271524849	0.264952735	0.272505755	0.283825478	0.294533112	0.304063651	0.311663755	0.317992597	0.323473931	0.328484427	0.332647415	0.337846246	0.342809658	0.348526411	0.357021107	0.365421635	0.36892545	0.366355457	0.358088326	0.351139541	0.353325011	0.366245582	0.388184835	0.412362705	0.43601428	0.459267874	0.483100962	0.506751175	0.524474293	0.533957748	0.535295709	0.530277359	0.519938542	0.503149247	0.487537046	0.480682429	0.477877017	0.475872029	0.472266193	0.466965341	0.462229494	0.459235747	0.453746559	0.44288982	0.429349306	0.415385037	0.405564142	0.400828295	0.398634972	0.396629985	0.393772585];  % y coordinates
w = max(x)-min(x);
n = 5;
%Minimize sum of squared differences
A = zeros(numel(x),2*n+1);
b = zeros(numel(x),1);
A(:,1) = 1.0;
for i = 1:numel(x)
  A(i,2:n+1) = cos(2*pi/w*(1:n)*x(i));
  A(i,n+2:2*n+1) = sin(2*pi/w*(1:n)*x(i));
end
b = y.';
coeffs = A\b;
coeffs1 = coeffs.';
F1 = @(t) coeffs1(1) + sum(coeffs1(2:n+1).*cos(2*pi/w*(1:n)*t)) + sum(coeffs1(n+2:2*n+1).*sin(2*pi/w*(1:n)*t));
%Minimize L2-norm
coeffs(1) = 2/w*trapz(x,y)/2;
for i = 1:n
  coeffs(i+1) = 2/w*trapz(x,y.*cos(2*pi/w*x*i));
  coeffs(n+1+i) = 2/w*trapz(x,y.*sin(2*pi/w*x*i));
end
coeffs2 = coeffs.';
F2 = @(t) coeffs2(1) + sum(coeffs2(2:n+1).*cos(2*pi/w*(1:n)*t)) + sum(coeffs2(n+2:2*n+1).*sin(2*pi/w*(1:n)*t));

t = linspace(min(x),max(x),150);
hold on
plot(x,y,'o')
plot(t,arrayfun(@(t)F1(t),t),'r')
plot(t,arrayfun(@(t)F2(t),t),'b')
hold off
grid on

I don't have much experience in this field - so I don't know if this is the common way to compute Fourier series coefficients from data. Most probably, fft with its associated algorithms is the usual method.

Answer 2

[Hassaan]

An initial idea:

% Corrected initialization of x and y arrays with multiline syntax
x = [0.015404238, 0.027389034, 0.03937383, 0.051358625, 0.063343421, ...
    0.075328217, 0.087313012, 0.099297808, 0.111282604, 0.1232674, ...
    0.135252195, 0.147236991, 0.159221787, 0.171206582, 0.183191378, ...
    0.195176174, 0.20716097, 0.219145765, 0.231130561, 0.243115357, ...
    0.255100152, 0.267084948, 0.278320694, 0.28880739, 0.298919562, ...
    0.308657208, 0.318394855, 0.329256076, 0.341240872, 0.353225667, ...
    0.365210463, 0.377195259, 0.389180054, 0.40116485, 0.413149646, ...
    0.425134442, 0.437119237, 0.449104033, 0.461088829, 0.473073624, ...
    0.48505842, 0.497043216, 0.509028012, 0.521012807, 0.532997603, ...
    0.544982399, 0.556967194, 0.56895199, 0.580936786, 0.592921582, 0.6]; 

y = [0.271524849, 0.264952735, 0.272505755, 0.283825478, 0.294533112, ...
    0.304063651, 0.311663755, 0.317992597, 0.323473931, 0.328484427, ...
    0.332647415, 0.337846246, 0.342809658, 0.348526411, 0.357021107, ...
    0.365421635, 0.36892545, 0.366355457, 0.358088326, 0.351139541, ...
    0.353325011, 0.366245582, 0.388184835, 0.412362705, 0.43601428, ...
    0.459267874, 0.483100962, 0.506751175, 0.524474293, 0.533957748, ...
    0.535295709, 0.530277359, 0.519938542, 0.503149247, 0.487537046, ...
    0.480682429, 0.477877017, 0.475872029, 0.472266193, 0.466965341, ...
    0.462229494, 0.459235747, 0.453746559, 0.44288982, 0.429349306, ...
    0.415385037, 0.405564142, 0.400828295, 0.398634972, 0.396629985, 0.393772585]; 

% Time range
t = linspace(0, 1, 1000);

% Number of Fourier coefficients to compute
numCoefficients = 5;

% Calculate the fundamental frequency
T = max(x) - min(x); % Period
w = 2 * pi / T; % Fundamental frequency

% Initialize coefficients arrays
a = zeros(1, numCoefficients);
b = zeros(1, numCoefficients);

% Calculate the coefficients
for n = 1:numCoefficients
    a(n) = 2 * sum(y .* cos(n * w * x)) / length(x);
    b(n) = 2 * sum(y .* sin(n * w * x)) / length(x);
end

% Generate the Fourier series curve
curve = zeros(size(t));
for n = 1:numCoefficients
    curve = curve + a(n) * cos(n * w * t) + b(n) * sin(n * w * t);
end
curve = curve + a(1)/2; % Adding a0/2 term

% Plot original data and fitted curve
figure;
plot(x, y, 'ro', 'DisplayName', 'Original Data');
hold on;
plot(t, curve, 'b-', 'DisplayName', 'Fitted Curve');
hold off;
xlabel('x');
ylabel('y');
title('Fourier Series Curve Fitting');
legend('Location', 'best');
grid on;

% Display Fourier coefficients
fprintf('Fourier Coefficients:\n');

Fourier Coefficients:

fprintf('a0/2 = %f\n', a(1)/2);

a0/2 = -0.016802

for n = 1:numCoefficients-1
    fprintf('a%d = %f, b%d = %f\n', n, a(n+1), n, b(n+1));
end

a1 = 0.021770, b1 = 0.012268
a2 = 0.013380, b2 = -0.024543
a3 = 0.021739, b3 = 0.017184
a4 = 0.013186, b4 = -0.000746

-----------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

It's important to note that the advice and code are based on limited information and meant for educational purposes. Users should verify and adapt the code to their specific needs, ensuring compatibility and adherence to ethical standards.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

Feel free to contact me.

[OLuxanna]

Thank you so much Mr. Hassaan for your quick response.

I did try the code you provided but it didnt work well when i tried to fit another curve which means when i changed the x,y values.

what i would like to have as an outcome was something like that:

and these parameters should fit this equation: vel=(a0/2)+(sum(a.*cos(w*n.*t')+b.*sin(w*n.*t'),2));

I have tried to do something with the one that i got years ago like in the example below :

t= 0:0.001:0.9; %% Adjust this value based on your requirements

n=1:15; % Adjust this value based on your requirements

a=[-2235.255	50.679	64.73	56.417	93.742	68.503	64.914	52.404	57.082	50.865	45.539	40.151	34.25	30.927	24.192];

b=[792.613	-1003.748	-72.055	-267.773	-44.991	-100.476	-17.655	-35.398	-11.197	-10.323	3.272	9.12	14.151	16.328	22.273];

w=7.108;

a0=27250.392;

vel=(a0/2)+(sum(a.*cos(w*n.*t')+b.*sin(w*n.*t'),2));

[val0,idx0] = min(vel) ;
% max
[val1,idx1] = max(vel) ; 

plot(t,vel)
hold on
plot(t(idx0),vel(idx0),'*r')
plot(t(idx1),vel(idx1),'*b')
plot(t,vel,'b','LineWidth',2)

but couldnt manage to modify it for any curve

Answer 3

[Paul]

Hi OLuxanna,

The code is attempting to implement the complex exponential Fourier series using the C_n coefficients.

Let's look at the code:

% Given x and y coordinates
x = [0.015404238	0.027389034	0.03937383	0.051358625	0.063343421	0.075328217	0.087313012	0.099297808	0.111282604	0.1232674	0.135252195	0.147236991	0.159221787	0.171206582	0.183191378	0.195176174	0.20716097	0.219145765	0.231130561	0.243115357	0.255100152	0.267084948	0.278320694	0.28880739	0.298919562	0.308657208	0.318394855	0.329256076	0.341240872	0.353225667	0.365210463	0.377195259	0.389180054	0.40116485	0.413149646	0.425134442	0.437119237	0.449104033	0.461088829	0.473073624	0.48505842	0.497043216	0.509028012	0.521012807	0.532997603	0.544982399	0.556967194	0.56895199	0.580936786	0.592921582	0.6];  % x coordinates
y = [0.271524849	0.264952735	0.272505755	0.283825478	0.294533112	0.304063651	0.311663755	0.317992597	0.323473931	0.328484427	0.332647415	0.337846246	0.342809658	0.348526411	0.357021107	0.365421635	0.36892545	0.366355457	0.358088326	0.351139541	0.353325011	0.366245582	0.388184835	0.412362705	0.43601428	0.459267874	0.483100962	0.506751175	0.524474293	0.533957748	0.535295709	0.530277359	0.519938542	0.503149247	0.487537046	0.480682429	0.477877017	0.475872029	0.472266193	0.466965341	0.462229494	0.459235747	0.453746559	0.44288982	0.429349306	0.415385037	0.405564142	0.400828295	0.398634972	0.396629985	0.393772585];  % y coordinates
% Define time range
t = 0:0.001:1;
% Number of Fourier coefficients to compute
numCoefficients = 5;
% Compute Fourier coefficients
coefficients = zeros(1, numCoefficients);

These lines are trying to compute the C_n Fourier series coefficients. C_n is defined in terms of an integral, which is Eqn (7) on the linked page. The equations are trying to approximate the integral with a rectangular integration. However, if you compare to Eqn (7) you'll see that the integrand (the term inside the sum) is not correct. It's missing a term that you'll have to compute from the x data. And the coefficient also needs to be divided by that same term. And those equations aren't correct for rectangular integration, anyway. It would be easier and more accurate to use trapz instead.

for k = 1:numCoefficients
    coefficients(k) = sum(y .* exp(-1i * 2 * pi * (k-1) * x));
    coefficients(k) = coefficients(k) / length(x);
end

These equations are trying to implement Eqn (3) on the linked doc page. But the exponential is missing a term. Also, that sum in Eqn (3) goes from -N to N, and this code only goes from 0 to N (or 1 to N+1 using Matlab's 1-based indexing). But for a real function, the coefficents satisfy C(-n) = conjugate(C(n)), which you can use to correct the equation for curve.

% Generate Fourier series curve
curve = zeros(size(t));
for k = 1:numCoefficients
    curve = curve + coefficients(k) * exp(1i * 2 * pi * (k-1) * t);
end
% Plot original data and fitted curve
plot(x, y, 'ro', 'DisplayName', 'Original Data');
hold on;

If implemented correctly, the imaginary part of curve should be (very, very close to) zero, which you can use for an error check on the code.

plot(t, real(curve), 'b-', 'DisplayName', 'Fitted Curve');
hold off;
xlabel('x');
ylabel('y');
title('Fourier Series Curve Fitting');
legend('Location', 'best');
grid on;

With all of the corrections above, I got this figure.

[Torsten]

@Paul

I think one has to be a little careful with

coefficients(i) = sum(...)/length(x) 

because the x-values are not equally spaced.

[Paul]

Agreed. I did say that the equations aren't correct for rectangular integration, which is what I think was the intention. Much easier to use trapz, as suggested above. I guess I could have been more clear to use the two argument form of trapz with the vector x as the first argument. Better yet, use integral to compute the coefficients.