fakedad - 2 years ago 191
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

`hilbert()`
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
`hilbert()`
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;
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
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download