I am trying to repeat an example I found in a paper.

I have to solve this ODE:

$25 \ddot{x} + 15 \dot{x} + 330000 x = p(t)$

where $p(t)$ is a white noise sequence band-limited into the 10-25 Hz range.

Then it says the system is simulated using a Runge-Kutta procedure. The sampling frequency is set to 1000 Hz and a Gaussian white noise is added to the data in such a way that the noise contributes to 5% of the signal r.m.s value (how do I use this last info?).

The main problem is related to the band-limited white noise. I followed the instructions I found here https://dsp.stackexchange.com/questions/9654/how-to-generate-band-limited-gaussian-white-noise and I wrote the following code:

% White noise excitation
% design FIR filter to filter noise in the range [10-25] Hz
Fs = 1000; % sampling frequency

% Time infos (for white noise)
T = 1/Fs;
tspan = 0:T:4; % I saw from the plots in the paper that the simulation last 4 seconds
tstart = tspan(1);
tend = tspan (end);

b = fir1(64, [10/(Fs/2) 25/(Fs/2)]);
% generate Gaussian (normally-distributed) white noise
noise = randn(length(tspan), 1);
% apply filter to yield bandlimited noise
p = filter(b,1,noise);

Now I have to define the function for the ode, but I do not know how to give it the $p$ (white noise):

I tried coding this way

odefun = @(t,u,m,k1,c,p)[u(2); 1/m*(-k1*u(1)-c*u(2)+p)];
q0 = [0;0]; % Initial conditions
[t,q,te,qe,ie] = ode45(@(t,u)odefun(t,u,m,k1,c,p),tspan,q0(:,1),options);

But I get this error message:

@(T,U)ODEFUN(T,U,M,K1,C,NB) returns a vector of length 4002, but the length of initial conditions vector is 2. The vector returned by @(T,U)ODEFUN(T,U,M,K1,C,NB) and the initial conditions vector must have the same number of elements.

I already had a look at the following questions but I do not understand where the problem could be.





I might have found a way to do it:

function [y] = p(t)
    b = fir1(64, [10/(Fs/2) 25/(Fs/2)]);
    % generate Gaussian (normally-distributed) white noise
    noise = randn(length(t), 1);
    % apply filter to yield bandlimited noise
    y = filter(b,1,noise);

odefun = @(t,u,m,k1,c,p)[u(2); 1/m*(-k1*u(1)-c*u(2)+p(t))];

[t,q,te,qe,ie] = ode45(@(t,u)odefun2(t,u,m,k2,c,@p),tspan,q0(:,1),options);

Does it look fine to you?

  • $\begingroup$ RK-integration looks like overkill to me. Your system is a linear filter with transfer function G(s)=1/(25*s^2+15*s+330000). So, one could just apply the filter as you do with the band limitation. But, maybe there is a special purpose of the simulation which I do not know. $\endgroup$
    – Tobias
    Nov 23, 2014 at 10:49
  • $\begingroup$ I am making this simulation to study the trend of the system restoring force as function of the displacement (the result should be a linear trend) $\endgroup$
    – Rhei
    Nov 23, 2014 at 11:05


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.