Нахождение смещения двух графиков (кросс-корреляция)
Есть график сигнала и точно такой же, сдвинутый по оси 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);
}
}
А также в экселе я построил графики для наглядности
b должно сместится к a полностью, но оно смещается на 94 правее, причём оригинал может быть смещён на любую величину, результат будет тот же для этой формы сигнала. Я бы мог это устранить костылём, подав в signal_c произведение одного из двух похожих сигналов и узнав, насколько он соврёт, но хотелось бы узнать причину, что я делаю не так.
P.S. Хочу добавить график кросс-корреляции, который я получаю после обратного преобразования в signal_c. Он имеет подозрительно большие значения, значительно превышающие исходные