Kramers Kronig versus Hilbert

Dear Matlab community,

A quick question about Kramers-Kronig.

I am aware that KK is used extensively in optics, but here the application is slightly different. I measure a complex impedance over a frequency range. Now I'd like to run the real part as well as the imaginary part of that impedance through KK and compare the calculate the error (residuals).

This is done to see if the measurement was in equilibrium and the values measured are sensible. There are some parts of code out there that I've tried but they all seem to fail at some point. I found that there is the built-in function 'hilbert' that might just do the trick with some adjustment, but it is exactly these adjustments, that I do not know.

As the original data I have an array with col1, col2, col3 being frequency, real part of the impedance, and imaginary part of the impedance, respectively.

Here is the KK result as plotted from a commercially available software ...

When running my code I get:

blue: KK

red: original data

I tend to believe that the commercial SW is working properly :-)

If that is true, then the overall shape of my results looks alright, but there seems to be a significant shift in the x-coordinate... (to lower values)

The two pieces of code are shown below.

Input values for the calculating the imKK value are:

the real part of the impedance, re

the angular frequency, omega (which is: freq. * 2 * pi())

Input values for the calculating the reKK value are:

the NEGATIVE imaginary part of the impedance, im (as over most of the frequency range ImZ<0)

the angular frequency, omega (which is: freq. * 2 * pi())

The imKK calc is done like this:

++++++++++++++++++++++++++++

function imkk = kkre2im(re,omega)

% Perform Kramers-Kronig transform of real data re and

% angular frequencies omega to get imaginary parts. Handle

% division by zero by letting integrand be the mean

% of the adjoining values.

imkk = zeros(size(re)); % set up vector for reconstructed imag parts

N = length(re); % number of data points

warning off Matlab:divideByZero % avoid warning, since division

% for each data point, do the transform

for k = 1:N;

rek = re(k);

omegak = omega(k);

% compute the integrand

inte = (re - rek)./(omega.^2 - omegak.^2);

% handle division by zero by taking mean and end-points

% by setting equal to adjoining point

switch k

case 1

inte(k) = inte(2);

case N

inte(k) = inte(N-1);

otherwise

inte(k) = ( inte(k-1) + inte(k+1) )/2;

end

% integrate

%intecur = integrate(omega,inte);

intecur = integrate(log10(omega),inte.*omega*log(10));

% scale by constant to get final imaginary parts

imkk(k) = 2*omegak/pi*intecur;

end

warning on Matlab:divideByZero

++++++++++++++++++++++++++++

++++++++++++++++++++++++++++

The reKK calc is done like this:

++++++++++++++++++++++++++++

function rekk = kkim2re(im,omega)

% Perform Kramers-Kronig transform of imaginary data im and

% angular frequencies omega to get real parts. Handle

% division by zero by letting integrand be the mean

% of the adjoining values.

rekk = zeros(size(im)); % set up vector for reconstructed real parts

N = length(im); % number of data points

warning off Matlab:divideByZero % avoid warning, since division

% for each data point, do the transform

for k = 1:N;

imk = im(k);

omegak = omega(k);

% compute the integrand

inte = (omega.*im - omegak*imk)./(omega.^2 - omegak.^2);

% handle division by zero by taking mean and end-points

% by setting equal to adjoining point

switch k

case 1

inte(k) = inte(2);

case N

inte(k) = inte(N-1);

otherwise

inte(k) = ( inte(k-1) + inte(k+1) )/2;

end

% integrate

intecur = integrate(log10(omega),inte.*omega*log(10));

% scale by constant to get final real parts

rekk(k) = -2/pi*intecur;

end

warning on Matlab:divideByZero

++++++++++++++++++++++++++++

++++++++++++++++++++++++++++

Any suggestions as to why the results are different and (in case you also believe the results shown by the commercial software) what am I doing wrong?

The original data is like this (Again, with col1, col2, col3 being frequency, real part of the impedance, and imaginary part of the impedance, respectively):

10078.1250000000 0.00130600000000000 0.00293600000000000

8296.87500000000 0.00128000000000000 0.00241600000000000

6755.51464800000 0.00125700000000000 0.00196100000000000

5671.87500000000 0.00124500000000000 0.00163200000000000

4641.54394500000 0.00123400000000000 0.00131700000000000

3814.33813500000 0.00121800000000000 0.00106100000000000

3170.95581100000 0.00120900000000000 0.000850000000000000

2619.48535200000 0.00120900000000000 0.000663000000000000

2159.92651400000 0.00118400000000000 0.000476000000000000

1792.27941900000 0.00123500000000000 0.000333000000000000

1459.31604000000 0.00125700000000000 0.000168000000000000

1216.94714400000 0.00128500000000000 9.50000000000000e-05

998.263855000000 0.00132900000000000 -1.70000000000000e-05

824.652771000000 0.00133600000000000 -9.50000000000000e-05

676.081726000000 0.00138800000000000 -0.000189000000000000

564.236084000000 0.00142800000000000 -0.000267000000000000

460.379456000000 0.00147600000000000 -0.000349000000000000

383.522705000000 0.00152400000000000 -0.000422000000000000

315.504822000000 0.00158000000000000 -0.000500000000000000

260.416656000000 0.00164300000000000 -0.000581000000000000

217.013885000000 0.00171000000000000 -0.000662000000000000

177.556824000000 0.00179700000000000 -0.000755000000000000

146.484375000000 0.00190100000000000 -0.000845000000000000

121.228455000000 0.00201500000000000 -0.000938000000000000

100.446426000000 0.00214900000000000 -0.00102400000000000

82.7205890000000 0.00230000000000000 -0.00110900000000000

68.2645650000000 0.00247300000000000 -0.00116800000000000

56.2500000000000 0.00265000000000000 -0.00122600000000000

46.8750000000000 0.00282800000000000 -0.00126100000000000

38.4221310000000 0.00301600000000000 -0.00128000000000000

31.6722980000000 0.00320600000000000 -0.00130800000000000

26.0416680000000 0.00339900000000000 -0.00132300000000000

21.7013890000000 0.00359200000000000 -0.00131800000000000

17.7556820000000 0.00376700000000000 -0.00129100000000000

14.7405660000000 0.00393000000000000 -0.00127200000000000

12.2070310000000 0.00409400000000000 -0.00122200000000000

9.93114500000000 0.00426600000000000 -0.00116900000000000

8.22368400000000 0.00440700000000000 -0.00110700000000000

6.85307000000000 0.00446700000000000 -0.00105100000000000

5.63401500000000 0.00462400000000000 -0.00100000000000000

4.68750000000000 0.00464700000000000 -0.000957000000000000

3.82965700000000 0.00471500000000000 -0.00101000000000000

3.15869300000000 0.00470100000000000 -0.00113200000000000

2.60127600000000 0.00484600000000000 -0.00126300000000000

2.14629100000000 0.00507800000000000 -0.00141500000000000

1.76753400000000 0.00525200000000000 -0.00152900000000000

1.45393900000000 0.00559300000000000 -0.00157000000000000

1.20936500000000 0.00589200000000000 -0.00159500000000000

0.999041000000000 0.00619300000000000 -0.00156500000000000

0.820641000000000 0.00652000000000000 -0.00129400000000000

0.675822000000000 0.00654400000000000 -0.00134600000000000

0.564759000000000 0.00680100000000000 -0.00107200000000000

0.468750000000000 0.00689900000000000 -0.000870000000000000

0.384221000000000 0.00717100000000000 -0.000851000000000000

0.316723000000000 0.00713500000000000 -0.000612000000000000

0.261872000000000 0.00715500000000000 -0.000508000000000000

0.216014000000000 0.00713200000000000 -0.000361000000000000

0.178232000000000 0.00721300000000000 -0.000297000000000000

0.146944000000000 0.00732800000000000 -0.000258000000000000

0.121438000000000 0.00730000000000000 -0.000120000000000000

0.100160000000000 0.00726600000000000 -7.70000000000000e-05

Answer 1

[Bruno Melo]

I cannot run your code. Not sure if you posted the whole code.

There is a very good tool which you can use for this purpose which is the LIN-KK. It is freely available as a standalone app and to my best knowledge it is one of the best out there.

I would also like to find one MATLAB code for KK transformations applied to impedance spectroscopy, but all codes I found did not work as good as LIN-KK

Hope it helps.

[Dmitry Kolosnitsyn]

How to evaluate the correctness of the algorithm?

The fact is that I processed the same impedance data using LIN-KK and using commercial software, the screenshots of which were given by the author (this is EC-Lab).

And the resulting data is very different.

The question is who to trust

Answer 2

[Michael Tsuk]

Doing a numerical integral of a function with singularities, such as the Hilbert kernel, requires care. I'm not sure what your "integrate" function is doing, but it doesn't look like you have any way to tell it about the singularities. It wouldn't surprise me if that's what's leading to your offset.

You'd have a different challenge using the "hilbert" function: it does the discrete Hilbert transform, which means it assumes data is on a uniform frequency grid, starting from zero. You'd have to interpolate your data onto such a grid, and the easy ways of doing interpolation introduce causality artifacts that might swamp the information you're trying to extract.

I'm not familiar with the LIN-KK tool that Bruno Melo recommended, but I looked at the webpage, and it seems to me that they understand the issues involved. So I'd be inclined to trust LIN-KK. And it looks pretty interoperable with MATLAB.