Я пытаюсь повторить пример, который нашел в статье.
Я должен решить эту ОДУ:
25 a + 15 v + 330000 x = p(t)
где p(t) — последовательность белого шума, ограниченная полосой частот в диапазоне 10–25 Гц; a — ускорение, v — скорость и x — перемещение.
Затем он говорит, что система моделируется с использованием процедуры Рунге-Кутты. Частота дискретизации установлена на 1000 Гц, и к данным добавляется гауссов белый шум таким образом, что шум составляет 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 Гц, поэтому я ожидаю увидеть резонанс в отклике, поскольку белый шум находится в диапазоне 10-25 Гц.
Я также просмотрел этот Q / A, но я все еще не могу заставить его работать:
Как решить систему обыкновенного дифференциального уравнения с параметрами, зависящими от времени
Решение ОДУ при функция задается дискретными значениями -matlab-
awgn
. Или см. этот вопрос. - person horchler   schedule 13.12.2014ode45
, попробовать жесткий решатель, такой какode15s
, либо реализовать фиксированный шаг- размер схемы Рунге-Кутта самостоятельно и использовать достаточно маленькие шаги (это никогда не хорошо для репликации, когда они даже не удосуживаются сообщить конкретный используемый решатель и размер шага, если это был фиксированный метод). - person horchler   schedule 13.12.2014