Clear Filters
Clear Filters

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

8 views (last 30 days)
% 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.
  2 Comments
OLuxanna
OLuxanna on 28 Mar 2024
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.

Sign in to comment.

Answers (3)

Torsten
Torsten on 28 Mar 2024
Edited: Torsten on 28 Mar 2024
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.

Hassaan
Hassaan on 28 Mar 2024
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.
  1 Comment
OLuxanna
OLuxanna on 28 Mar 2024
Moved: Torsten on 28 Mar 2024
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

Sign in to comment.


Paul
Paul on 29 Mar 2024
Edited: Paul on 29 Mar 2024
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.
  2 Comments
Torsten
Torsten on 29 Mar 2024
Edited: Torsten on 29 Mar 2024
I think one has to be a little careful with
coefficients(i) = sum(...)/length(x)
because the x-values are not equally spaced.
Paul
Paul on 29 Mar 2024
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.

Sign in to comment.

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!