FFT fitting of a damped sinusoid

Hello, I have the following code. For some reason I can't figure out why I can't fit my data with an FFT that has the correct period, phase, and amplitude. I can't use standard fourier series because I want to obtain the spectral graphs. I can only fit a single period without problems.

tspan=[0:0.1:30];
x0=[0.2;0.5];


[tf,xf]=ode45(@func,tspan,x0);

xout=xf(:,1);
xdot=xf(:,2);


figure
[x1,x2]=meshgrid(4:0.09:9,-3:0.09:6);
k=8;
c=2;
m=5;
U=x2;
V=-(c/m)*x2-(k/m)*x1+10;
L=sqrt(U.^2+V.^2);
quiver(x1,x2,U./L,V./L,0.7,'b');
xlabel('x(m)')
ylabel('xdot (m/s)')
title('Phase Portrait')
hold on

plot(xout,xdot,'r','linewidth',2)

xlim([4.8 8]) 
ylim([-2 1.5])

hold off

dt=tspan(end)/length(tspan);

T=5;
NSP= floor(T/dt);

yfit = xout(1:6*NSP);

coefs=fft((1/NSP)*yfit);

tr=0:0.01:T;

An=-imag(coefs); 

Bn=real(coefs);

ynew=0;

for N=1:NSP


ynew=ynew+An(N)*sin((N-1)*2*pi/T*tr);
ynew=ynew+Bn(N)*cos((N-1)*2*pi/T*tr);



end




figure
plot(tr,ynew,"r",'linewidth',2)
title('Mass 2 Response')
xlabel('time(s)')
ylabel('amplitude(m)')

hold on

plot(tspan,xout,"bo")

hold off


% figure
% plot(abs(coefs))
% xlabel('Frequency (Hz)')
% ylabel('Magnitude')
% title('Mass 2 Magnitude Plot')
% 

function xdot=func(t,x)
k=8;
c=2;
m=5;

xdot(1)=x(2);
xdot(2)=-(c/m)*x(2)-(k/m)*x(1)+10;

xdot=xdot';
end

Answer 1

[Sudarsanan A K]

Hi Ibrahim,

It is my understanding that you are trying to fit your data using the Fast Fourier Transform (FFT) to obtain the spectral graphs and facing issue with plotting the reconstructed signal based on the Fourier coefficients ('An' and 'Bn').

In this regard, I have made some slight modifications to your code as follows so as to address the desired use-case:

• Modified 'T' Calculation: I modified the calculation of 'T' to use the full time span 'tspan(end)'. This ensures that 'T' represents the correct period of the signal.
• Removed Subset Selection for 'yfit': In your code, there is a line 'yfit = xout(1:6*NSP);' that selected a subset of the data for fitting. I removed this line to ensure that the fitted signal is generated using the entire 'xout' data.
• Revised 'tr' Time Vector: I modified the variable 'tr' by using 'linspace' to create a linearly spaced vector that covers the full time span. This ensures that the fitted signal is plotted correctly over the entire time range.

Please find the modified code attached below.

tspan = 0:0.1:30;
x0 = [0.2; 0.5];

[tf, xf] = ode45(@func, tspan, x0);

xout = xf(:, 1);
xdot = xf(:, 2);

figure
[x1, x2] = meshgrid(4:0.09:9, -3:0.09:6);
k = 8;
c = 2;
m = 5;
U = x2;
V = -(c/m) * x2 - (k/m) * x1 + 10;
L = sqrt(U.^2 + V.^2);
quiver(x1, x2, U./L, V./L, 0.7, 'b');
xlabel('x(m)')
ylabel('xdot (m/s)')
title('Phase Portrait')
hold on
plot(xout, xdot, 'r', 'linewidth', 2)
xlim([4.8 8])
ylim([-2 1.5])
hold off

dt = tspan(end) / length(tspan);
T = tspan(end);    % modified T to use the full time span tspan(end)
NSP = floor(T / dt);
coefs = fft((1 / NSP) * xout);
tr = linspace(0, T, length(tspan)); % used 'linspace' to create a linearly spaced vector that covers the full time span
An = -imag(coefs);
Bn = real(coefs);
ynew = zeros(size(tr));
for N = 1:NSP
    ynew = ynew + An(N) * sin((N-1) * 2*pi/T * tr) + Bn(N) * cos((N-1) * 2*pi/T * tr);
end

figure
plot(tr, ynew, 'r', 'linewidth', 2)
title('Mass 2 Response')
xlabel('time(s)')
ylabel('amplitude(m)')
hold on
plot(tspan, xout, 'bo')
grid on
hold off

function xdot = func(~, x)
    k = 8;
    c = 2;
    m = 5;

    xdot(1) = x(2);
    xdot(2) = -(c/m) * x(2) - (k/m) * x(1) + 10;

    xdot = xdot';
end

I hope this helps you to resolve your issue.