The problem with using a two-element vector for ‘tspan’ is that the MATLAB ODE integrators do not produce regularly-spaced time values with it, and the fft function requires that they must be regularly-spaced (constant sampling intervals and constant sampling frequency). To get them to be regularly-spaced, it is necessary to define them as such.
Try this —
% Parameters
c1 = 0.3;
k1 = 3;
k2 = 0.5;
w = 1.7; % rad/s
A = 0.5;
% Initial conditions
y0 = [0; 0]; % y1(0) = 0, y2(0) = 0
%Time span
% tspan = [0, 50];
Fs = 1000; % Sampling Frequency (Hz)
L = 60; % Signal Length (s)
tspan = linspace(0, L*Fs, L*Fs+1)/Fs % Time Vector
% ODE Solver
[t, y] = ode45(@duffing, tspan, y0, [], c1, k1, k2, w, A);
% Y1 = fft(y(:,1)); % y1 signal Fourier transform
[Y1,Fv] = FFT1(y(:,1),t);
[Y1max,idx] = max(mag2db(abs(Y1)));
Maximum_Magnitude_dB = Y1max
Peak_Frequency_Hz = Fv(idx)
figure
subplot(2,1,1)
semilogx(Fv, mag2db(abs(Y1)))
ylabel('Magnitude (dB)')
grid
axis('padded')
subplot(2,1,2)
semilogx(Fv, rad2deg(unwrap(angle(Y1))))
grid
xlabel('Frequency (Hz)')
ylabel('Phase (°)')
axis('padded')
% Duffing equation function
function dydt = duffing(t,y,c1,k1,k2,w,A)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -c1*y(2) - k1*y(1) - k2*y(1)^3 + A*sin(w.*t);
end
function [FTs1,Fv] = FFT1(s,t)
s = s(:);
t = t(:);
L = numel(t);
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv);
end
.
