1

我编写了一些代码，使我能够从 DICOM 文件中获取数据并将数据分离为单独的 FID 信号。

from pylab import *
import dicom
import numpy as np
import scipy as sp

plan=dicom.read_file("1.3.46.670589.11.38085.5.22.3.1.4792.2013050818105496124")
all_points = array(plan.SpectroscopyData)
cmplx_data = all_points[0::2] -1j*all_points[1::2]
frames = int(plan.NumberOfFrames)
fid_pts = len(cmplx_data)/frames
del_t = plan.AcquisitionDuration / (frames * fid_pts)

fid_list = []
for fidN in arange(frames):
    offset = fidN * fid_pts 
    current_fid = cmplx_data[offset:offset+fid_pts]
    fid_list.append(current_fid)

我现在想量化其中的一些数据，因此，在应用傅立叶变换和移位后，我尝试使用 Scipy 的 quad 函数并遇到以下错误：

spec = fftshift(fft(fid_list[0])))
sp.integrate.quad(spec, 660.0, 700.0)



error                                     Traceback (most recent call last)
/home/dominicc/Experiments/In Vitro/Glu, Cr/Phantom 1: 10mM Cr, 5mM Glu/WIP_SV_PRESS_ME_128TEs_5ms_spacing_1828/<ipython-input-107-17cb50e45927> in <module>()
----> 1 sp.integrate.quad(fid, 660.0, 700.0)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
243     if type(args) != type(()): args = (args,)
244     if (weight is None):
--> 245         retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
246     else:
247         retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
307     if points is None:
308         if infbounds == 0:
--> 309             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
310         else:
311             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

error: First argument must be a callable function.

有人可以建议一种使该对象可调用的方法吗？在阅读了 Python 中的“可调用”是什么？. 读完之后，我尝试了

def fid():
    spec = fftshift(fft(fid_list[0])) 
    return spec

哪个返回了相同的错误。

任何帮助将不胜感激，谢谢。

4

1 回答 1

4

由于您正在积分仅在网格上定义的函数（即，数字数组），您需要使用“积分，给定固定样本”中的例程之一。

事实上，第一个参数quad() 必须是一个函数，你可以调用的东西（“可调用”）。例如，您可以这样做：

>>> from scipy.integrate import quad
>>> def f(x):
...     return x**2
... 
>>> quad(f, 0, 1)
(0.33333333333333337, 3.700743415417189e-15)

这可以通过采样来完成f()：

>>> from scipy.integrate import simps
>>> x_grid = numpy.linspace(0, 1, 101)  # 101 numbers between 0 and 1
>>> simps(x_grid**2, dx=x_grid[1]-x_grid[0])  # x_grid**2 is an *array*.  dx=0.01 between x_grid values
0.33333333333333337

您的情况就像在第二个示例中一样：您集成了一个采样函数，因此使用一个 or cumtrapz()，simps()or是很自然的romb()。

于 2013-05-09T13:50:50.717 回答