我试图重复我在一篇论文中找到的一个例子。
我必须解决这个 ODE:
25 a + 15 v + 330000 x = p(t)
其中p(t)是频带限制在 10-25 Hz 范围内的白噪声序列;a是加速度,v是速度,x是位移。
然后它说系统是使用 Runge-Kutta 程序模拟的。采样频率设置为 1000 Hz,并将高斯白噪声添加到数据中,使噪声占信号 rms 值的 5%(如何使用最后的信息?)。
主要问题与带限白噪声有关。我按照我在这里找到的说明https://dsp.stackexchange.com/questions/9654/how-to-generate-band-limited-gaussian-white-noise并编写了以下代码:
% 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);
现在我必须为 ode 定义函数,但我不知道如何给它p(白噪声)......我试过这样:
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);
end
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,options);
事实是输入激励不能正常工作:方程的自然频率约为 14 Hz,因此我希望看到响应中的共振,因为白噪声在 10-25 Hz 范围内。
我也看过这个 Q/A,但我仍然无法让它工作: