CN104836547B - 一种短群延时数字滤波方法 - Google Patents
一种短群延时数字滤波方法 Download PDFInfo
- Publication number
- CN104836547B CN104836547B CN201510305222.9A CN201510305222A CN104836547B CN 104836547 B CN104836547 B CN 104836547B CN 201510305222 A CN201510305222 A CN 201510305222A CN 104836547 B CN104836547 B CN 104836547B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- cic
- comp
- digital
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
本发明公开了一种短群延时数字滤波方法。该方法将CIC(Cascaded Integrator‑Comb)数字滤波后的时域信号傅立叶变换到频域进行补偿，截取谱图中间部分后反傅立叶变换得到滤波后的时域数据。本发明采用CIC数字滤波器结合软件补偿，相比于传统的CIC数字滤波器结合FIR(Finite Impulse Response)数字滤波器补偿的方法，获得了更短的数字滤波器群延时。
Description
技术领域
本发明涉及一种数字滤波方法，具体涉及一种短群延时数字滤波方法，适用于核磁共振仪器以及其它需要进行数字滤波的信号处理系统。
背景技术
现代核磁共振仪器系统为了消除模拟滤波器所带来的失真，普遍采用了基于过采样的数字抽取滤波技术。现有的数字抽取滤波技术主要由CIC(Cascaded Integrator-Comb)数字滤波器与FIR(Finite Impulse Response)数字滤波器级联而成。CIC滤波器由于不需要乘法操作，具有硬件资源消耗小、执行效率高、群延时小的优点，但是其不平坦的通带幅频响应使得谱峰的相对幅度严重失真；而FIR滤波器具有可变的幅频特性和严格的线性相频特性，现代核磁共振仪器通常在CIC滤波器后级联高阶FIR作为补偿滤波器，得到了较好的通带平坦度、带外抑制度和过渡带陡峭度。然而，采用FIR作为补偿滤波器也存在一定的不足之处，主要是FIR滤波器带来较长的群延时。为了同时得到好的通带平坦度、带外抑制度和较窄的过渡带，FIR滤波器需要较大的阶数，阶数越大，群延时越大，信号向后推迟的程度就越大，核磁信号有效数据点的损失也就越大，从而一定程度上造成了信号的失真，这种失真在信号采集点数较少时则更加明显。
针对核磁共振仪器等系统中数字滤波器带来的群延时问题，科研工作者们提出了一些改进方法。Jon G.Wurl等(美国专利“DIGITAL FILTER PRE-CHARGING”US005652518)基于改变系统初始条件来优化滤波器瞬态响应的思想，提出核磁仪器接收部分在核磁信号前面根据第一个点的值增加一组平滑数据，经过数字滤波后，可实时的得到一组不含群延时的核磁共振FID(Free Induction Decay)信号。这种方法对核磁信号的第一个点的依赖性较大，而核磁信号的第一个点比较容易出现偏差，故该方法在核磁信号第一个点出现偏差时会产生谱图的基线扭曲；Petr SADOVSKY等从理论上对这种改变系统初始条件来优化滤波器瞬态响应的方法进行了分析与探讨(OPTIMISATION OF THE TRANSIENT RESPONSE OFA DIGITAL FILTER.Radioengineering 2000,9(2))，但第一个点出现偏差时谱图基线扭曲的问题并没有得到解决。
本专利为改善以上数字滤波器群延时带来的相关问题，提出一种短群延时数字滤波方法。具体内容为：先对过采样数据采用CIC实时数字滤波得到FID数字信号，然后在频域进行CIC的频域补偿并反傅立叶变换得到补偿后的FID信号。从而在不损失数字滤波器性能的前提下，满足了核磁共振实验的需求。
发明内容
本发明的目的是针对现有技术存在的问题，提出的一种短群延时数字滤波方法。
为了实现上述目的，本发明采用了以下技术措施：
一种短群延时数字滤波方法，包含以下步骤：
步骤1、核磁共振仪器载入设定的脉冲序列和采样参数，并按照设定的过采样频率对FID信号进行采样；
步骤2、将步骤1过采样得到的数字信号通过CIC数字滤波器进行CIC数字滤波；
步骤3、根据步骤2中CIC数字滤波器的幅频响应曲线Hcic(f)，计算CIC数字滤波器的频域幅度补偿曲线Hcomp(f)，将步骤2中CIC数字滤波后的数据进行傅立叶变换得到谱数据Scic，并将谱数据Scic与频域幅度补偿曲线Hcomp(f)逐点相乘，得到补偿后的谱数据Scomp；
步骤4、记步骤3中CIC补偿后的谱数据Scomp长度为L，截取Scomp居中的L/P个数据点，P>1，记为谱数据S1/P；
步骤5、对步骤4中得到的谱数据S1/P进行反傅立叶变换，得到经过CIC补偿后的FID信号。
如上所述的步骤3包括以下步骤：
步骤3.1、计算CIC数字滤波器的总幅频响应Hcic(f)：
其中，m表示级联CIC数字滤波器的个数，m>＝1，Mm，Rm，Nm，fm分别表示第m级CIC数字滤波器的延时因子、抽取率、阶数、频率与主瓣宽度的归一化值；
步骤3.2、根据步骤3.1获得的CIC数字滤波器的总幅频响应Hcic(f)计算频域幅度补偿曲线Hcomp(f)：
步骤3.3、将步骤2中CIC数字滤波得到的数据进行傅立叶变换得到谱数据Scic；
步骤3.4、将步骤3.3中的谱数据Scic与步骤3.2中的频域幅度补偿曲线Hcomp(f)逐点相乘，得到补偿后的谱数据Scomp。
本发明相对于现有技术，具有以下有益效果：
本发明采用CIC数字滤波器结合软件补偿，相比于传统的CIC数字滤波器结合FIR(Finite Impulse Response)数字滤波器补偿的方法，获得了更短的数字滤波器群延时。
附图说明
图1为本发明方法的总流程图。
图2为经过CIC数字滤波器后的FID信号。
图3为经过CIC数字滤波器后FID信号对应的谱图。
图4为CIC数字滤波器对应的频域幅度补偿曲线。
图5为CIC补偿完成后截取谱图中间1/4的谱图。
图6对截取后的谱图进行反傅立叶变换得到的FID信号。
图7-1和图7-2相同实验条件下CIC+FIR滤波器(图7-1)与本发明方法(图7-2)
的群延时对比图。
具体实施方式
下面通过实施例，结合附图，对本发明的技术方案作进一步的详细描述。
实施实例1：
一种短群延时数字滤波方法，总体流程图如图1所示，该方法总体包括以下步骤：
步骤1、核磁共振仪器载入设定的脉冲序列和采样参数，并按照设定的过采样频率对FID信号进行采样得到数字信号(在本实施例中：脉冲序列为单脉冲，采样参数主要包括谱宽sw＝10kHz，采样时间acqutime＝0.2048s，采样点数为2048，过采样频率Freqov＝60MHz)。
步骤2、将步骤1中采集得到的数字信号通过三级CIC数字滤波器进行CIC数字滤波，滤波后的FID数据如图2所示。(在本实施例中：共使用三级CIC滤波器，第一级滤波器的微分延时因子M1＝2，抽取率R1＝2，阶数N1＝3，第二级滤波器的微分延时因子M2＝2，抽取率R2＝3，阶数N2＝5，第三级滤波器的微分延时因子M3＝2，抽取率R3＝250，阶数N3＝5)
步骤3、根据三级CIC数字滤波器的幅频响应曲线Hcic(f)，计算三级CIC数字滤波器的频域幅度补偿曲线Hcomp(f)，将步骤2中三级CIC数字滤波后的FID数据进行傅立叶变换得到谱数据Scic(如图3所示)，并将谱数据Scic与CIC频域幅度补偿曲线Hcomp(f)逐点相乘，得到补偿后的谱数据Scomp。
(本实施例中：采用三级CIC滤波器，三级CIC滤波器的补偿曲线Hcomp(f)如图4所示，三级CIC滤波器的补偿曲线Hcomp(f)为：
其中
f1＝[-0.5/(R2*R3),0.5/(R2*R3),],f2＝[-0.5/R3,0.5/R3],f3＝[-0.5,0.5])
步骤4、记步骤3中CIC补偿后的谱数据Scomp长度为L，截取Scomp居中的L/P个数据点(P>1)，记为谱数据S1/P(在本实施例中：P＝4)，如图5所示。
步骤5、对步骤4中得到的谱数据S1/P进行反傅立叶变换，得到经过CIC补偿后的FID信号，如图6所示。
本方法与传统CIC+FIR方法的群延时对比：
如图7-1和图7-2所示，相同采样谱宽、采样时间、采样点数等实验条件下，传统的CIC结合256阶FIR数字滤波器方案得到的群延时为3.6ms(图7-1)，而采用本方法得到的群延时为0.2ms(图7-2)。
实施实例2：
步骤1、核磁共振仪器载入设定的脉冲序列和采样参数，并按照设定的过采样频率对FID信号进行采样得到数字信号(在本实施例中：脉冲序列为单脉冲，采样参数主要包括谱宽sw＝10kHz，采样时间acqutime＝0.2048s，采样点数为2048，过采样频率Freqov＝60MHz)。
步骤2、将步骤1中采集得到的数字信号通过三级CIC数字滤波器进行CIC数字滤波，滤波后的FID如图2所示。(在本实施例中：共使用三级CIC滤波器，第一级滤波器的微分延时因子M1＝2，抽取率R1＝2，阶数N1＝3，第二级滤波器的微分延时因子M2＝2，抽取率R2＝3，阶数N2＝5，第三级滤波器的微分延时因子M3＝2，抽取率R3＝250，阶数N3＝5)
步骤3、将步骤2中CIC数字滤波后的数据进行傅立叶变换，得到谱数据Scic(如图3所示)，记谱数据Scic长度为Lcic，截取谱数据Scic居中的Lcic/P个数据点(P>1)，记为谱数据S1/P(在本实施例中：P＝4)。
步骤4、根据步骤2中三级CIC数字滤波器的幅频响应曲线Hcic(f)，计算CIC数字滤波器的频域幅度补偿曲线Hcomp(f)，并将谱数据S1/P与CIC频域幅度补偿曲线Hcomp(f)逐点相乘，得到补偿后的谱数据Scomp。
(本实施例中：采用三级CIC滤波器，三级CIC滤波器的补偿曲线Hcomp(f)如图4所示，三级CIC滤波器的补偿曲线Hcomp(f)为：
其中
f1＝[-0.5/(4*R2*R3),0.5/(4*R2*R3),],f2＝[-0.5/(4*R3),0.5/(4*R3)],f3＝[-0.5/4,0.5/4])
步骤5、对步骤4中得到的谱数据Scomp进行反傅立叶变换，得到经过CIC补偿后的FID信号。
实施实例2是将实施实例1中步骤3(CIC补偿)和步骤4(谱数据截取)的先后顺序互换后的结果，两者得到补偿后的时域数据以及对应的群延时效果是一致的。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
Claims (2)
1.一种短群延时数字滤波方法，其特征在于，包含以下步骤：
步骤1、核磁共振仪器载入设定的脉冲序列和采样参数，并按照设定的过采样频率对FID信号进行采样；
步骤2、将步骤1过采样得到的数字信号通过CIC数字滤波器进行CIC数字滤波；
步骤3、根据步骤2中CIC数字滤波器的幅频响应曲线Hcic(f)，计算CIC数字滤波器的频域幅度补偿曲线Hcomp(f)，将步骤2中CIC数字滤波后的数据进行傅立叶变换得到谱数据Scic，并将谱数据Scic与频域幅度补偿曲线Hcomp(f)逐点相乘，得到补偿后的谱数据Scomp；
步骤4、记步骤3中CIC补偿后的谱数据Scomp长度为L，截取Scomp居中的L/P个数据点，P>1，记为谱数据S1/P；
步骤5、对步骤4中得到的谱数据S1/P进行反傅立叶变换，得到经过CIC补偿后的FID信号。
2.根据权利要求1所述的一种短群延时数字滤波方法，其特征在于，所述的步骤3包括以下步骤：
步骤3.1、计算CIC数字滤波器的总幅频响应Hcic(f)：
<mrow>
<msub>
<mi>H</mi>
<mi>cic</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>f</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mrow>
<mo>|</mo>
<mfrac>
<mrow>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mi>&pi;</mi>
<msub>
<mi>M</mi>
<mn>1</mn>
</msub>
<msub>
<mi>f</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mi>&pi;</mi>
<msub>
<mi>f</mi>
<mn>1</mn>
</msub>
<mo>/</mo>
<msub>
<mi>R</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>|</mo>
</mrow>
<msub>
<mi>N</mi>
<mn>1</mn>
</msub>
</msup>
<mo>*</mo>
<msup>
<mrow>
<mo>|</mo>
<mfrac>
<mrow>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mi>&pi;</mi>
<msub>
<mi>M</mi>
<mn>2</mn>
</msub>
<msub>
<mi>f</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mi>&pi;</mi>
<msub>
<mi>f</mi>
<mn>2</mn>
</msub>
<mo>/</mo>
<msub>
<mi>R</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>|</mo>
</mrow>
<msub>
<mi>N</mi>
<mn>2</mn>
</msub>
</msup>
<mo>*</mo>
<mo>.</mo>
<mo>.</mo>
<mo>.</mo>
<mo>.</mo>
<mo>.</mo>
<mo>.</mo>
<mo>*</mo>
<msup>
<mrow>
<mo>|</mo>
<mfrac>
<mrow>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mi>&pi;</mi>
<msub>
<mi>M</mi>
<mi>m</mi>
</msub>
<msub>
<mi>f</mi>
<mi>m</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>sin</mi>
<mrow>
<mo>(</mo>
<mi>&pi;</mi>
<msub>
<mi>f</mi>
<mi>m</mi>
</msub>
<mo>/</mo>
<msub>
<mi>R</mi>
<mi>m</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>|</mo>
</mrow>
<msub>
<mi>N</mi>
<mi>m</mi>
</msub>
</msup>
</mrow>
其中，m表示级联CIC数字滤波器的个数，m>＝1，Mm，Rm，Nm，fm分别表示第m级CIC数字滤波器的延时因子、抽取率、阶数、频率与主瓣宽度的归一化值；
步骤3.2、根据步骤3.1获得的CIC数字滤波器的总幅频响应Hcic(f)计算频域幅度补偿曲线Hcomp(f)：
<mrow>
<msub>
<mi>H</mi>
<mi>comp</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>f</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msub>
<mi>H</mi>
<mi>cic</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>f</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
步骤3.3、将步骤2中CIC数字滤波得到的数据进行傅立叶变换得到谱数据Scic；
步骤3.4、将步骤3.3中的谱数据Scic与步骤3.2中的频域幅度补偿曲线Hcomp(f)逐点相乘，得到补偿后的谱数据Scomp。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510305222.9A CN104836547B (zh) | 2015-06-05 | 2015-06-05 | 一种短群延时数字滤波方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510305222.9A CN104836547B (zh) | 2015-06-05 | 2015-06-05 | 一种短群延时数字滤波方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104836547A CN104836547A (zh) | 2015-08-12 |
| CN104836547B true CN104836547B (zh) | 2017-09-19 |
Family
ID=53814223
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510305222.9A Active CN104836547B (zh) | 2015-06-05 | 2015-06-05 | 一种短群延时数字滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104836547B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108011615B (zh) * | 2017-12-25 | 2020-06-02 | 北京怡和嘉业医疗科技股份有限公司 | 一种信号处理的方法和装置 |
| CN111966322A (zh) * | 2020-08-31 | 2020-11-20 | 广州视源电子科技股份有限公司 | 音频信号处理方法、装置、设备及存储介质 |
| CN113945878B (zh) * | 2021-10-13 | 2023-06-20 | 中国科学院精密测量科学与技术创新研究院 | 一种四核素同步磁共振成像和图像重建方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1588113A (zh) * | 2004-07-23 | 2005-03-02 | 华东师范大学 | 一种核磁共振成像信号的接收方法 |
| EP2667508A2 (en) * | 2012-05-21 | 2013-11-27 | STMicroelectronics Inc | Method and apparatus for efficient frequency-domain implementation of time-varying filters |
-
2015
- 2015-06-05 CN CN201510305222.9A patent/CN104836547B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1588113A (zh) * | 2004-07-23 | 2005-03-02 | 华东师范大学 | 一种核磁共振成像信号的接收方法 |
| EP2667508A2 (en) * | 2012-05-21 | 2013-11-27 | STMicroelectronics Inc | Method and apparatus for efficient frequency-domain implementation of time-varying filters |
Non-Patent Citations (1)
| Title |
|---|
| 数字化磁共振成像谱仪;徐勤;《中国优秀博硕士学位论文全文数据库》;20061015;第35-45、50-54页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104836547A (zh) | 2015-08-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104836547B (zh) | 一种短群延时数字滤波方法 | |
| CN103837884B (zh) | 基于时域分析的数字核脉冲信号梯形成形算法 | |
| CN103424183B (zh) | 机械振动信号检测异常干扰消除方法 | |
| CN104730580B (zh) | 地震资料异常振幅压制方法 | |
| CN102692650A (zh) | 一种具有假频压制功能的井筒波分离方法 | |
| CN105745706A (zh) | 用于扩展频带的装置、方法和程序 | |
| CN102055435A (zh) | 一种窄带数字滤波器 | |
| CN104614767A (zh) | 基于分段延拓的时变地震子波相位校正方法 | |
| CN105093282B (zh) | 基于频率约束的能量置换面波压制方法 | |
| CN104635264B (zh) | 叠前地震数据的处理方法及设备 | |
| CN104091591B (zh) | 一种音频处理方法及装置 | |
| CN104502702B (zh) | 检测电力信号的频率的方法和系统 | |
| CN105915193A (zh) | 一种用于多相滤波器的改进生成方法 | |
| CN108872402B (zh) | 超声波巴特沃斯、汉宁窗组合带阻滤波方法 | |
| CN103777221B (zh) | 基于窗函数法的数字核脉冲信号高斯成形方法 | |
| CN104579239B (zh) | 一种滤波系统的过滤方法 | |
| CN103941280B (zh) | 基于冲激响应不变法的数字核脉冲高斯成形方法 | |
| CN105811921B (zh) | 一种抑制工频谐波干扰的方法及滤波器 | |
| CN106356069B (zh) | 一种信号处理方法和装置 | |
| CN104392093A (zh) | 一种基于分数阶复合积分算子的零相位滤波器及其滤波方法 | |
| CN114894478A (zh) | 一种滚动轴承微弱故障特征提取方法 | |
| Li et al. | Compensation method for the CIC filter in digital down converter | |
| CN102571031B (zh) | 一种用于实现小波变换的宽带低插损声表面波滤波器组 | |
| CN112596097A (zh) | 一种基于权冲击函数的核信号前端处理系统 | |
| Wei et al. | De-noising the speech signal with fir filter based on matlab |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| EXSB | Decision made by sipo to initiate substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20190821 Address after: 430071 Wuchang District, Wuhan City, Hubei Province, Xiaohong Shanxi 30, Wuhan Institute of Biology, Chinese Academy of Sciences, No. 002 Patentee after: Wuhan Zhongke Yunchu Technology Co., Ltd. Address before: 430071 Wuchang, Wuhan, Hubei small Hong Kong District, No., No. 30 Patentee before: Wuhan Inst. of Physics and Mathematics, Chinese Academy of Sciences |