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.
