ОПТИМАЛЬНАЯ ИНТЕРПОЛЯЦИЯ
Просеков О. В. (sc2@pisem.net)
2 ноября 2004 г.
Санкт-Петербургский городской семинар «Всплески (wavelets) и их приложения». Секция «Дискретный гармонический анализ»: http://www.math.spbu.ru/user/dmp/dha/
Contents
Постановка задачи
Пусть N = mn, где n > 1, и r - натуральное число. Рассматривается экстримальная задача:
В ней требуется построить возможно более гладкий сигнал принимающий в узлах ln заданные значания z(l). Гладкость характеризуется квадратом нормы конечной разности r-го порядка. Обычно r = 2.
Исходные данные
r = 2;
n = 8;
m = 8;
N = m*n;
rand('state', 0);
z = rand(m, 1);
Инициализация локальных данных
зависимость от 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
Вычисление промежуточного массива констант
зависимость от r:
A = 4 * sin(pi/N * A) .^ (-2*r); % n x m-1 A = ifft(A).'; % n x m-1
Вычисление массивов констант
lambda = 1 ./ A(:,1); % m-1 x 1 D = A(:,2:end) .* omega; % m-1 x n-1
Оптимальная интерполяция
frc = isreal(z); x = zeros(m, n); B = zeros(m, n-1);
x(:,1) = z; Z = fft(z);
Z(2:end) = lambda .* Z(2:end);
B(1,:) = Z(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(:,2:end) = B; x = reshape(x.', [], 1); if frc, x = real(x); end
Для сравнения воспользуемся втроенной в MATLAB функцией interpft
y = interpft(z, N);
Графики сигналов
j = 0:N-1; l = 0:m-1; figure hold on stem(j, y, 'm', '-o', 'MarkerSize', 5) stem(j, x, 'g', '-o', 'MarkerSize', 5) stem(l*n, z, 'b', '-o') legend('interpft', 'optinterp', 'original') hold off
ЛИТЕРАТУРА
Малозёмов В. Н., Машарский С. М. Основы дискретного гармонического анализа. Часть 1. СПб.: НИИММ, 2003. С. 51-57.