1
$\begingroup$

I am processing audio and trying to obtain a 90º phase shift when processing block by block.

I have implemented a Hilbert Transformer (from Lyon's 3rd. Ed. & MATLAB). MATLAB code below (I believe this is similar to what MATLAB hilbert does).

function result = my_hilbert(s)
    fprintf("Using My Hilbert\n");
    n = length(s)

    %fftLength = 2 ^ nextpow2(2 * n - 1);
    fftLength = n

    H = zeros(1, fftLength);
    H(1) = 1;
    H(2:n / 2) = 2;
    H(n / 2 + 1) = 1;

    x = zeros(1, fftLength);
    x(1:n) = s
    X = fft(x)

    Y = X .* H
    y = ifft(Y)

    result = y(1:n);
end

This implementation seems to assume one is processing the complete audio stream in a single call.

I need to be able to achieve a phase shift but will be processing one block (block size between 1024 and 2048 samples) at a time.

When I simply call this implementation on each block I receive audible pops around the block boundaries.

I have considered overlap+add but when using an FFT length other than length of the input, the results are no longer the same as MATLAB's hilbert() function.

How can one obtain a 90º phase shift when processing block by block without these discontinuities?

$\endgroup$
  • $\begingroup$ I think this answer could be helpful. $\endgroup$ – Matt L. Mar 18 at 21:58
  • $\begingroup$ What should be the length of the kernel? Same length as the input? So, should FFT length be 2 * length(s) - 1 rounded up to the next power of 2? I think when I did this the overlap was all 0. If you're both saying that length H should be the same as s then I'll revisit the overlap to confirm all 0's or not. $\endgroup$ – Random Mar 18 at 23:06
  • $\begingroup$ @MattL. 's "$N$" is the same as my "$L$" (the length of the FIR) and Matt's "$L$" is the same as my "$B$" (the number of samples processed per block). i like to, whenever an FFT or DFT is being used, reserve the term "$N$" for the size of the DFT. $\endgroup$ – robert bristow-johnson Mar 18 at 23:11
  • $\begingroup$ @robertbristow-johnson this makes perfect sense to me - I have implemented other code to do fast convolutions and block processing and it works. However, what I'm struggling with here is that using any DFT length (N in your terms) other than the length of the input and the results are no longer the same as MATLAB's hilbert() function. $\endgroup$ – Random Mar 18 at 23:20
1
$\begingroup$

The key to block (or stream) filtering in Matlab is using the filter's command ability to return and reuse its internal state.

To illustrate, first create a filter that implements the Hilbert transform.

hilb = designfilt('hilbertfir','FilterOrder',60, ...
       'TransitionWidth',400,'SampleRate',1e4);

This example is taken directly from the documentation (type doc hilbert in Matlab).

Now, we need this form of the filter command (see doc filter):

[y, zo] = filter(b,a,x,zi)

which uses initial conditions zi and returns final conditions in zf.

To put everything together, say you have incoming data and store it in variable block. Then, you can do something like:

block = acquire_block_of_data();
[output, zo] = filter(hilb, 1, block);
while 1
    block = acquire_block_of_data();
    zi = zo;
    [output, zo] = filter(hilb, 1, block, zi)
end

Note that we need an initial call to filter without state; alternatively, you can set let zi have an initial value equal to a vector of zeros with the same length as the filter order.

$\endgroup$
  • 2
    $\begingroup$ and the primary issue is still the old overlap-add or overlap-save. when you use the FFT or DFT to filter samples in blocks, your filter must be an FIR, and the the H[k] must be the DFT of an FIR h[n]. and you must follow the rules about overlap-add or overlap-save which is that the size of the FFT, $N$, must be (virtually) as large as the sum of the length of the FIR $L$ and the number of samples, $B$, processed in each block. $$ N \ge L+B-1 $$ you must do that (and the overlap-adding or overlap-saving) to avoid the clicks between blocks. $\endgroup$ – robert bristow-johnson Mar 18 at 19:30
  • 1
    $\begingroup$ since the Hilbert transformer has an impulse response that is not FIR, you must approximate it with an FIR. the transfer function H[k] cannot be that simply of a Hilbert transformer. $\endgroup$ – robert bristow-johnson Mar 18 at 19:33
  • $\begingroup$ @robertbristow-johnson Good points -- I focused just on the Matlab implementation. $\endgroup$ – MBaz Mar 18 at 20:40
  • 1
    $\begingroup$ it's not the length of the frequency response, H[k], that you need to worry about, it's the length of the impulse response h[n] that you need to worry about. the impulse response must be finite in length (hence an "FIR") and that finite length i am calling "$L$". the length of (linearly) convolving a block of samples of length $B$ with an FIR of length $L$ is $B+L-1$. make sure your FFT is at least that big. your FIR of length $L$ will be padded with $N-L$ zeros and, for overlap-add, your block of $B$ samples is padded with $N-B$ zeros. (overlap-save does not pad this block). $\endgroup$ – robert bristow-johnson Mar 18 at 21:46
  • 3
    $\begingroup$ well MATLAB's hilbert() is not for a continuously-running filter in real time. if what you have is a finite-sized signal and you want to filter that whole thing without splitting into blocks, MATLAB's hilbert()is the way to do it. but if you want to break apart a very long signal into blocks and filter the blocks and then stitch together the results, that is what we do with fast convolution. then your practical Hilbert transformer is an FIR approximation to the ideal Hilbert that is delayed a little to make the FIR causal. so the results won't be exactly the same. $\endgroup$ – robert bristow-johnson Mar 19 at 2:17

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.