2016, Vol. 33文章信息
- 郭新宇, 肖亮
- GUO Xin-yu, XIAO Liang
- 基于数据采集卡和GPU的MRI信号接收设计
- Design of A Receiver for Magnetic Resonance Imaging Based on Data Acquisition Card and GPU
- 波谱学杂志, 2016, 33(4): 581-589
- Chinese Journal of Magnetic Resonance, 2016, 33(4): 581-589
- http://dx.doi.org/10.11938/cjmr20160407
-
文章历史
收稿日期: 2016-01-25
收修改稿日期: 2016-10-24
DOI:10.11938/cjmr20160407
谱仪[1, 2]是磁共振成像（Magnetic Resonance Imaging，MRI）系统的核心设备，负责行使成像序列的运行、射频脉冲的发生、梯度波形的发生和回波信号的接收等功能. 磁共振图像质量与信号采集和处理过程密切相关.因此设计高性能的信号接收方法具有重要意义.
传统的模拟正交解调，存在两个正交通道不匹配的问题，易导致图像出现伪影[3].
为了解决这一问题，数字化的接收机成为了主要的发展方向[4].随着高速模数转换器（Analog Digital Converter，ADC）技术的发展，在低场MRI系统中，出现了直接采样数字接收系统[5-7].它对经过前置放大器放大的磁共振信号直接进行模数转换，使用专用数字下变频器（Digital Down Converter，DDC）实现数字混频、数字滤波和抽取.这种磁共振信号接收方式具有较好的一致性，并且简化了设计.但DDC的信号处理过程不够灵活，有限字长效应还会造成信号精度损失，最终导致重建图像信噪比降低[8].
本文设计了一种基于高速数据采集卡和图形处理器（Graphics Processing Unit，GPU）的数字接收机，该接收机不再使用专用DDC，而是使用计算机上的高速数据采集卡对回波信号直接采样（采样速率为40 MHz），再使用GPU并行编程实现数字混频、多级数字滤波抽取和图像重建.GPU通常含有数百至上千个计算单元，浮点运算性能可以达到中央处理器（Central Processing Unit, CPU）的数十倍，因此本文设计的接收方法可以显著提高信号处理的速度.由于整个接收算法完全采用软件编程实现，因此有利于设计多种数字滤波抽取方案，提高了整个接收机的灵活性.同时，信号处理过程全部采用双精度浮点格式运算，可以进一步改善信号运算精度，提高图像信噪比.
1 数字接收机设计本文设计的数字接收机主要由高速数据采集卡、GPU和PC机组成.高速数据采集卡和GPU分别通过外部设备互联（Peripheral Component Interconnect, PCI）总线和外部设备高速互联（PCI Express, PCI-E）总线与内存进行数据交互.高速数据采集卡采用ADLINK公司PCI-9846[9]，其具有4路采集通道和16-bit ADC，采样速率可达40 MHz.
GPU采用NVIDIA Geforce GTX TITAN BLACK，其具有2 880个统一计算设备架构（Compute Unified Device Architecture, CUDA）处理器核心（CUDA cores），全局存储器容量为6 144 MB.CPU为Intel Core i7-4770，内存为16 GB.
数据采集卡在谱仪输出的触发脉冲的控制下对前置放大器输出的信号进行采样，数据直接传输至PC机内存进行存储.数字接收程序运行时，将数据采集卡采样的数据通过PCI-E总线传输至GPU的全局存储器中，并行实现数字混频、多级数字滤波抽取和图像重建.图像重建使用cuFFT[10]中的库函数完成二维快速傅里叶变换（Fast Fourier Transform，FFT）.cuFFT是CUDA[11]提供的经过优化的GPU并行FFT函数库，可以快速完成图像重建.最后，通过PCI-E总线将重建的磁共振图像数据传输到PC机内存进行显示.整个接收系统的结构如图 1所示.
同时，为了保证射频发射和接收相位相干，本设计采用了以下三点思路：第一，数据采集卡的采样时钟与射频合成的时钟来源于同一时钟，频率为40 MHz.第二，控制数据采集卡采样的触发脉冲由谱仪产生，其发生时刻在脉冲序列中的相对位置是固定不变的.第三，信号接收时，数字混频用到的数字本振由GPU程序预先计算，其设置的相位是精确的.
|
| 图 1 接收机设计框图 Fig. 1 Block diagram of the proposed digital receiver |
数字混频并不改变同相和正交（I和Q）两路数字信号的速率，很难直接对其进行处理，因此需要对该数字信号进行抽取.以采样速率为40 MHz、成像带宽为25 kHz为例，抽取的倍数为40 MHz/25 kHz = 1 600.将1 600分解为2×2×2×2×2×2×5×5，可以使用这些因子的组合设计抽取的倍数.根据信号处理理论，在抽取之前需要进行滤波以防止频谱混叠.考虑到采用一级滤波无法设计出通带和阻带性能都满足要求的滤波器，所以采用三级滤波抽取算法.
对于使用专用数字下变频器进行信号接收的系统，其通常采用级联积分梳状滤波器（Cascade Integrate Comb，CIC）和有限长冲激响应滤波器（Finite Impulse Response，FIR）的组合.以文献[12]为例，其采用DDC AD6620[13]作为数字接收系统的核心器件实现数字混频和数字滤波抽取.AD6620包含二阶级联积分梳状滤波器（CIC2），五阶级联积分梳状滤波器（CIC5）和可编程FIR滤波器.作为对比，这里首先将CIC2、CIC5和FIR滤波器组成的抽取滤波器组应用于本设计的信号接收中.实验表明该方法由于使用双精度浮点数据类型，具有很高的精度，不存在专用DDC中每一级数字滤波后都要对数据进行动态范围压缩的问题，相比于使用AD6620的接收系统所成的图像，信噪比有所提高.
专用DDC片内物理资源有限，例如AD6620，只能采用CIC滤波器和FIR滤波器的组合.但是由于CIC滤波器单级衰减只有13 dB，通常需要多级级联以增加阻带衰减.然而与此同时，通带纹波也会随之增加.对于N级CIC滤波器，通带纹波也增加N倍.因此，本文设计了完全由FIR滤波器组成的抽取滤波器组.每一级FIR滤波器均为128阶，前两级FIR滤波器使用Hamming窗设计，第三级滤波器采用等波纹滤波器.前文中提到以采样速率为40 MHz、成像带宽为25 kHz为例，抽取倍数应为1 600，综合考虑每一级滤波器的性能，本设计中第1级至第3级FIR滤波器的抽取系数分别为16、25和4.以该抽取系数设计的抽取滤波器的特性如图 2所示.第1级FIR滤波器的截止频率为1.25 MHz，通带纹波为0.04 dB；第2级FIR滤波器的截止频率为0.05 MHz，通带纹波为0.009 dB；第3级FIR滤波器的截止频率为12.5 kHz，通带纹波为0.02 dB.
|
| 图 2 抽取滤波器的幅频特性和相频特性. (a) 第1级FIR滤波器的幅频特性；(b) 第1级FIR滤波器的相频特性；(c) 第2级FIR滤波器的幅频特性；(d) 第2级FIR滤波器的相频特性；(e) 第3级FIR滤波器的幅频特性；(f) 第3级FIR滤波器的相频特性 Fig. 2 Magnitude frequency response and phase frequency response of the decimating filters. (a) Magnitude frequency response of the first stage FIR filter; (b) Phase frequency response of the first stage FIR filter; (c) Magnitude frequency response of the second stage FIR filter; (d) Phase frequency response of the second stage FIR filter; (e) Magnitude frequency response of the third stage FIR filter; (f) Phase frequency response of the third stage FIR filter |
CUDA是NVIDIA公司推出的通用平行计算平台和编程模型，使GPU可以方便的用于通用计算.CUDA架构下，代码分为在CPU上运行的主机代码和在GPU上运行的设备代码（内核函数）.主机代码与设备代码通过PCI-E总线进行数据交互.在GPU上执行的内核函数需要用户为其分配合理的线程块和线程.每个内核函数都在一个网格中执行，网格由线程块组成，而每个线程块包含了多个线程，通过CUDA API为内核函数指定网格和线程块的维度.
经过数据采集卡采样的磁共振信号以双精度浮点格式存储.又以采样速率为 40 MHz、成像带宽为25 kHz为例，平均次数（NS）为2，256个相位编码步.为重建一帧256×256像素的图像，每行需采样256次，每次采集时间为40 μs，每行总计采样409 600个点，共需要采样512行.考虑到滤波器的死时间，每行采样500 000个点.将该数据存储在大小为512×500 000的矩阵中，使用GPU编程并行实现对每行数据的数字混频、多级滤波抽取，最后完成图像重建.在GPU的全局存储器中开辟大小为512×500 000的空间，用于存储数字混频后的I、Q数据.进行数字滤波抽取时，内核函数线程块中的每一个线程对应于存储数字混频后的数据矩阵的一个元素，负责计算每个元素滤波后的值，并依据线程索引完成抽取.所有线程计算完成后，调用cuFFT中的二维FFT库函数进行图像重建，最后将重建的结果从GPU全局存储器中传输至PC机内存.数字混频、多级滤波抽取和图像重建的流程如图 3所示.
|
| 图 3 数字混频、多级滤波抽取和图像重建的流程 Fig. 3 Process of digital mixing, multi-stage filter and decimation and image reconstruction |
表 1为GPU与C语言和MATLAB的计算速度对比，给出了实现全部数字混频，多级滤波抽取和图像重建的时间.可见，GPU的计算速度相比MATLAB能提高大约10~30倍，相比C语言能提高18~43倍.
| t (GPU)/s | t (MATLAB)/s | t (C语言)/s | |
| CIC2+CIC5+FIR | 1.48 | 14.31 | 26.71 |
| FIR+FIR+FIR | 1.96 | 64.65 | 85.76 |
首先，给出对高速数据采集卡采样的512行数据中的一行进行信号处理的结果，如图 4所示.图 4(a)为该行的采样数据.图 4(b)、(c)和(d)分别为该行数据经过每一级FIR滤波抽取后的时域图形.图 4(e)为滤波抽取后待重建信号的频谱.
|
| 图 4 高速数据采集卡采样的512行数据中的一行进行信号处理的结果. (a)中间一行采样数据；(b)第1级滤波后的波形；(c)第2级滤波后的波形；(d)第3级滤波后的波形；(e)中间一行信号处理后的频谱 Fig. 4 Signal processing results of a line. (a) Sampling data of a line; (b) Output signal of the first stage filter; (c) Output signal of the second stage filter; (d) Output signal of the third stage filter; (e) Output spectrum of a line |
其次，为进一步验证本文设计的数字接收机的性能，在0.3 T永磁系统中进行了成像试验.图 5为水模的自旋回波T1加权图像.图 6为人脑的自旋回波T1加权图像. 实验参数如下：脉冲间隔时间（TR）= 450 ms，回波时间（TE）= 16 ms，视野（FOV）= 240 mm×240 mm，扫描次数（NS）为2，中心频率为12.78 MHz.
图 5(a)、(b)为本文设计的数字接收机所成的图像，其信噪比（Signal to Noise Ratio, SNR）分别为67.93、67.75.图 5(c)为使用DDC AD6620所成的图像，信噪比为57.04. 相比于图 5(c)，图 5(a)、(b)的信噪比分别提高了19.1%和18.8%.从图 6也可以看到，使用本文设计的数字接收机得到的人脑组织成像结构较清楚，没有干扰线与直流亮点.
|
| 图 5 水模的自旋回波T1加权图像. (a)使用CIC2+CIC5+FIR滤波器的成像结果；(b)使用三级FIR滤波器的成像结果；(c)使用DDC AD6620的成像结果 Fig. 5 T1-weighted spin echo images of phantom using (a) CIC2+CIC5+FIR filters, (b) three-stage FIR filters and (c) DDC AD6620 |
|
| 图 6 人脑的自旋回波T1加权图像. (a)使用CIC2+CIC5+FIR滤波器的成像结果；(b)使用三级FIR滤波器的成像结果 Fig. 6 T1-weighted spin echo images of brain using (a) CIC2+CIC5+FIR filters and (b) three-stage FIR filters |
本文提出了基于高速数据采集卡和GPU的磁共振信号接收方法.数据采集卡在谱仪输出的触发脉冲控制下直接采样前置放大器输出的信号，由GPU并行编程实现数字混频、多级数字滤波抽取和图像重建.由于采用商业化的高速数据采集卡和GPU，简化了系统设计，采用软件编程实现信号接收结合双精度浮点数据格式，提高了运算速度和重建图像的信噪比，增强了系统的灵活性.MRI结果验证了本设计的可行性.
| [1] | Liu Min(刘敏), Qiu Wen-qi(邱雯绮), Sun Hui-jun(孙惠军) et al . Progress in the portable NMR spectrometer(便携式核磁共振谱仪的研究进展)[J]. Chinese J Magn Reson(波谱学杂志) , 2014, 31 (4) : 504-514 |
| [2] | Li G Y, Jiang Y, Yan X L . Digital nuclear magnetic resonance spectrometer[J]. Rev Sci Instrum , 1994, 72 : 4460-4463 |
| [3] | Kasal M, Halamek J, Husek V et al . Signal processing in transceiver for nuclear magnetic resonance and imaging[J]. Rev Sci Instrum , 1994, 65 : 1897-1994 DOI:10.1063/1.1144840 |
| [4] | Liu Chao-yang(刘朝阳), Tan Ping(谭萍), Qiu Jian-qing(裘鉴卿) et al . The principle & implement of digital receiver for NMR(核磁共振谱仪数字接收机的原理与实现)[J]. Chinese J Magn Reson(波谱学杂志) , 2000, 17 (2) : 161-165 |
| [5] | Xiao L, Sun B H, Ma Y P et al . A digital magnetic resonance imaging spectrometer using digital signal processor and field programmable gate array[J]. Rev Sci Instrum , 2013, 84 (5) : 054702 DOI:10.1063/1.4803007 |
| [6] | Ren Jie-jing(任洁静), Xu Qin(徐勤), Li Geng-ying(李鲠颖) . A digital receiver based on direct RF sampling for low-field MRI scanner(用于低场磁共振成像的直接射频采样数字接收机)[J]. Chinese J Magn Reson(波谱学杂志) , 2007, 24 (1) : 27-33 |
| [7] | Hu Peng(胡鹏), Yang Liang(杨亮), Hu Shao-bin(胡少斌) et al . Design and implementation of NMR receiver(核磁共振波谱仪接收机的设计与实现)[J]. Chinese J Magn Reson(波谱学杂志) , 2011, 28 (1) : 42-49 |
| [8] | Sun Hong-yu(孙宏宇), Tang Wei-nan(汤伟男), Wang Wei-min(王为民) . Improving signal-to-noise ratio in magnetic resonance imaging using dynamic receiver gain(利用动态接收增益提高磁共振成像信噪比的方法)[J]. Chinese J Magn Reson(波谱学杂志) , 2014, 31 (4) : 515-522 |
| [9] | PCI-9846 User's Manual, ADLINK Technology Inc[EB]. 2010. |
| [10] | CUFFT Library User's Guide, NVIDIA Corporation[EB]. 2014. |
| [11] | CUDA C Programming Guide, NVIDIA Corporation[EB]. 2014. |
| [12] | Xiao L, Xin T, Tang W N et al . A high-field magnetic resonance imaging spectrometer using an oven-controlled crystal oscillator as the local oscillator of its radio frequency transceiver[J]. Rev Sci Instrum , 2014, 85 (9) : 094705 DOI:10.1063/1.4894657 |
| [13] | AD6620 Data Sheet, Analog Devices Inc[EB]. 2001. |
本作品采用知识共享署名 4.0 国际许可协议进行许可。