Re: Hilbert transform

[Ben Abbott]

On Jul 6, 2012, at 10:19 PM, Sergei Steshenko wrote:

> ----- Original Message -----
>> From: Juan Pablo Carbajal <address@hidden>
>> To: Sergei Steshenko <address@hidden>
>> Cc: Ben Abbott <address@hidden>; "address@hidden" <address@hidden>
>> Sent: Saturday, July 7, 2012 1:49 AM
>> Subject: Re: Hilbert transform
>> 
> [snip]
>> 
>> You are working with a discrete transformation very sensitive to the
>> length of the signal. If you check with a longer signal everything is
>> fine.
>> 
>> t  = linspace (0,1,100);
>> y  = sin (2*pi*8*t) - sin (2*pi*10*t);
>> yg= imag ( hilbert ( imag ( hilbert (y))));
>> plot (t,y,t,-yg,'o')
>> 
>> Cheers
>> 
>> 
> 
> OK, the trouble continues :).
> 
> I am trying to run Hilbert transform in order to calculate a minimum phase 
> system phase response.
> 
> According to various sources, e.g. 
> https://en.wikipedia.org/wiki/Minimum_phase#Relationship_of_magnitude_response_to_phase_response
>  ,
> 
> phase(omega) = -Hilbert(log(abs(freq_response(omega))))
> .
> 
> Of course, the above is accurate up to an additive constant, e.g. the 
> transform can not know whether the system is inverting or not.
> 
> Here is a simple piece of code using Butterworth filter as a test case (see 
> it also in the attachment):
> 
> #-----------------------------------------------------------------------------
> N = 1000;
> 
> 
> frequency_range = 0:N - 1;
> zfrequency_range = exp(i * pi * frequency_range / N);
> 
> 
> lower_zcutoff = 0.01;
> 
> zorder = 2;
> figure_number = 1;
> 
> 
> [bfzeros, bfpoles, bfgain] = butter(zorder, lower_zcutoff);
> 
> znum = ones(1, length(frequency_range));
> for zero_number = 1:length(bfzeros)
>   znum = znum .* (bfzeros(zero_number) - zfrequency_range);
> endfor
> 
> zdenom = ones(1, length(frequency_range));
> for pole_number = 1:length(bfpoles)
>   zdenom = zdenom .* (bfpoles(pole_number) - zfrequency_range);
> endfor
> 
> frequency_response = bfgain * znum ./ zdenom;
> 
> figure(figure_number++);
> semilogx(frequency_range(2:end), abs(frequency_response(2:end)));
> title("magnitude response");
> 
> min_phase = -imag(hilbert(log(abs(frequency_response))));
> 
> 
> figure(figure_number++);
> semilogx(frequency_range(2:end), 180 * unwrap(arg(frequency_response(2:end)), 
> pi) / pi);
> title("measured phase response in degrees");
> 
> figure(figure_number++);
> semilogx(frequency_range(2:end), 180 * unwrap(min_phase(2:end), pi) / pi);
> title("min phase response in degrees");
> #============================================================================
> .
> 
> Both magnitude response and measured phase response look as expected for a 
> 2-nd order Butterworth LPF, but its minimum phase reconstituted through 
> Hilbert transform looks really wacky.
> 
> Increasing N makes it look even whackier.
> 
> I performed this test because minimum phase reconstituted through Hilbert 
> transform for my real experimental data looks whacky, though, admittedly, 
> less whacky than in this simple test case.
> 
> Any ideas ? Can anybody run this on both Octave and Matlab ?
> 
> Thanks,
>   Sergei.<test_butterworth_plus_hilbert.m>

Sergei,

The hilbert function relies upon the fft(). Which implies the signal being 
transformed should be periodic.

You'll need to add the negative frequencies to the frequency_response.

        frequency_response = [frequency_response, 
fliplr(frequency_response(2:end))]

And then use the Hilbert transform to recover the phase information.

Ben