Re: Hilbert transform

[Sergei Steshenko]

----- 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
Description: Text Data