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
