論文で見つけた例を再現しようとしています。
この ODE を解決する必要があります。
25 a + 15 v + 330000 x = p(t)
ここで、p(t)は 10 ～ 25 Hz の範囲に帯域制限されたホワイト ノイズ シーケンスです。aは加速度、vは速度、xは変位です。
次に、システムはルンゲクッタ手順を使用してシミュレートされていると表示されます。サンプリング周波数は 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);
ここで、頌歌の関数を定義する必要がありますが、 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 も見ましたが、まだ機能させることはできません。