4

Есть график сигнала и точно такой же, сдвинутый по оси 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. Он имеет подозрительно большие значения, значительно превышающие исходныевведите сюда описание изображения

4
  • Почему вы считает, что в результате данных фокусов с преобразованием Фурье получится сдвиг фазы? Если использовали какую-то статью или учебник дайте ссылку. Сейчас выглядит просто как очень странный код, поэтому на ошибку в логике указать невозможно.
    – Chorkov
    24 мая 2023 в 14:40
  • Если вы ищите максиму свертки сигналов, то нужно брать максимум абсолютной величине комплексного значения, а не его действительной части.
    – Chorkov
    24 мая 2023 в 15:44
  • @Chorkov на выходе после обратного преобразования фурье ожидается значение корреляции двух векторов a и b, где по оси x цикличное смещение b относительно a. Но благодаря вашему замечанию я написал функцию вычисления конусного расстояния векторов и обнаружил, что выходной график во первых лежит в диапазоне от -1 до 1, а во вторых совсем немного отличается по форме и его пик находится в нужном месте, в отличии от сдвинутого результата вычисления кросс-корреляции при помощи кода выше. Я лишь пытался опираться на алгоритм быстрой кросс-корреляции. 24 мая 2023 в 17:32
  • @Chorkov почему я так считаю не знаю, но тут researchgate.net/post/How-to-perform-cross-correlation-by-FFT цитата 1) Put real signals in complex arrays --> real value as real part and set imaginary part to zero 2) compute complex FFT of both complex arrays 3) multiplication of complex FFT_1 and conjugate complex FFT_2 4) compute inverse FFT of multiplication 24 мая 2023 в 17:41

1 ответ 1

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];
}
1
  • но не забывайте, что длина результата должен быть не менее суммы длин обеих сигналов минус 1, для этого буфера сигналов удлиняют и заполняют нулями в конце 24 мая 2023 в 18:24

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge you have read our privacy policy.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.