Re: Hilbert transform

[Sergei Steshenko]

----- Original Message -----
> From: Ben Abbott <address@hidden>
> To: Sergei Steshenko <address@hidden>
> Cc: Juan Pablo Carbajal <address@hidden>; "address@hidden" <address@hidden>
> Sent: Saturday, July 7, 2012 7:11 AM
> Subject: Re: Hilbert transform
> 
[snip]
> 
> 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
>

Thanks a lot, that _really_ helps with the Butterworth filter example.

However, it makes things worse with my experimental data. The data may be 
problematic - it intentionally has measured data up to _half_ Nyquist frequency 
(sample_rate/4), so data above it is noise. Maybe this noise baffles the 
transform.

Anyway, please find modified code below and also attached - I think it can be a 
useful demo/learning example.

Regards,
  Sergei.


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");

test_butterworth_plus_hilbert.m
Description: Text Data