MATLAB Answers

How to obtain FFT of nonuniformly spaced data?

159 views (last 30 days)
lekhani gaur
lekhani gaur on 22 Jul 2021
Edited: Paul on 3 Sep 2021 at 0:32
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

Sign in to comment.

Accepted Answer

Paul
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
Paul on 3 Sep 2021 at 0:30
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,[])
y0_new = 4×1
0.0048 0 0 0
yp0_new = 4×1
1.0e+03 * 0.0000 2.6954 -0.0000 0.6220
resrnm1 = 4.6932e-10
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

Sign in to comment.

More Answers (2)

Mateo González Hurtado
Mateo González Hurtado on 25 Aug 2021 at 9:13
Hi,
If it is not your intention to recalculate the response with ODE: have you tried to interpolate the non-uniform data?
The answer posted on this link might be helpful.
Regards,

S Priyadharshini
S Priyadharshini on 2 Sep 2021 at 9:24
doc ode15i
  1 Comment
Walter Roberson
Walter Roberson on 2 Sep 2021 at 9:41
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by