Сжатое зондирование и восстановление сигнала (OMP)

В этом примере мы рассматриваем классическую задачу обработки сигналов: имеется сигнал большой длины n = 256, но известно, что он разрежен — только k = 20 его отсчётов отличны от нуля, остальные равны нулю. Вместо того чтобы измерять все 256 значений, мы делаем лишь m = 128 линейных измерений вида y = A·x, где A — случайная гауссова матрица размера m × n. Задача — по «сжатым» измерениям y восстановить исходный сигнал x.

В чём суть технологии. Классическая теория (теорема Котельникова) требует числа измерений порядка длины сигнала. Однако если сигнал разрежен, то теория сжатого зондирования (Compressed Sensing) гарантирует точное восстановление уже при m ≈ k·log(n), то есть при m ≪ n. Ключевая идея — использовать нелинейный разреженный решатель вместо обычного МНК.

Алгоритм OMP (Orthogonal Matching Pursuit, ортогональное преследование согласованием) — один из самых простых жадных алгоритмов восстановления. На каждой итерации он:

находит столбец матрицы A, наиболее коррелирующий с текущим остатком r = y − A·x;
добавляет его индекс в поддерживаемое множество;
решает задачу наименьших квадратов только на выбранных столбцах;
обновляет остаток.

Через k итераций алгоритм возвращает разреженное решение с k ненулевыми компонентами. В нашем примере он восстанавливает сигнал с относительной ошибкой порядка 10⁻¹⁵ — то есть практически идеально, что наглядно демонстрируется графиками исходного и восстановленного сигналов.

Реализация

В примере используются три библиотеки:

LinearAlgebra — стандартная библиотека линейной алгебры: матричные операции, нормы, решение линейных систем (используется для norm, A', Asub \ y).
SparseArrays — стандартная библиотека для работы с разреженными векторами и матрицами (используется для sparsevec, findnz, эффективного хранения сигнала с малым числом ненулевых элементов).
Random — стандартная библиотека генерации случайных чисел (используется для randn, randperm и Random.seed! для воспроизводимости).

Pkg.add(["LinearAlgebra", "SparseArrays", "Random"])
using LinearAlgebra, SparseArrays, Random

omp(A, y, k) — реализация алгоритма Orthogonal Matching Pursuit (OMP) для жадного восстановления разреженного вектора x из уравнения y = A·x, где k — ожидаемое число ненулевых элементов.

Алгоритм:

Инициализирует: r = y, support = [], x = zeros(n).
Повторяет k раз:
- вычисляет корреляции A' * r и берёт их модуль;
- обнуляет уже выбранные индексы;
- выбирает индекс j с максимальной корреляцией и добавляет его в support;
- строит подматрицу Asub = A[:, support];
- решает МНК-задачу Asub \ y, получает коэффициенты coef;
- обновляет x: обнуляет и записывает coef в позиции support;
- пересчитывает остаток r = y - Asub * coef;
- если norm(r) < 1e-12, прерывает цикл.
Возвращает x.

function omp(A, y, k)
    m, n = size(A)
    r = copy(y)
    support = Int[]
    x = zeros(n)
    for _ in 1:k
        c = abs.(A' * r)
        c[support] .= 0
        j = argmax(c)
        push!(support, j)
        Asub = A[:, support]
        coef = Asub \ y
        x .= 0
        x[support] = coef
        r = y - Asub * coef
        if norm(r) < 1e-12
            break
        end
    end
    return x
end

omp (generic function with 1 method)

Подготовка данных для сжатого зондирования:

Задаём параметры: n = 256, m = 128, k = 20, где m < n обеспечивает сжатие.
Фиксируем генератор: Random.seed!(42) для воспроизводимости.
Создаём разреженный сигнал: x_true = zeros(n), выбираем k случайных индексов через randperm(n)[1:k] и записываем туда randn(k).
Строим матрицу измерений A = randn(m, n) и нормируем её столбцы: A = A ./ sqrt.(sum(A.^2, dims=1)) — чтобы OMP не давал преимущество столбцам с большей энергией.
Вычисляем сжатые измерения: y = A * x_true.

n = 256                # длина сигнала
m = 128                # число измерений (m < n)
k = 20                 # разреженность (число ненулевых элементов)

Random.seed!(42)
x_true = zeros(n)
support = randperm(n)[1:k]
x_true[support] = randn(k)
A = randn(m, n)
A = A ./ sqrt.(sum(A.^2, dims=1))
y = A * x_true;

Вычисляем восстановленный сигнал: x_rec = omp(A, y, k), параметры:

A — матрица измерений;
y — вектор измерений;
k — разреженность.

Функция возвращает x_rec длины n с не более чем k ненулевыми компонентами. Далее сравниваем x_rec с x_true.

x_rec = omp(A, y, k);

После этого вычисляем относительную ошибку восстановления, выводим её, число ненулевых элементов в исходном и восстановленном сигналах и проверяем совпадение их поддержек.

err = norm(x_true - x_rec) / norm(x_true)
println("Относительная ошибка восстановления: ", err)
println("Число ненулевых элементов в исходном сигнале: ", k)
println("Число ненулевых элементов в восстановленном: ",
        count(!iszero, x_rec))
println("Поддержка восстановлена верно: ",
        Set(findall(!iszero, x_true)) == Set(findall(!iszero, x_rec)))

Относительная ошибка восстановления: 3.772586544601344e-16
Число ненулевых элементов в исходном сигнале: 20
Число ненулевых элементов в восстановленном: 20
Поддержка восстановлена верно: true

Далее строим график исходного сигнала x_true и поверх него пунктиром красного цвета график восстановленного сигнала x_rec для их сравнения.

plot(x_true, label="Исходный", ylabel="Амплитуда",
          title="Сравнение сигналов", linewidth=2)
plot!(x_rec, label="Восстановленный",
      linewidth=2, linestyle=:dash, color=:red)

Вывод

В этом примере успешно решена задача сжатого зондирования: по неполным измерениям y = A·x (где m = 128 < n = 256) восстановлен разреженный сигнал x с k = 20 ненулевыми компонентами. Алгоритм OMP показал высокую точность:

относительная ошибка восстановления составила 3.77e-16, что близко к машинной точности;
число ненулевых элементов в исходном и восстановленном сигналах совпадает (20);
поддержка (позиции ненулевых элементов) восстановлена верно.

Таким образом, пример наглядно подтверждает основные положения теории сжатого зондирования: при достаточной разреженности сигнала его можно точно восстановить по существенно меньшему числу измерений, а жадный алгоритм OMP является эффективным инструментом для решения такой задачи.