Вычислительный аналитический сигнал с использованием FFT в C++

Question

В настоящее время я потерял попытку выяснить, как реализовать эквивалентную версию функции MATLAB hilbert() в C++. Я очень новичок в обработке сигналов, но, в конечном счете, я хотел бы выяснить способ фазового сдвига любого заданного сигнала на 90 градусов. Я пытался следовать методу, предложенному в this question on MATLAB central, который, по-видимому, работает на основе тестов с использованием GNU Octave.

У меня есть то, что я считаю рабочим выполнением как БПФ, так и обратного БПФ, и я попытался реализовать метод, описанный в ответе this post, чтобы вычислить аналитический сигнал. Я попытался сделать это, применив БПФ, установив верхнюю половину массива в ноль, а затем применив обратный БПФ, но, основываясь на графиках, которые я сделал на выходе из теста, должна быть проблема с тем, как я осуществил поиск аналитического сигнала.

Что было бы подходящим способом реализации функции hilbert() от MATLAB в C++ с учетом рабочей реализации БПФ и обратного БПФ? Есть ли лучший способ добиться сдвига фазы на 90 градусов?

Answer 1

При проверке реализации MATLAB следующее должно возвращать тот же результат, что и функция hilbert. Очевидно, вам придется изменить его, чтобы он соответствовал вашей конкретной реализации. Я предполагаю, что существует класс signal.

signal hilbert(const signal &x) 
{ 
    int limit1, limit2; 
    signal xfreq = fft(x); 
    if (x.numel % 2 == 0) { 
     limit1 = x.numel/2; 
     limit2 = limit1 + 1; 
    } else { 
     limit1 = (x.numel + 1)/2; 
     limit2 = limit1; 
    } 
    // multiply the first half by 2 (except the first element) 
    for (int i = 1; i < limit1; ++i) { 
     xfreq[i].real *= 2; 
     xfreq[i].imag *= 2; 
    } 
    for (int i = limit2; i < x.numel; ++i) { 
     xfreq[i].real = 0; 
     xfreq[i].imag = 0; 
    } 
    return ifft(xfreq); 
} 

Редактировать: Забыл установить вторую половину в нули. Edit2: Исправлена ​​логическая ошибка. Я закодировал следующее в MATLAB, которое соответствует гильберту.

function h = hil(x) 
    n = numel(x); 
    if (mod(n,2) == 0) 
     limit1 = n/2; 
     limit2 = limit1 + 2; 
    else 
     limit1 = (n+1)/2; 
     limit2 = limit1+1; 
    end 

    xfreq = fft(x); 

    for i = 2:limit1 
     xfreq(i) = xfreq(i)*2; 
    end 
    for i = limit2:n 
     xfreq(i) = 0; 
    end 

    h = ifft(xfreq); 
end