简介：FID.zip是一份基于Bloch方程求解核磁共振自由感应衰减（FID）信号的Matlab代码包，面向NMR教学实验、脉冲序列设计以及弛豫机制分析等场景，适合物理、生物医学工程背景的学生与研究者使用。压缩包共6个文件，其中5个为m脚本、1个为txt说明文档，整体仅2KB，结构紧凑；脚本中既包含四阶Runge-Kutta法数值积分器（rk4sys），也有实现不同NMR脉冲激发与采集流程的多个子程序，txt数据文件则给出了三种不同初始条件或参数设置下的FID模拟结果，便于对照验证。已有202人学习下载。通过修改代码内T1、T2弛豫时间、拉莫尔频率及射频脉冲等参数，可系统探索FID信号的衰减形态与相位变化，理解T2弛豫对信号消失快慢的决定性作用。这套代码既适合课堂演示，也可作为科研中快速验证NMR模型的起点，能帮助使用者将抽象的Bloch方程转化为可视化的动态过程，提升对核磁共振物理图像与数值方法的掌握。
1. 为什么我用 80 行 Matlab 重写了一遍 Bloch 方程数值解
拿到FID.zip的时候，里面最值钱的不是那几个.m脚本，而是「FID 信号为什么长这样」这个问题本身。Bloch 方程把宏观磁化矢量 M 当作一个在旋转坐标系里被外场推着走的三维矢量，预cession、T1/T2 弛豫、射频脉冲全塞进一个四阶常微分方程组里，看似只有三个分量，实际写进代码后每个分量的耦合项都很容易被抄错一个符号，结果是曲线能画出来、衰减常数却对不上。这个包用rk4sys.m做四阶 Runge-Kutta 积分，SZCS.m到SZCS3.m对应三种输入情形，外加一份三种输入的情况(FID).txt用来核对数值结果。适合手上要做 NMR 弛豫建模、或者想搞懂 CPMG 这类序列为什么能压住 T2* 的人直接拿来当脚手架，比自己从零推公式快得多。
2. Bloch 方程的旋转坐标系形式与 RK4 积分骨架
2.1 从实验室系到旋转系的方程改写
教科书上的 Bloch 方程通常写成实验室系形式，B1 射频场以 Larmor 频率振荡，直接数值积分就得把步长压到纳秒级，不然一个周期都采不到几个点。实操里都会换到以 ω 旋转的坐标系，把射频场变成静态的 B1，方程简化为：
dMx/dt = -Ω·My - Mx/T2 dMy/dt = Ω·Mx - My/T2 + ω1·Mz dMz/dt = -ω1·My - (Mz-M0)/T1其中 Ω = ω0 - ω 是偏共振量，ω1 = γB1 是射频场强度，M0 是热平衡磁化强度。三种输入情形的差别就落在这三个参数的变化上：第一种是纯自由进动（ω1=0），第二种加一个恒定 B1，第三种带偏共振 Ω 不为零。写代码时最容易翻车的地方是 Ω 和 ω1 的单位——如果 γ 用 rad/(s·T)、场用 T，那时间轴就是秒；要是把 γ 按 MHz/T 的习惯填进去，T2 又给成毫秒，跑出来的衰减能差三个数量级。
2.2 rk4sys.m 的通用积分结构
包里的rk4sys.m是个通用四阶 RK4 求解器，输入是导数函数句柄、初值向量、时间向量和参数结构体，返回每个时间点的状态。它的核心是这个循环，我按自己的写法复现一下：
function Y = rk4sys(f, t, y0, p) % f: 导数函数句柄 @(t,y,p) % t: 时间向量，y0: 初值 [Mx My Mz]' % p: 参数结构体，含 T1 T2 omega1 Omega M0 n = numel(t); Y = zeros(n, numel(y0)); Y(1,:) = y0(:)'; for k = 1:n-1 h = t(k+1) - t(k); % 当前步长，允许非均匀 yk = Y(k,:)'; k1 = f(t(k), yk, p); k2 = f(t(k)+h/2, yk+h/2*k1, p); k3 = f(t(k)+h/2, yk+h/2*k2, p); k4 = f(t(k)+h, yk+h*k3, p); Y(k+1,:) = (yk + h/6*(k1 + 2*k2 + 2*k3 + k4))'; end end逐项拆开看：k1用当前步左端点的斜率，k2、k3取半步中点斜率，k4用右端点，最后按 1:2:2:1 加权求平均，截断误差是 O(h^5)，单步误差 O(h^6)。这意味着它比欧拉法能容忍大得多的步长。参数h按t的差分取，所以支持非均匀时间轴，处理「脉冲期间细采、自由进动期间粗采」很方便。导数函数 f 要写成独立文件或匿名函数，注意它的签名顺序必须是(t, y, p)，否则k2、k3的调用会直接报参数不足。
2.3 参数表与步长选择
把三种输入情形的参数落成表，照着改就能跑：
| 参数 | 含义 | 情形一（FID） | 情形二（加 B1） | 情形三（偏共振） |
|---|---|---|---|---|
| T1 | 纵向弛豫 | 1000 ms | 1000 ms | 1000 ms |
| T2 | 横向弛豫 | 50 ms | 50 ms | 50 ms |
| M0 | 平衡磁化 | 1 | 1 | 1 |
| omega1 | 射频强度 | 0 | 2π·1 kHz | 2π·1 kHz |
| Omega | 偏共振 | 0 | 0 | 2π·200 Hz |
| 时长 | 总采样 | 500 ms | 500 ms | 500 ms |
步长的经验值是 T2 的 1/50 到 1/100。T2=50 ms 时 h 取 0.5~1 ms 就够，再大曲线会明显偏。想验证精度，把 h 减半再跑一次，如果两条曲线肉眼重合，说明当前步长已经进入收敛区。
提示：RK4 是显式方法，遇到强射频 ω1 很大时步长必须压到 1/ω1 以下，否则会数值发散。这时候换隐式方法或直接用矩阵指数解更稳。
3. SZCS.m 与 SZCS1~3.m 三种输入情形的代码对照
3.1 导数函数怎么写
先给一份能直接跑的导数函数，把 Bloch 三个分量写全：
function dy = bloch_deriv(~, y, p) % y = [Mx; My; Mz]，p 含 T1 T2 omega1 Omega M0 Mx = y(1); My = y(2); Mz = y(3); dy = zeros(3,1); dy(1) = -p.Omega*My - Mx/p.T2; % 横向分量，含偏共振与T2衰减 dy(2) = p.Omega*Mx - My/p.T2 + p.omega1*Mz;% B1 把纵向磁化翻到横向 dy(3) = -p.omega1*My - (Mz - p.M0)/p.T1; % 纵向分量，T1 回平衡 end三个方程里，第一、二行的-M/T2是横向弛豫，第三行的(Mz-M0)/T1是纵向弛豫，omega1项负责把 Mz 和 My 互相耦合——这就是脉冲作用的数学来源。SZCS.m大概率就是这套结构的原版，SZCS1~3.m改的是p里的 omega1 和 Omega 初始值，以及t的构建方式。
3.2 三种情形的入口脚本差异
情形一最干净，ω1=0、Ω=0，初值取[0; 1; 0]（脉冲刚打完、磁化躺在横向面上），跑出来就是一条 exp(-t/T2) 包络的自由进动：
p = struct('T1',1000,'T2',50,'omega1',0,'Omega',0,'M0',1); t = 0:0.5:500; Y = rk4sys(@bloch_deriv, t, [0;1;0], p); signal = Y(:,1) + 1i*Y(:,2); % 复信号，实部虚部对应两路正交检波 plot(t, abs(signal)); xlabel('t / ms'); ylabel('|M_{xy}|');情形二把omega1换成2*pi*1（kHz 量级），你会看到 Mz 被反复翻转、|Mxy| 出现振荡，看起来像一段幅度调制的波形。情形三再把Omega设成2*pi*0.2，进动频率偏离参考频率，复信号的相位开始随时间线性增加，angle(signal)是一条斜率等于 Ω 的斜线——这条斜线就是化学位移的来源。
三种输入的情况(FID).txt的作用是给数值结果做交叉验证。常见做法是把 txt 里的参考值按列读进来，和自己的 Y 对齐相减，算一下最大绝对误差：
ref = load('三种输入的情况(FID).txt'); err = max(abs(ref(:,2) - abs(signal))); fprintf('最大偏差 = %.3e\n', err);如果 err 在 1e-3 量级以内，说明方程、单位、步长都没写错。要是偏大，优先查 T1/T2 是秒还是毫秒、omega1 是 Hz 还是 rad/s，这两个坑我见过太多次。
3.3 初值与单位的三类误用
第一类是把[0;1;0]写成[0;0;1]，那画出来的是纵向恢复曲线，形态是 (1-exp(-t/T1))，跟 FID 完全不沾边。第二类是 T2 填 50 但时间轴单位是秒，信号在 0.05 s 处就衰完，图上只剩一条平线。第三类是 omega1 用 Hz 直接代入，没乘 2π，导致进动周期算大 6.28 倍，包络形状对但尺度全错。这三种错误的共同点是曲线「看起来合理」，只有跟 txt 参考值对不上时才会暴露，所以每次改完参数都建议跑一遍误差比对。
4. 从 FID 信号反推 T2 与 CPMG 序列的数值验证
4.1 用数值解拟合衰减常数
FID 的核心信息是 T2，数值模拟后反向拟合一遍才算闭环。把包络取对数做线性回归，斜率就是 -1/T2：
env = abs(signal); idx = env > 0.05*max(env); % 丢掉噪声段，只拟合有效区 x = t(idx)'; y = log(env(idx)); coef = polyfit(x, y, 1); T2_fit = -1/coef(1); fprintf('拟合 T2 = %.2f ms, 设定值 = %.2f ms\n', T2_fit, p.T2);阈值0.05*max不是随便取的：FID 尾部接近数值噪声时取对数会被放得很大，把回归斜拉平。拟合值跟设定值差在 1% 以内，说明步长和单位都没问题。这个套路可以推广到任意序列——只要能拿到信号包络，就能反推等效弛豫时间。
4.2 CPMG 序列为什么能绕开 T2*
单次 FID 测到的是 T2*，它把 T2 和静磁场不均匀性混在一起，通常比真实 T2 短一大截。CPMG 用一串 180° 重聚焦脉冲把相位散开又掰回来，回波峰值只受 T2 支配。用数值法模拟这套序列，只需在时间轴上分段：每隔 TE/2 把 omega1 打开一小段、符号按 180° 脉冲翻转，其余时间 omega1=0。验证办法是改 Ω 看回波峰是否保持不变——如果峰高随 Ω 增大而稳，说明重聚焦逻辑写对了；如果峰高跟着 Ω 掉，多半是脉冲相位或时长设错了。
注意：180° 脉冲的 omega1×tp 必须等于 π。设 omega1=2π·1 kHz，tp 就得取 0.5 ms，少一点就是 167° 而非 180°，回波会出现残余相位。
4.3 常见异常与排查表
| 现象 | 可能原因 | 快速验证 |
|---|---|---|
| 曲线发散到 inf | 步长过大或显式 RK 遇强场 | h 减半重跑 |
| 衰减太快 | T2 单位按秒填成毫秒 | 查 p.T2 与 t 单位 |
| 相位不随时间走 | Omega 没进导数函数 | 打印 dy(1) 看有无 Omega 项 |
| 回波不出现 | 脉冲翻转角不是 180° | 检查 omega1*tp = π |
| 与 txt 偏差大 | omega1 少乘 2π | 对比角频率与频率 |
排查基本遵循「先单位、后公式、最后步长」的顺序，绝大多数问题在第一层就能定位。
本文还有配套的精品资源，点击获取
