MATLAB Answers

How to obtain FFT of nonuniformly spaced data?

조회 수: 141(최근 30일)
lekhani gaur
lekhani gaur 2021년 7월 22일
편집: Paul 2021년 9월 3일 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
lekhani gaur
lekhani gaur 2021년 7월 23일
Thanks a lot. It worked.

댓글을 달려면 로그인하십시오.

채택된 답변

Paul
Paul 2021년 7월 24일
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
Paul
Paul 2021년 9월 3일 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

댓글을 달려면 로그인하십시오.

추가 답변(2개)

Mateo González Hurtado
Mateo González Hurtado 2021년 8월 25일
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 2021년 9월 2일 9:24
doc ode15i
  댓글 수: 1
Walter Roberson
Walter Roberson 2021년 9월 2일 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.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by