function [x, f] = optinterp(z, nn, rr)
% ОПТИМАЛЬНАЯ ИНТЕРПОЛЯЦИЯ
% Просеков О. В. (sc2@pisem.net)
% 2 ноября 2004 г.
%
% Санкт-Петербургский городской семинар «Всплески (wavelets) и их приложения».
% Секция «Дискретный гармонический анализ»: http://www.math.spbu.ru/user/dmp/dha/
% 
%% ОПИСАНИЕ %
% OPTINTERP (optimal interpolation) 
%
%   [x, f] = optinterp(z, n)
%   [x, f] = optinterp(z, n, r)
%   [x, f] = optinterp(z)
%
% z - интерполируемый сигнал    %  m x 1
% n - число точек интерполяции  %  1 x 1
% r - степень гладкости (r>1)   %  1 x 1
%
% x - результат интерполяции    % mn x 1
% f - значение целевой функции  %  1 x 1
%
%% ЗАДАЧА %
% Пусть $N = m\,n$, где $n \ge 2$, и $r$~--- натуральное число.
% Рассматривается экстримальная задача:
% \begin{gather*}
% f(x) :=  \|\Delta^r(x)\|^2 \to \min\,, \\
% x(l\,n) = z(l)\,, \quad l \in 0:m-1\,; \quad x \in \mathbb{C}_N\,.
% В ней требуется построить возможно более гладкий сигнал принимающий в
% узлах $l\,n$ заданные значания $z(l)$. Гладкость характеризуется
% квадратом нормы конечной разности $r$-го порядка. По умолчанию $r = 2$.
%
%% ПРИМЕР %
% r = 2;
% n = 8;
% m = 8;
% 
% rand('state', 0);
% z = rand(m, 1);
% 
% [x, f] = optinterp(z, n, r);
%
%% ЛИТЕРАТУРА %
% Малозёмов В. Н., Машарский С. М. Основы дискретного гармонического анализа.
% Часть 1. СПб.: НИИММ, 2003. С. 51-57.

%% Локальные данные
persistent m n r lambda D

% Инициализация параметров алгоритма
if (nargin == 2)
    rr = 2;
end

if nargin ~= 1
    m = length(z);
    n = nn;
    r = rr;
    [lambda, D] = initialization(m, n, r);
end

% Оптимальная интерполяция
[x, z] = processing(z, lambda, D, m, n);

if nargout == 2
    % Вычисление значения целевой функции
    % $$ f(x_*) = \frac 1 m \sum_{p=1}^{m-1} \lambda_p \, \|Z(p)\|^2 $$
    f = 1/m * sum(lambda .* abs(z(2:end)).^2);
end


%% Локальные функции


function [lambda, D] = initialization(m, n, r)
% Инициализация локальных данных
    % --- зависимость от m и n ---
N = m*n;
q = (0:n-1)';
p = 1:m-1;
s = 1:n-1;
omega_N = exp(pi/N*2i);
omega = omega_N .^ (p'*s);                 % m-1 x n-1
A = repmat(p, n, 1) + repmat(q*m, 1, m-1); % n x m-1
% Вычисление промежуточного массива констант
% $$ A[p,q] = 4 \sin\Big(\frac \pi N (p+qm)\Big), \qquad p = 1:m-1, \quad q = 0:n-1 $$
    % --- зависимость от r ---
A = 4 * sin(pi/N * A) .^ (-2*r);           % n x m-1
A = ifft(A).';                             % n x m-1
% Вычисление массивов констант
% $$ \lambda_p = 1/A[p,1]\,, \quad D[s,p] = A[p,s] \, \omega_N^{p\,s}\,, \qquad p \in 1:m-1\,, \quad s = 1:m-1\,. $$
lambda = 1 ./ A(:,1);                      % m-1 x 1
D = A(:,2:end) .* omega;                   % m-1 x n-1


function [x, Z] = processing(Z, lambda, D, m, n)
% Оптимальная интерполяция
frc = isreal(Z);
x  = zeros(m, n);
B  = zeros(m, n-1);
% $$ x(ln) = z(l), \qquad l = 0:m-1 $$
x(:,1) = Z;
Z = fft(Z);
% $$ \tilde{Z}(p) = \lambda_p Z(p), \qquad p = 1:m-1 $$
Z(2:end) = lambda .* Z(2:end);
% $$ B[s,0] = Z(0), \qquad s = 1:n-1 $$
B(1,:) = Z(1);
% $$ B[s,p] = \tilde{Z}(p) \, D[s,p], \qquad s = 1:n-1, \quad p = 1:m-1 $$
B(2:end,:) = repmat(Z(2:end), 1, n-1) .* D; % m-1 x n-1
B = ifft(B);                                % m x n-1
% $$ x(s+ln) = B[s,l], \qquad l = 0:m-1, \quad s = 1:n-1 $$
x(:,2:end) = B;
x = reshape(x.', [], 1);
if frc, x = real(x); end