How to obtain FFT of nonuniformly spaced data?

6 vues (au cours des 30 derniers jours)
lekhani gaur
lekhani gaur le 22 Juil 2021
Modifié(e) : Paul le 3 Sep 2021
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 commentaires
dpb
dpb le 22 Juil 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.
lekhani gaur
lekhani gaur le 23 Juil 2021
Modifié(e) : lekhani gaur le 2 Août 2021
Thanks a lot. It worked.

Connectez-vous pour commenter.

Réponse acceptée

Paul
Paul le 24 Juil 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 commentaires
lekhani gaur
lekhani gaur le 2 Sep 2021
@Paul Thank you so much.
1. This method worked for certain range of the parameters, specially for those where i was facing issues of inconsistent initial conditions. However, ODE45 could not solve for some values of parameters involved in my equations, where ODE15i is giving results.
2. I am keeping y values (y1 and y3) fixed for my next set of equation allowing first and second derivative (yp(1) and yp(2)) to change by using 1s and 0s in decic function. Also, It is needed to check whether the consistent initial condtions obtained from decic functions are realistic. It does not always give realistic values while fixing certain yp values. it is required to check for certain combinations of 1s and 0s in y and yp.
3. Yes, the Modified Equations are an attempt to force ode15i to return the values of yp(2) and yp(4) in y(5) and y(6). It did not yield any solution (all y are zero), as you may check with the code.
Paul
Paul le 3 Sep 2021
Modifié(e) : Paul le 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,[])
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

Connectez-vous pour commenter.

Plus de réponses (2)

Mateo González Hurtado
Mateo González Hurtado le 25 Août 2021
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 le 2 Sep 2021
doc ode15i
  1 commentaire
Walter Roberson
Walter Roberson le 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.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by