Re: Hilbert transform

[John B. Thoo]

On Jul 7, 2012, at 7:01 AM, John B. Thoo wrote:

> here is another routine for the Hilbert transform that avoids using "imag".
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin hilb.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> function out = hilb(in)
> % HILBERT computes real hilbert transform of real input
> % Has symbol -i sgn(k)
> lx = length (in);
> in = fft (in);
> in (1:(lx/2+1)) = -i*in (1:(lx/2+1)); 
> in (lx/2+2:lx) = i*in (lx/2+2:lx); 
> out = real (ifft(in));
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%% begin hilb.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Hi,  When I tried to use this with Sergei's latest example, replacing

min_phase = unwrap(-imag(hilbert(log(abs([frequency_response, 
fliplr(frequency_response(2:end))])))), pi); # thanks to Ben Abbot for 'fliplr' 
part

with

min_phase = unwrap(-hilb(log(abs([frequency_response, 
fliplr(frequency_response(2:end))]))), pi); # thanks to Ben Abbot for 'fliplr' 
part

I got the error

octave-3.2.3:5> test_butterworth_plus_hilbert
error: subscript indices must be either positive integers or logicals.
error: called from:
error:   /Users/jbthoo/Desktop/hilb.m at line 7, column 16
error: evaluating argument list element number 1
error: evaluating argument list element number 1
error:   /Users/jbthoo/Desktop/test_butterworth_plus_hilbert.m at line 35, 
column 11
octave-3.2.3:5>

Modifying the routine above to

%%%%%%%%%%%%%%%%%%%%%%%%%%% begin fixed hilb.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = hilb (in)
% HILBERT computes real hilbert transform of real input
% Has symbol -i sgn(k)
lx = length (in);
in = fft (in);
in (1:floor (lx/2)+1) = -i*in (1:floor (lx/2)+1); 
in (floor( lx/2)+2:lx) = i*in (floor (lx/2)+2:lx); 
out = real (ifft (in));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% end fixed hilb.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

then I think I get the same results.  I'm sorry for taking up bandwidth on this 
list, but I wanted to fix that.  Thanks for your patience.

---John.


-------------------------cut here: Sergie's example-------------------------
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");


phase_response = unwrap(arg(frequency_response), pi);

min_phase = unwrap(-imag(hilbert(log(abs([frequency_response, 
fliplr(frequency_response(2:end))])))), pi); # thanks to Ben Abbot for 'fliplr' 
part


figure(figure_number++);
semilogx(frequency_range(2:end), 180 * phase_response(2:end) / pi);
title("measured phase response in degrees");

figure(figure_number++);
semilogx(frequency_range(2:end), 180 * min_phase(2:numel(frequency_range)) / 
pi);
title("min phase response in degrees");

figure(figure_number++);
semilogx(frequency_range(2:end), 180 * (phase_response(2:end) - 
min_phase(2:numel(frequency_range))) / pi);
title("phase responses difference in degrees");