How to obtain FFT of nonuniformly spaced data?
5 views (last 30 days)
Show older comments
Greetings !
I have obtained displacement time history of a sytem using ODE15i. I need to obtain FFT curve of this response. As ODE solvers use adaptive algorithm, the time interval between steps is not constant. I have tried both hanning windowing option and zero padding, but not successful. It is giving broad spectrum instead of peaks at certain frequencies. I tried NUFFT command too, but it doesn't seem to help it out. If I fix the time step while solving ODE, then only I get good FFT curve. I guess, It is not good option to fix the time step in adaptive algorithms. Here, I am attaching the code and the data for your reference which contain both fixed step and variable step data.
Kindly suggest what should be done in this scenario. Thank you.
2 Comments
dpb
on 22 Jul 2021
For an analytic solution from an ODE solver the solution is smooth between sampled output points so just resample the time signal to a uniform sample rate and use the normal FFT.
Only if you had a discrete time sampled signal would the above not be true, but even there if there's something going on between sampled points faster than the slowest sample rate/longest sample interval, you're missing it entirely, anyway, so there's nothing to be learned/gained nor information lost by resampling it as well.
Accepted Answer
Paul
on 24 Jul 2021
You can try using the tspan input to od15i to specify equally spaced sample points in time at which the solution is desired. Check
doc ode15i
for more information.
26 Comments
Paul
on 3 Sep 2021
Edited: Paul
on 3 Sep 2021
1. I only used ode45 to illustrate that an explicit solver is feasible. Another solver (along with appropriate options) might be more appropriate for your problem. This link (Choose an ODE Solver) may be of interest if you haven't read it already.
2. Without seeing an actual example, it's hard to say anything more about resetting the initial conditions at the transitions, other than to say it sounds like you're holding y(1) and y(3) fixed, and allowing y(2), y(4) and yp(1:4) to float.
3. Here is my attempt to have ode15i return the second derivatives.
global t0 A1 A2 a B1 B2 C1 C2 D1 D2
% define parameters
A1 = 300;
A2 = 300;
a = 0.3*A1;
B1 = 1E6;
B2 = 1E6;
C1 = 1E2;
C2 = 1E2;
D1 = 1E3;
D2 = 1E3;
t0 = 0.1;
% initial condition for BASIC EQUATIONS
% y(1) and y(3) must be zero
y0=[0;0;0;0];
yp0=[0;0;0;0];
tspan_a = [0 1];
Even though the stated desire is to have y(1) and y(3) initialized to zero, the call to decic allows them to vary. The result of the call to decic is that y(1) ~=0. Maybe use the fixed_y0 input to decic to hold y(1) and y(3) at zero.
[y0_new,yp0_new,resrnm1] = decic(@diff,0,y0,[],yp0,[])
The call to ode15i should use y0_new as the initial condition as that was the output from decic().
%[t,y]=ode15i(@diff,tspan_a,y0,yp0_new,options); % this should be y0_new!
odefunz = @(t,z,zp) ([diff(t,z(1:4),zp(1:4)) ; z(5:6)-zp([2 4])]);
options = odeset('RelTol',1e-6,'AbsTol',1e-8,"Jacobian",@Fjac2,"Events",@dydt1Events);
Might want to consider setting MaxSep and InitialStep.
z0 = [y0_new; yp0_new(2); yp0_new(4)];
zp0 = [yp0_new; 0; 0];
[t,z] =ode15i(odefunz,tspan_a,z0,zp0,options);
plot(t,z(:,1:4)),grid
This result looks identical(ish) to the result in the comment above using ode45
Now check that z(5:6) returned from ode15i are the second derivatives of y(1) and y(3), or the first derivatives of y(2) and y(4), i..e,. yp(2) and yp(4)
subplot(211);plot(t,z(:,5),t,gradient(z(:,2),t),'o'),grid
subplot(212);plot(t,z(:,6),t,gradient(z(:,4),t),'o'),grid
Looks very reasonable.
function dydt = diff(t,y,yp)
global t0 A1 A2 a B1 B2 C1 C2 D1 D2
% define load
if t<=t0
F = 1E6*(1-t/t0);
else
F=0;
end
% second order coupled differential equations (BASIC EQUATION)
% yp(1) is 1st derivative of y1 and written as y(2) ;
% yp(2) is 2nd derivative of y1 ;
% yp(3) is 1st derivative of y3 and written as y(4) ;
% yp(4) is 2nd derivative of y3 ;
%% A1*d2/dt2(y1)+a*(d2/dt2(y1)-d2/dt2(y3))+B1*y1+D1*(y1-y3)+C1*(d/dt(y1)-d/dt(y3))=F
%% A2*d2/dt2(y3)+a*(d2/dt2(y3)-d2/dt2(y1))+B2*y3+D2*(y3-y1)+C1*(d/dt(y3)-d/dt(y1))=0
dydt=[yp(1)-y(2);
yp(3)-y(4);
A1*yp(2)+a*(yp(2)-yp(4))+B1*y(1)+D1*(y(1)-y(3))+C1*(y(2)-y(4))-F;
A2*yp(4)+a*(yp(4)-yp(2))+B2*y(3)+D2*(y(3)-y(1))+C2*(y(4)-y(2))];
end
% jacobian for (BASIC EQUATION)
function [dfdy, dfdp] = Fjac1(t,y,yp)
global A1 A2 a B1 B2 C1 C2 D1 D2
dfdy=[0 -1 0 0;
0 0 0 -1;
B1+D1 C1 -D1 -C1;
-D2 -C2 B2+D2 C2];
dfdp=[1 0 0 0;
0 0 1 0;
0 A1+a 0 -a;
0 -a 0 A2+a];
end
function [dfdz, dfdzp] = Fjac2(t,z,zp)
[dfdy, dfdyp] = Fjac1(t,z(1:4),zp(1:4));
dfdz = blkdiag(dfdy,eye(2));
dfdzp = [dfdyp zeros(4,2); -[0 1 0 0 0 0; 0 0 0 1 0 0] ];
end
function [value,isterminal,direction] = dydt1Events(t,y,yp)
value = y(2);
isterminal = 1;
direction = -1;
end
More Answers (2)
Mateo González Hurtado
on 25 Aug 2021
Hi,
If it is not your intention to recalculate the response with ODE: have you tried to interpolate the non-uniform data?
Regards,
0 Comments
S Priyadharshini
on 2 Sep 2021
doc ode15i
1 Comment
Walter Roberson
on 2 Sep 2021
Please be more specific about what the poster should be looking for in that documentation? Paul already referred to that documentation more than a month ago -- and followed up with reference to particular points of information to pay attention to.
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!