信息处理 - 函数的双指数（双指数）卷积 - 吾爱随笔录

函数的双指数（双指数）卷积

信息处理离散信号卷积反卷积

2022-01-31 09:56:04

概括

我正在尝试对最初从反卷积计算的一些数据进行卷积（反之亦然）。但是我没有得到预期的图表。

蓝色是预期的，红色是预期的插值版本。然后菱形线是各种卷积，两个半衰期中的一个或两个都在卷积中活跃。

问题

我是否正确使用了双指数公式？我认为是这样，因为它与所绘制的两个斜率的书相匹配。
ISR 是根据胰岛素浓度的反卷积计算的，然后在 ISR 上运行卷积应该给出胰岛素浓度，对吗？我相信这是正确的，因为“计算卷积运算的逆运算称为反卷积”。来自https://en.wikipedia.org/wiki/Convolution
我是否正确使用了卷积函数？我假设不是。这是函数“窗口”的问题，试图将其应用于整个 Y 集，而不是一次将其应用于 5 分钟的 ISR 数据？

数据

我有以下数据是从一篇论文的（https://doi.org/10.2337/diabetes.51.2007.S258）图表中提取的（使用https://apps.automeris.io/wpd/）。

数据是：

实验收集的以 mU/ml 为单位的胰岛素浓度。
从第一个数据计算的以 mU/ml/min 为单位的胰岛素分泌率 (ISR)。

图来自https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only

ISR是根据论文计算的，

一种基于胰岛素双指数消失率的反卷积方法，假设胰岛素的半衰期为 2.8 和 5 分钟，其中部分缓慢成分为 28% (19)

参考文献 19 是https://pubmed.ncbi.nlm.nih.gov/11134098/

和图2的描述

ISR 通过对血浆胰岛素浓度的反卷积分析进行估计。

因此，因为 ISR 是根据反卷积计算的，所以在 ISR上运行卷积应该给出胰岛素浓度，对吗？

公式

使用https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4似乎公式是（其他来源支持1、2、3和第 40 页公式 4 和 5）

y = a e^{- b t} + c e^{- d t}

$y=ae^{-bt} + ce^{-dt}$

所以我插入了论文提到的具体值（我不是 100% 知道“分数慢分量”是什么，请参阅答案）。

注意：第一个 0.28 是分数慢分量，第二个是 2.8 分钟。

y = 72 % e^{- l o g (2) t / 5.0} + 28 % e^{- l o g (2) t / 2.8}

$y=72\%e^{-log(2)t/5.0} + 28\%e^{-log(2)t/2.8}$

所以

y = 0.72 e^{- l o g (2) t / 5.0} + 0.28 e^{- l o g (2) t / 2.8}

$y=0.72e^{-log(2)t/5.0} + 0.28e^{-log(2)t/2.8}$ 放在一起

首先，我想确保我理解了公式。所以我把它画在一个半逻辑图上。

这是基于 Joseph T. DiPiro 所著的临床药代动力学概念第 78 页上的图 6-9

所以看来我有正确的公式我认为

运行卷积

我得到以下三张图，第一张是原图。其次是从图中提取的ISR数据，最终的图有，

蓝色的原始胰岛素含量
红色插值版本到 110 个数据点
三个卷积具有两个双指数斜率或只有一个是红色的
以及用宝蓝色对原始数据进行卷积。

完整代码

import numpy as np
import pandas
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import splrep, splev
from scipy.optimize import curve_fit
import urllib.request
import scipy as sp

# set matplotlib display properties
import matplotlib as mpl
mpl.rcParams['lines.linewidth'] = 2
font = {'family' : 'normal',
        'weight' : 'bold',
        'size'   : 22}

mpl.rc('font', **font)


# data from figure 2
# https://diabetes.diabetesjournals.org/content/51/suppl_1/S258.figures-only

#data_ins = pandas.read_csv("Insulin.dat", header = None, delimiter = '\t')

#data_isr = pandas.read_csv("ISR.dat", header = None, delimiter = '\t')

insulin_json = '{"0":{"0":4.143,"1":13.954,"2":23.984,"3":34.014,"4":44.044,"5":54.073,"6":64.103,"7":74.133,"8":83.944,"9":93.974,"10":104.004,"11":113.816,"12":123.845,"13":133.875,"14":144.123,"15":153.935,"16":163.964,"17":173.994,"18":184.024,"19":194.054,"20":203.865,"21":213.895},"1":{"0":12.821,"1":12.919,"2":3.649,"3":1.381,"4":1.381,"5":4.635,"6":24.951,"7":37.081,"8":29.586,"9":27.219,"10":15.878,"11":8.481,"12":18.639,"13":24.26,"14":12.032,"15":12.426,"16":15.582,"17":23.57,"18":16.765,"19":9.172,"20":5.03,"21":8.383}}'
isr_json = '{"0":{"0":1.746,"1":3.71,"2":5.675,"3":7.639,"4":9.821,"5":11.786,"6":13.968,"7":15.714,"8":17.897,"9":19.861,"10":21.825,"11":23.79,"12":25.754,"13":27.718,"14":29.683,"15":31.865,"16":33.829,"17":35.794,"18":37.758,"19":39.722,"20":41.687,"21":43.651,"22":45.833,"23":47.798,"24":49.98,"25":51.726,"26":53.909,"27":55.873,"28":57.837,"29":59.802,"30":61.766,"31":63.73,"32":65.913,"33":67.877,"34":69.841,"35":71.806,"36":73.77,"37":75.734,"38":77.698,"39":79.772,"40":81.627,"41":83.81,"42":85.774,"43":87.738,"44":89.921,"45":91.885,"46":93.849,"47":95.813,"48":97.778,"49":99.742,"50":101.706,"51":103.671,"52":105.853,"53":107.817,"54":109.673,"55":111.746,"56":113.71,"57":115.893,"58":117.639,"59":119.603,"60":121.786,"61":123.75,"62":125.714,"63":127.897,"64":129.861,"65":131.935,"66":133.899,"67":135.972,"68":137.718,"69":139.683,"70":141.647,"71":143.611,"72":145.903,"73":147.758,"74":149.94,"75":151.905,"76":153.651,"77":155.615,"78":157.798,"79":159.762,"80":161.726,"81":163.909,"82":165.873,"83":167.837,"84":169.802,"85":171.984,"86":173.948,"87":175.913,"88":177.877,"89":179.841,"90":181.806,"91":183.988,"92":185.952,"93":187.917,"94":189.772,"95":191.845,"96":193.81,"97":195.992,"98":197.956,"99":199.921,"100":201.885,"101":204.067,"102":206.032,"103":207.887,"104":209.96,"105":212.143,"106":214.107,"107":215.853,"108":218.036,"109":220.0},"1":{"0":1.68,"1":2.651,"2":2.533,"3":2.84,"4":2.959,"5":2.036,"6":1.491,"7":2.58,"8":0.757,"9":0.828,"10":0.379,"11":1.964,"12":1.254,"13":0.331,"14":1.112,"15":0.97,"16":1.302,"17":0.663,"18":1.68,"19":0.734,"20":0.237,"21":1.42,"22":2.746,"23":1.349,"24":0.355,"25":1.42,"26":1.893,"27":1.302,"28":2.331,"29":2.036,"30":2.864,"31":2.012,"32":11.811,"33":4.166,"34":2.438,"35":2.935,"36":5.87,"37":5.278,"38":5.515,"39":2.97,"40":2.13,"41":1.538,"42":1.325,"43":5.751,"44":5.728,"45":2.462,"46":3.243,"47":1.657,"48":3.598,"49":1.112,"50":1.491,"51":1.112,"52":2.651,"53":3.243,"54":0.734,"55":0.521,"56":0.781,"57":3.479,"58":0.805,"59":0.805,"60":5.065,"61":5.254,"62":1.325,"63":3.148,"64":1.728,"65":5.479,"66":1.882,"67":4.923,"68":1.586,"69":3.053,"70":0.592,"71":0.45,"72":2.391,"73":1.018,"74":1.207,"75":2.331,"76":2.982,"77":1.373,"78":3.574,"79":1.538,"80":2.225,"81":3.74,"82":1.799,"83":2.84,"84":1.633,"85":7.669,"86":2.036,"87":1.562,"88":1.728,"89":4.639,"90":1.041,"91":1.716,"92":2.769,"93":0.852,"94":3.574,"95":0.734,"96":0.568,"97":2.746,"98":0.663,"99":0.805,"100":1.444,"101":0.71,"102":2.225,"103":0.876,"104":2.201,"105":2.296,"106":1.988,"107":0.615,"108":2.367,"109":2.947}}'
data_ins = pandas.read_json(insulin_json)
data_isr = pandas.read_json(isr_json)


def main():

    new_length = 110
    new_x = np.linspace(data_ins.iloc[:,0].min(), data_ins.iloc[:,0].max(), new_length)
    new_y = sp.interpolate.interp1d(data_ins.iloc[:,0], data_ins.iloc[:,1], kind='cubic')(new_x)

    # function of a biexponential decay
    # https://swharden.com/blog/2020-09-24-python-exponential-fit/
    # https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_exponential_decay_2phase.htm
    # https://pharmacy.ufl.edu/files/2013/01/two-compartment-model.pdf
    # http://websites.umich.edu/~elements/07chap/html/07prof5.htm
    #
    # "This type of behavior is observed, for example, in the radioactive decay 
    # of a mixture of two nuclides with different half-lives, one short 
    # lived and the other relatively longer-lived."
    # y=ae^(-bt) +ce^(-dt)
    # https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470126714.app4
    #
    # A short introduction to pharmacokinetics
    # R. URSO, P. BLARDI, G. GIORGI
    # https://www.europeanreview.org/wp/wp-content/uploads/6.pdf
    # https://www.certara.com/knowledge-base/simplifying-deconvolution/
    def biExp(x, a, b, c, d):
        return (a * np.exp(x*b)) + (c * np.exp(x*d))



    ################################################################
    # 2.8 and 5.0 minutes and 0.28 percent from                    #
    # Ultradian Oscillations of Insulin Secretion in Humans        #
    # -------------------------------------------------------------#
    # "A deconvolution method based on a biexponential             #
    # disappearance rate of insulin, assuming half-lives           #
    # for insulin of 2.8 and 5 min with a                          #
    # fractional slow component of 28% (19)"                       #
    # ref 19 is https://pubmed.ncbi.nlm.nih.gov/11134098/          #
    # https://doi.org/10.2337/diabetes.51.2007.S258                #
    # -------------------------------------------------------------#
    # Direct measurement of pulsatile insulin secretion from       #
    # the portal vein in human subjects                            #
    # -------------------------------------------------------------#
    # "a biexponential insulin disappearance model in the          #
    # systemic circulation, consisting of earlier directly         #
    # estimated half-lives of 2.8 and 5.0 min and a                #
    # fractional slow component of 0.28 in healthy fasting humans" #
    # https://pubmed.ncbi.nlm.nih.gov/11134098/ (ref 19 above)     #
    # -------------------------------------------------------------#
    # In humans at least 75% of insulin secretion arises from      #
    # punctuated insulin secretory bursts                          #
    # -------------------------------------------------------------#
    # "insulin kinetics of 2.8 min (first half-life),              #
    # 5.0 min (second half-life), and                              #
    # a fractional slow component of 0.28"                         #
    # https://pubmed.ncbi.nlm.nih.gov/9374676/                     #
    # -------------------------------------------------------------#

    a = 1 - 0.28 # I assume from graphpad link that speaks of percent
    b = -np.log(2)/5.0 # 1/minutes
    c = 0.28 # fractional slow component of 0.28?
    d = -np.log(2)/2.8 # 1/minutes

    # See if the exponential function looks correct compared to 
    # Concepts in Clinical Pharmacokinetics
    # By Joseph T. DiPiro
    # Page 78 Figure 6-9
    fig, ax = plt.subplots(1, figsize=(14,10))
    ax.semilogy()
    ax.title.set_text("Both components in blue, and lines showing each constituent slope")
    ax.yaxis.set_label_text("Insulin (mU/ml) or ISR (mU/ml/min)")
    ax.xaxis.set_label_text("Time (minutes)")

    x_window = 100
    x_values = np.linspace(-x_window, x_window)
    # graph both
    ax.plot(x_values, biExp(x_values, a,b,c,d), 
        marker = '_', color = 'blue')

    # only graph the fast part (first)
    ax.plot(x_values, biExp(x_values, a,b,0,d), 
        marker = 'x', color = 'green')

    # only graph the slow part (second)
    ax.plot(x_values, biExp(x_values, 0,b,c,d), 
        marker = 'x', color = 'orange')



    # now calculate the actual convolution
    fig, ax = plt.subplots(3, figsize=(14,24))

    # create a file-like object from the url
    f = urllib.request.urlopen("https://diabetes.diabetesjournals.org/content/diabetes/51/suppl_1/S258/F2.large.jpg?width=800&height=600&carousel=1")

    # turn off axis since they really don't add much
    ax[0].set_axis_off()
    # read the image file in a numpy array
    img = plt.imread(f, format='jpg')
    ax[0].imshow(img, cmap='gray')


    # display ISR which was calculated orginally from the deconvolution
    # of measured Insulin Concentration
    ax[1].title.set_text("Insulin Secretion Rate")
    ax[1].yaxis.set_label_text("ISR (mU/ml/min)")
    ax[1].xaxis.set_label_text("Time (minutes)")
    ax[1].plot(data_isr.iloc[:,0], data_isr.iloc[:,1])


    # display non convolution version
    # (needs some smoothing)
    ax[2].title.set_text("Insulin Content")
    ax[2].yaxis.set_label_text("INSULIN (mU/ml)")
    ax[2].xaxis.set_label_text("Time (minutes)")
    ax[2].plot(new_x, new_y, color = "red")
    ax[2].plot(data_ins.iloc[:,0], data_ins.iloc[:,1], color = "blue")

    # convolution of ISR
    ax[2].plot(data_isr.iloc[:,0], np.convolve(data_isr.iloc[:,1],
        biExp(data_isr.iloc[:,1], a,b,c,d), mode = "same")/5, # bi exponential function
    color = "royalblue", marker = "d")

    ax[2].plot(new_x, np.convolve(new_y,
        biExp(new_y, a,b,c,d), mode = "same") / 10, # bi exponential function
    color = "maroon", marker = "d")

    # break them down to components
    ax[2].plot(new_x, np.convolve(new_y,
        biExp(new_y, 0,b,1,d), mode = "same") / 10, # bi exponential function
    color = "firebrick", marker = "d")


    ax[2].plot(new_x, np.convolve(new_y,
        biExp(new_y, 1,b,0,d), mode = "same") / 10, # bi exponential function
    color = "tomato",  marker = "d")


if __name__ == "__main__":
    main()

CSV 数据导出

Time (M),Insulin (mU/ml),Time (M),ISR (mU/ml/min)
4.143,12.821,1.746,1.680
13.954,12.919,3.710,2.651
23.984,3.649,5.675,2.533
34.014,1.381,7.639,2.840
44.044,1.381,9.821,2.959
54.073,4.635,11.786,2.036
64.103,24.951,13.968,1.491
74.133,37.081,15.714,2.580
83.944,29.586,17.897,0.757
93.974,27.219,19.861,0.828
104.004,15.878,21.825,0.379
113.816,8.481,23.790,1.964
123.845,18.639,25.754,1.254
133.875,24.260,27.718,0.331
144.123,12.032,29.683,1.112
153.935,12.426,31.865,0.970
163.964,15.582,33.829,1.302
173.994,23.570,35.794,0.663
184.024,16.765,37.758,1.680
194.054,9.172,39.722,0.734
203.865,5.030,41.687,0.237
213.895,8.383,43.651,1.420
,,45.833,2.746
,,47.798,1.349
,,49.980,0.355
,,51.726,1.420
,,53.909,1.893
,,55.873,1.302
,,57.837,2.331
,,59.802,2.036
,,61.766,2.864
,,63.730,2.012
,,65.913,11.811
,,67.877,4.166
,,69.841,2.438
,,71.806,2.935
,,73.770,5.870
,,75.734,5.278
,,77.698,5.515
,,79.772,2.970
,,81.627,2.130
,,83.810,1.538
,,85.774,1.325
,,87.738,5.751
,,89.921,5.728
,,91.885,2.462
,,93.849,3.243
,,95.813,1.657
,,97.778,3.598
,,99.742,1.112
,,101.706,1.491
,,103.671,1.112
,,105.853,2.651
,,107.817,3.243
,,109.673,0.734
,,111.746,0.521
,,113.710,0.781
,,115.893,3.479
,,117.639,0.805
,,119.603,0.805
,,121.786,5.065
,,123.750,5.254
,,125.714,1.325
,,127.897,3.148
,,129.861,1.728
,,131.935,5.479
,,133.899,1.882
,,135.972,4.923
,,137.718,1.586
,,139.683,3.053
,,141.647,0.592
,,143.611,0.450
,,145.903,2.391
,,147.758,1.018
,,149.940,1.207
,,151.905,2.331
,,153.651,2.982
,,155.615,1.373
,,157.798,3.574
,,159.762,1.538
,,161.726,2.225
,,163.909,3.740
,,165.873,1.799
,,167.837,2.840
,,169.802,1.633
,,171.984,7.669
,,173.948,2.036
,,175.913,1.562
,,177.877,1.728
,,179.841,4.639
,,181.806,1.041
,,183.988,1.716
,,185.952,2.769
,,187.917,0.852
,,189.772,3.574
,,191.845,0.734
,,193.810,0.568
,,195.992,2.746
,,197.956,0.663
,,199.921,0.805
,,201.885,1.444
,,204.067,0.710
,,206.032,2.225
,,207.887,0.876
,,209.960,2.201
,,212.143,2.296
,,214.107,1.988
,,215.853,0.615
,,218.036,2.367
,,220.000,2.947

2个回答

好吧，这里有几个问题需要解决。

ISR 是根据胰岛素浓度的反卷积计算的，然后在 ISR 上运行卷积应该给出胰岛素浓度，对吗？

也许。反卷积是一个定义不明确的问题。通常它只有在系统和数据都表现得相当好的情况下才有效，即使这样，信息和数据也通常会丢失很多。

我是否正确使用了双指数公式？我认为是这样，因为它与所绘制的两个斜率的书相匹配。

我不相信你会。指数衰减可以写成请注意，常数的单位是“时间上的一个”或频率。如果半衰期为 5 分钟，则相应的衰减曲线为半衰期在分母中，您需要的因子使其实际变为“一半”而不是

h (t) = e^{- a t}

$h(t) = e^{-at}$

h_{5} (t) = e^{- \log (2) \frac{t}{5 m i n}}

$h_5(t) = e^{-\log(2)\frac{t}{5min}}$

\log (2)

$\log (2)$

e^{- 1}

$e^{-1}$

您还有一个明显的缩放问题。IRS 在 10 左右达到峰值，而胰岛素上升到 40。单靠双指数模型无法解释这种增益差异。

一旦你修复了时间常数和缩放比例，事情看起来会好一些

模拟的胰岛素仍然更加摇摆不定，但这也是实际胰岛素被下采样 5 倍的结果。一般来说，这种类型的数据的采样率感觉非常低，并且可能存在明显的混叠也是。

这是 Matlab 代码。我通过将文本文件复制粘贴到电子表格中并将其保存为 .CSV 文件来导入数据

%% script
[a,b,c] = xlsread('foo.csv');  % read CSV file
%% eliminate the empty cells
t1 = a(:,1); % time
t1(isnan(t1)) = [];
x1 = a(:,2); % insulin
x1(isnan(x1)) = [];
t2 = a(:,3);  % time
x2 = a(:,4);  % IRS

%% buid the impules response and filter the data
h = .72*exp(-log(2)/5.*t2)+.28*exp(-log(2)/2.8.*t2);
% filter with the first 40 taps, that seems plenty
y = filter(h(1:40),1,x2);


%% plot it
clf;
plot(t2,[x2 y]);
hold on
% scale factor: match the RMS
x1s = x1./rms(x1)*rms(y);
plot(t1,x1s);
grid on
xlabel('Time in minutes');
legend('ISR','Simulated','Actual');
h = get(gca,'Children'); set(h,'Linewidth',2);

我确实想指出，我弄清楚了分数慢速组件是什么。根据参考资料

Ultradian Oscillations of Insulin Secretion in Humans
Chantal Simon and Gabrielle Brandenberger
Ref 19
-----------------------
Direct measurement of pulsatile insulin secretion from the portal vein in human subjects
S H Song  1 , S S McIntyre, H Shah, J D Veldhuis, P C Hayes, P C Butler
Ref 19
https://watermark.silverchair.com/jcem4491.pdf
Ref 6
-----------------------
1995: Pulsatile insulin secretion accounts for 70% of total insulin secretion during fasting
 NIELS PORKSEN, STEPHEN MUNN, JEFFERY STEERS, STEPHEN VORE, JOIIANNES VELDHUIS, AND PETER BUTLER
Ref 30
https://pubmed.ncbi.nlm.nih.gov/7573425/
-----------------------
1987: The pituitary gland secretes in bursts: Appraising the nature of glandularsecretoryimpulsesbysimultaneousmultiple-parameter deconvolution of plasma hormone concentrations
JOHANNES D. VELDHUIS*t, MARK L. CARLSONt, AND MICHAEL L. JOHNSON
https://www.pnas.org/content/84/21/7686
https://www.pnas.org/content/pnas/84/21/7686.full.pdf

在 1987 年的最后一篇论文中，它有“TheE(t-z) 项可以通过具有两个半衰期 HL1 和 HL2 的一个或两个分量模型和一个描述它们的相对关系的分数 f 来合理地描述贡献：”

E (t - z) = f e^{- (0.693 (t - z) / H L 1])} + (1 - f) e^{- 0.693 (t - z) / H L 2}

$E(t-z) = fe^{-(0.693(t-z)/HL1])} + (1-f)e^{-0.693(t-z)/HL2}$

这引用了参考文献 3，即 Lakowicz,JR(1983)Principles of Fluorescence Spectroscopy (Plenum, New York), pp. 65-87。

其它你可能感兴趣的问题

上一篇BPSK 与 QPSK 是否具有相同的性能？下一篇在频谱/频域中执行转置卷积？