Is it possible to implement a block-wise Hilbert transformer using FFT

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


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.

• Did you implement this as overlap add ? May 3, 2022 at 10:36
• @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. May 3, 2022 at 11:38
• 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). May 3, 2022 at 16:03
• @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? May 4, 2022 at 0:42
• @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. May 4, 2022 at 0:49

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

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:

• 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. May 3, 2022 at 11:35
• @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? May 3, 2022 at 11:39
• @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. May 3, 2022 at 11:42
• Please see the edit. May 3, 2022 at 11:44
• y0 is not equal to y1 in my example even when I use MATLAB's hilbert. That's my main problem. May 3, 2022 at 11:49