Я написал код для реализации алгоритма, который принимает на вход вектор q действительных чисел и возвращает на выходе сложную матрицу R. Приведенный ниже код Matlab создает график, показывающий входной вектор q и выходную матрицу R.
Учитывая только выход сложной матрицы R, я хотел бы получить входной вектор q. Могу ли я сделать это с помощью оптимизации методом наименьших квадратов? Поскольку в коде есть рекурсивная текущая сумма (rs_r и rs_i), вычисление столбца выходной матрицы зависит от вычисления предыдущего столбца.
Возможно, можно настроить нелинейную оптимизацию для перекомпоновки входного вектора q из выходной матрицы R?
Глядя на это с другой стороны, я использовал алгоритм для вычисления матрицы R. Я хочу запустить алгоритм «в обратном порядке», чтобы получить входной вектор q из выходной матрицы R.
Если нет возможности перекомпоновать начальные значения из вывода, тем самым рассматривая проблему как «черный ящик», тогда, возможно, при оптимизации можно использовать математику самой модели? Программа оценивает следующее уравнение:
Утильда (тау, омега) — это выходная матрица R. Переменная тау (время) содержит столбцы матрицы отклика R, тогда как переменная омега (частота) содержит строки матрицы отклика R. Интеграция выполняется как рекурсивная текущая сумма от tau = 0 до текущего временного шага tau.
Вот графики, созданные программой, размещенной ниже:
Вот полный код программы:
N = 1001;
q = zeros(N, 1); % here is the input
q(1:200) = 55;
q(201:300) = 120;
q(301:400) = 70;
q(401:600) = 40;
q(601:800) = 100;
q(801:1001) = 70;
dt = 0.0042;
fs = 1 / dt;
wSize = 101;
Glim = 20;
ginv = 0;
R = get_response(N, q, dt, wSize, Glim, ginv); % R is output matrix
rows = wSize;
cols = N;
figure; plot(q); title('q value input as vector');
ylim([0 200]); xlim([0 1001])
figure; imagesc(abs(R)); title('Matrix output of algorithm')
colorbar
Вот функция, которая выполняет расчет:
function response = get_response(N, Q, dt, wSize, Glim, ginv)
fs = 1 / dt;
Npad = wSize - 1;
N1 = wSize + Npad;
N2 = floor(N1 / 2 + 1);
f = (fs/2)*linspace(0,1,N2);
omega = 2 * pi .* f';
omegah = 2 * pi * f(end);
sigma2 = exp(-(0.23*Glim + 1.63));
sign = 1;
if(ginv == 1)
sign = -1;
end
ratio = omega ./ omegah;
rs_r = zeros(N2, 1);
rs_i = zeros(N2, 1);
termr = zeros(N2, 1);
termi = zeros(N2, 1);
termr_sub1 = zeros(N2, 1);
termi_sub1 = zeros(N2, 1);
response = zeros(N2, N);
% cycle over cols of matrix
for ti = 1:N
term0 = omega ./ (2 .* Q(ti));
gamma = 1 / (pi * Q(ti));
% calculate for the real part
if(ti == 1)
Lambda = ones(N2, 1);
termr_sub1(1) = 0;
termr_sub1(2:end) = term0(2:end) .* (ratio(2:end).^-gamma);
else
termr(1) = 0;
termr(2:end) = term0(2:end) .* (ratio(2:end).^-gamma);
rs_r = rs_r - dt.*(termr + termr_sub1);
termr_sub1 = termr;
Beta = exp( -1 .* -0.5 .* rs_r );
Lambda = (Beta + sigma2) ./ (Beta.^2 + sigma2); % vector
end
% calculate for the complex part
if(ginv == 1)
termi(1) = 0;
termi(2:end) = (ratio(2:end).^(sign .* gamma) - 1) .* omega(2:end);
else
termi = (ratio.^(sign .* gamma) - 1) .* omega;
end
rs_i = rs_i - dt.*(termi + termi_sub1);
termi_sub1 = termi;
integrand = exp( 1i .* -0.5 .* rs_i );
if(ginv == 1)
response(:,ti) = Lambda .* integrand;
else
response(:,ti) = (1 ./ Lambda) .* integrand;
end
end % ti loop
