fakedad fakedad - 4 months ago 61
C++ Question

Computing analytical signal using FFT in C++

I'm currently lost trying to figure out how to implement an equivalent version of MATLAB's

function in C++. I'm very new to signal processing, but, ultimately, I would like to figure out a way to phase shift any given signal by 90 degrees. I was attempting to follow the method suggested in this question on MATLAB central, which appears to work based on tests using GNU Octave.

I have what I believe to be a working implementation of both FFT and the inverse FFT, and I have tried implementing the method described in the answer to this post in order to compute the analytical signal. I have tried doing this by applying the FFT, setting the upper half of the array to zero, and then applying the inverse FFT, but, based on graphs I made of output from a test, there must be a problem with the way I have implemented finding the analytical signal.

What would be a suitable way to implement the
function from MATLAB in C++ given a working implementation of FFT and inverse FFT? Is there a better way of achieving the 90 degree phase shift?


Checking the MATLAB implementation the following should return the same result as the hilbert function. You'll obviously have to modify it to match your specific implementation. I'm assuming a signal class of some sort exists.

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);

Edit: Forgot to set the second half to zeros. Edit2: Fixed a logical error. I coded the following up in MATLAB which matches hilbert.

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

    xfreq = fft(x);

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

    h = ifft(xfreq);