2
$\begingroup$

I'm implementing a discrete Hilbert transformer and I know that an ideal Hilbert transformer is anti-causal and has infinite length so we can only make approximation. There are some FIR and IIR implementations, while MATLAB's hilbert() uses FFT to modify the spectrum of input.

What I've done is to keep the DC and Nyquist components, double the components between them and make the negative frequencies zero. But I got wrong results when processing signal block by block. A MATLAB example using windowing and overlap add is given below.

fs = 48e3;
Ts = 1/fs;
t = 0:Ts:1;
w = 2*pi*10;
sig = sin(w*t).';
N = length(sig);
y0 = hilbert(sig); % true result using full-length FFT
y1 = zeros(N, 1); % block-wise processing

blockSize = 1024;
hopSize = blockSize / 2;
startIdx = 1;
endIdx = blockSize;
win = hann(blockSize, 'periodic');

while endIdx < N
    x = sig(startIdx:endIdx) .* win; % windowing
    xa = hilbert(x);
    y1(startIdx:endIdx) = xa + y1(startIdx:endIdx); % overlap add
    startIdx = startIdx + hopSize; % 50% overlap
    endIdx = endIdx + hopSize;
end

figure; 
subplot(211)
plot(t, real(y0), t, imag(y0)); legend('original', 'hilbert')
subplot(212)
plot(t, real(y1), t, imag(y1)); legend('original', 'hilbert block-wise')

enter image description here

FIR and IIR filters have filter states saved in the memory which deal with this case, but I don't know how to solve it.

$\endgroup$
6
  • 1
    $\begingroup$ Did you implement this as overlap add ? $\endgroup$
    – Hilmar
    May 3 at 10:36
  • $\begingroup$ @Hilmar No. I take an FFT of N point input, modify the FFT results and do an IFFT, I got N point output for each block. It's not a convolution so I don't get a longer output to overlap add. $\endgroup$
    – DSP novice
    May 3 at 11:38
  • $\begingroup$ A modification in the frequency domain is a convolution in the time domain. So you will need to overlap (by the length of the impulse response that is above your desired noise floor or error level). $\endgroup$
    – hotpaw2
    May 3 at 16:03
  • $\begingroup$ @hotpaw2 As far as I know, overlap-add is used in two way -- one is for fast convolution using FFT, which is an LTI system, another one is for time-variant frequency-domain modification, usually required windowing to reduce frequency leakage. In my current code I used OLA in the second way. Do you mean that I should treat it as the first way, just like a convolution? $\endgroup$
    – DSP novice
    May 4 at 0:42
  • $\begingroup$ @hotpaw2 When I do y0 = hilbert(sig); the default FFT point is equal to the length of sig, while calculating linear convolution using FFT requires a FFT length of N+M-1, where N and M are the lengths of sig and the unknown impulse response respectively. So I don't think this process is a convolution. $\endgroup$
    – DSP novice
    May 4 at 0:49

1 Answer 1

1
$\begingroup$

If the issue is in implementing the Hilbert block by block in either case (using MATLAB's hilbert or own method) then this would properly be done using overlap-add or overlap-save block processing techniques. Here in this case the combination of windowing with the overlap and add (as commonly done) reduces the edges of greatest distortion in the Hibert implementations. This is implemented by following the pseudo-code in the above link for the implementation of linear convolution by replacing the FIR impulse response in that link with the Hilbert impulse response as described further below. The windowing I mention would be the windowing of the ideal time delayed and truncated Hilbert impulse response as would properly be done in the implementation of an FIR Hilbert using the windowing design approach.

If the issue is with replicating the Matlab hilbert implementation with FFT processing, below is an FFT implementation of the Hilbert in MATLAB that I confirmed produces an identical result as the hilbert function:

N = length(sig);
sig_spectrum = fft(sig);
hilbert_spectrum = zeros(1, N);
hilbert_spectrum(1) = sig_spectrum(1);
hilbert_spectrum(2: ceil(N/2) - 1) = 2 * sig_spectrum(2:ceil(N/2) - 1);
if mod(N, 2) == 0
  hilbert_spectrum(N/2) = sig_spectrum(N/2);
endif
hilbert_time = ifft(hilbert_spectrum);

Another approach is to simply implement the filter directly as an FIR filter since we know the coefficients are given by the impulse response, which for the Hilbert is given as:

$$h(t) = \frac{1}{\pi t}$$

As the OP mentioned this is non-causal and thus for implementation as FIR filter coefficients, the above impulse response is sampled, delayed and truncated (and windowed to reduce sidelobes).

Hilbert Implementation

One "trick" if we already have a half-band filter in our pocket, is we can easily convert it to a Hilbert by taking the absolute value of the coefficients, setting the largest center tap to 0 and negating all the taps left of center:

half band to hilbert

$\endgroup$
13
  • $\begingroup$ Yes but sorry Dan, I got the same result just as hilbert when taking the full-length FFT for the long input signal. But things went wrong when I do it block by block. $\endgroup$
    – DSP novice
    May 3 at 11:35
  • $\begingroup$ @DSPnovice What do you mean "Block by Block"--- you take the Hilbert of a block and then you do it with the FFT approach of that same block and you don't get the same result? What is the FFT in that case that is different from a "full-length FFT"? Isn't this the same thing-- just the total number of samples changes? Note the careful handling of the Nyquist bin depending on if the block length is even or odd, is that perhaps your issue? $\endgroup$ May 3 at 11:39
  • $\begingroup$ @DSPnovice Please try my code to see if you get a different result for your smaller block-- if not, maybe if you provide an example in your question the issue might be clearer. $\endgroup$ May 3 at 11:42
  • $\begingroup$ Please see the edit. $\endgroup$
    – DSP novice
    May 3 at 11:44
  • $\begingroup$ y0 is not equal to y1 in my example even when I use MATLAB's hilbert. That's my main problem. $\endgroup$
    – DSP novice
    May 3 at 11:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.