Нахождение смещения двух графиков (кросс-корреляция)

Рейтинг: 4Ответов: 1Опубликовано: 24.05.2023

Есть график сигнала и точно такой же, сдвинутый по оси x на указанную в коде величину. Мне нужно найти насколько один график смещён относительно другого. Дело в том, что по какой-то магической причине вычисляемый shift всегда на какое-то число больше, чем заданный кодом. Причём, если вычислить shift при нулевом смещении, условно пусть будет 94, то теперь отнимая 93 от shift при любом смещении мы будем получать изначально заданную величину точно. Т.е. отклонение (ошибка) вычисляемого shift не зависит от изначально заданного смещения. Если сделать N в 2 раза больше, то и ошибка станет ровно в 2 раза больше. А также ошибка хаотичным образом зависит от формы изначального сигнала. С чем это связано? На всякий случай скину много кода.

const int N = 4096;
std::vector<Complex> signal_a(N);
std::vector<Complex> signal_b(N);
std::vector<Complex> signal_c(N);

// ... тут код заполнения signal_a

for (int i = 0; i < N; i++) // заполняем signal_b со смещением относительно signal_a
    signal_b[(i + 500) % signal_b.size()] = signal_a[i];

fft(signal_a);
fft(signal_b); 

for (int i = 0; i < N; i++) { // комплексное произведение
    signal_c[i] = signal_a[i] * signal_b[i];
}

ifft(signal_c); // обратное БПФ

auto max_it = std::max_element(signal_c.begin(), signal_c.end(), [](Complex a, Complex b) {
    return a.real < b.real;
}); // находим итератор на максимальный в массиве элемент

int shift = std::min(
    std::distance(signal_c.begin(), max_it),
    -std::distance(signal_c.end(), max_it)
); // вычисляем индекс максимального элемента относительно края
// БПФ без рекурсии
void fft(std::vector<Complex>& x, bool inverse) {
    int n = x.size();

    std::vector<int> rev(n);
    for (int i = 0; i < n; ++i) {
        rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0);
        if (rev[i] < i) {
            std::swap(x[i], x[rev[i]]);
        }
    }

    for (int k = 2; k <= n; k <<= 1) {
        double angle = 2 * PI / k * (inverse ? 1 : -1);
        Complex wn(cos(angle), sin(angle));
        for (int i = 0; i < n; i += k) {
            Complex w(1);
            for (int j = 0; j < k / 2; ++j) {
                Complex u = x[i + j];
                Complex t = w * x[i + j + k / 2];
                x[i + j] = u + t;
                x[i + j + k / 2] = u - t;
                w = w * wn;
            }
        }
    }

    if (inverse) {
        for (int i = 0; i < n; ++i) {
            x[i] = Complex(x[i].real / n, x[i].imag / n);
        }
    }
}
// прямое
void fft(std::vector<Complex>& x) {
    fft(x, false);
}
// обратное
void ifft(std::vector<Complex>& x) {
    fft(x, true);
    for (int i = 0; i < x.size(); i++) {
        x[i] = Complex(x[i].real, -x[i].imag);
    }
}

А также в экселе я построил графики для наглядности 1

b должно сместится к a полностью, но оно смещается на 94 правее, причём оригинал может быть смещён на любую величину, результат будет тот же для этой формы сигнала. Я бы мог это устранить костылём, подав в signal_c произведение одного из двух похожих сигналов и узнав, насколько он соврёт, но хотелось бы узнать причину, что я делаю не так. P.S. Хочу добавить график кросс-корреляции, который я получаю после обратного преобразования в signal_c. Он имеет подозрительно большие значения, значительно превышающие исходныевведите сюда описание изображения

Ответы

▲ 1

Ошибка в том, что перед вычислением произведения рядов фурье нужно вычислить комплексное сопряжение, для std::Complex это действие выполняет std::conj, иначе можно записать как Complex(real, -imag).

for (int i = 0; i < N; i++) {
    signal_c[i] = Complex(signal_a[i].real, -signal_a[i].imag) * signal_b[i];
}