简介：针对声子晶体声传播特性分析的一维传递矩阵法MATLAB程序包，面向声学、凝聚态物理及结构振动领域的研究人员与学习者，用于快速模拟周期性复合结构的声波传播特性。压缩包仅包含单个m文件，体积约7KB，结构精简，便于直接调用与二次修改。程序基于一维传递矩阵法，将晶体逐段分解构建总体传递矩阵，并结合边界条件求解不同频率下的声波传递率，据描述其计算正确率可达98%，可作为科研与工程设计的辅助工具。已有246人学习下载，使用者可通过修改单元长度、材料密度、声速和周期数等参数，获得频率-传递率曲线，直观分析禁带位置与传输衰减特性，为声屏障、噪声控制及新型声学器件设计提供数据支撑。程序实现包含结构定义、传递矩阵串联、边界条件设置与结果可视化等关键环节，适合希望掌握计算原理并快速上手的入门与中级用户。

1. 拿到 huisang_v13.zip 先别急着跑：传递矩阵法算声子晶体到底在算什么

如果你像我一样收到过huisang_v13.zip这种命名接近随机数的压缩包，大概率是同行或学生丢来的声子晶体计算程序。这里面的传递矩阵法，是声子晶体能带和传递率计算里最常用的一类半解析方法：把每一层介质的状态向量沿厚度方向乘起来，再用 Bloch 周期条件提取色散关系；要算有限周期结构的透射率，就把单胞矩阵连乘 N 次。它能解决什么？不需要建模、不需要划网格，几毫秒扫完一条 0–20 kHz 的传递率曲线，特别适合参数扫描找禁带。适合谁：做超材料隔声、周期杆/周期板减振、声学超表面的人。第一次打开这个压缩包，先别急着点 run，花十分钟搞清楚它用的是物理状态基还是幅值基——这两个基的矩阵长得像、符号习惯差得远，后面很多翻车都是从这里开始的。

2. 传递矩阵法核心表达式：从位移连续到特征矩阵的推导与参数选择

TMM 在一维分层结构里的地位相当于查表算法：频率定下来，材料定下来，每一层就是一个 2×2 矩阵，从头乘到尾，边界条件就自动满足了。先别急着往压缩包里找代码，把这两条路径在脑子里过一遍，后面改参数才不会瞎。

2.1 为什么传递矩阵法能在一维声子晶体里横着走

一维声子晶体最常见的模型是周期叠层：A 层和 B 层交替堆叠，纵向波垂直入射。对这种结构，每一层内部的波动解就是两个方向的行波叠加，不需要像有限元那样把每层再切成一大堆单元。常见做法是用声阻抗 Z=ρc 和相位 φ=ωd/c 概括整层，层与层之间只要满足位移连续和应力连续两个条件就够了。

有限元能算任意 2D/3D 结构，但做参数扫描时成本感人：一个频率一个模型，扫 2000 个频点就得解 2000 次非稀疏矩阵。而传递矩阵法扫完一整条传递率曲线只需要一个循环，每个频点只做几个 2×2 矩阵乘起来再加一次特征值分解，这个成本差异就是它在一维周期结构里“横着走”的原因。局限也很明确：它只能吃分层结构和平面波垂直入射这盘菜；要算二维周期格子或弹性波斜入射，这包就只能当参考了。

2.2 层传递矩阵的两种写法：物理状态基与幅值基

先说物理状态基。取位移 u 和正应力 σ 组成状态向量，单层介质内部解是 u(x)=Ae^{-jkx}+Be^{jkx}，应力和位移通过弹性模量 E 关联。定义声阻抗 Z=ρc、相位 φ=kd，单层从左到右的传递矩阵是：

[[cos φ, sin φ / Z], [-Z sin φ, cos φ]]

这个矩阵行列式为 1，两个特征值是 e^{±jφ}，分别对应右行波和左行波穿过这一层的相位积累。注意左下角的负号：时间因子取 e^{+jωt} 还是 e^{-iωt}，负号位置会互换。你拿到的老代码十有八九最多只按其中一种约定写成，直接照抄很容易把透射率的相位弄反。

再说幅值基。每层内只保留前进波和后退波振幅 A、B。穿过一个厚度为 d 的均匀层，振幅变化只有相位项，传播矩阵写成对角阵：

P = diag(e^{-jφ}, e^{jφ})

穿过一个阻抗从 Z1 变到 Z2 的界面，由位移连续和应力连续可以解出界面矩阵：

I(1→2) = [[(1+r)/2, (1-r)/2], [(1-r)/2, (1+r)/2]]，其中r=Z1/Z2

当两测阻抗相等时，I 退化成单位阵，物理意义正好是“没有反射没有透射损失”。用幅值基算透射率，最后只需要取全局矩阵 M 的 M22 的倒数，这是最干净的一条路。

2.3 Bloch 条件与能带映射：trace(Q) 的临界值就是禁带

无限周期结构的核心是 Floquet 定理：在周期边界处，状态向量只差一个复数倍 λ=e^{jqa}，a 是单胞厚度，q 是 Bloch 波数。把幅值基下的单胞矩阵 Q 定义为：

Q = I(B→A) @ P(B) @ I(A→B) @ P(A)

这里的坐标取在 A 层起点，也就是每个胞刚进入 A 层的位置。Q 把当前胞 A 层起点处的振幅向量映射到下一个胞 A 层起点。因为 Q 是 2×2 且行列式为 1，特征值是 λ 和 1/λ，于是：

trace(Q) = λ + 1/λ = 2cos(qa)

当|trace(Q)| > 2时，q 变成复数，Bloch 波在传播方向上指数衰减，这就是禁带。画色散图时只要把每个频率的 trace(Q)/2 画在纵轴上，越过 ±1 的区间就是禁带；画传递率曲线时，同一频率范围对应透射低谷。这一条能直接验证你手上压缩包里的能带图和传递率曲线对不对。

3. 用 huisang_v13.zip 里的核心函数跑通单胞传递矩阵：最小复现流程与代码解读

理论立住之后，动手第一步是认文件，第二步是找一个能跑的最小例子，第三步是验证输出。别一上来就换成自己的材料参数，先用钢/环氧这种经典组合跑一遍，出了问题容易定位。

3.1 先认文件与输入：v13 常见目录结构和主函数入口

压缩包解压进去，常见做法是先看有没有 README 或注释头。我收到的这类包通常分成三类文件：主扫频脚本、单胞矩阵构造函数、画图脚本。文件名不一定叫transfer_matrix还是tmm_main，但输入基本逃不出这五个量：

如果压缩包里没有说明文档，先找主函数入口里最外层那个 for 循环：循环里每轮算一个频率点，内部调用的就是单胞矩阵函数。用 v13 这类编号很高的包还有个额外好处，量大说明迭代多，核心公式大概率是被验证过的；坏处是变量名可能改得和论文对不上，所以不要只看函数名，要看矩阵乘法顺序。

3.2 用钢/环氧双层声子晶体跑通传递矩阵：核心代码与参数表

下面我用 Python 和 numpy 复现最核心的幅值基流程。宿主介质取钢，单胞由钢层和环氧层组成，N=8。代码里保留注释，方便对照前面的公式。

import numpy as np # ==== 材料与几何输入 ==== rho = np.array([7800.0, 1180.0]) # 钢, 环氧树脂 (kg/m3) c = np.array([5200.0, 2200.0]) # 纵波声速 (m/s) d = np.array([0.02, 0.02]) # 单层厚度 (m) N = 8 # 周期数 Z = rho * c # 各层特征阻抗 Z0 = Z[0] # 宿主介质用钢的阻抗 a = d[0] + d[1] # 单胞总厚度 # ==== 幅值基基础矩阵 ==== def interface(Z_left, Z_right): r = Z_left / Z_right return np.array([[(1 + r) / 2, (1 - r) / 2], [(1 - r) / 2, (1 + r) / 2]], dtype=complex) def propagation(phi): return np.diag([np.exp(-1j * phi), np.exp(1j * phi)]) def cell_transfer(omega): I_host_A = interface(Z0, Z[0]) # 宿主 -> 钢 I_A_B = interface(Z[0], Z[1]) # 钢 -> 环氧 I_B_A = interface(Z[1], Z[0]) # 环氧 -> 下一胞的钢 I_B_host = interface(Z[1], Z0) # 环氧 -> 宿主 P_A = propagation(omega * d[0] / c[0]) P_B = propagation(omega * d[1] / c[1]) Q = I_B_A @ P_B @ I_A_B @ P_A # 单胞 Bloch 矩阵 M_N = I_B_host @ P_B @ I_A_B @ P_A @ np.linalg.matrix_power(Q, N-1) @ I_host_A return Q, M_N # ==== 扫频 ==== freqs = np.linspace(100, 20000, 2000) omega = 2 * np.pi * freqs trans_rate = np.zeros_like(freqs) gap_flag = np.zeros_like(freqs) for i, om in enumerate(omega): Q, M = cell_transfer(om) t = 1.0 / M[1, 1] # 透射系数 trans_rate[i] = np.abs(t) ** 2 # 能量透射率 gap_flag[i] = 1.0 if np.abs(np.trace(Q)) > 2.0 else 0.0

逻辑说明：interface按左介质到右介质的顺序拼矩阵，这正是前面推导的 I(1→2)。propagation里的-1j对应 e^{+jωt} 约定下的右行波相位积累；如果你手上的项目用的是 e^{-iωt}，把两个指数符号整体互换即可，计算|t|^2的结果不变。cell_transfer里 Q 的乘法顺序举足轻重：必须先 P_A 再界面再 P_B，方向反了 trace 不变但 M_N 会错。最后一个关键点在M[1, 1]：M 是把左侧宿主里的振幅向量映射到右侧宿主里，右侧没有背向波，所以是[A_tr; 0] = M @ [A_in; B_ref]，解出来透射系数就是 M22 的倒数。

3.3 输出解读：色散曲线、禁带区间与传递率低谷对照

跑完这段代码，把trans_rate画成 dB 曲线，把gap_flag画成阴影区间叠在同一个图里。你会看到两个结果：在|trace(Q)| > 2的频率区间，传递率曲线往下掉得又快又深；在通带内，传递率曲线会来回波动，这是有限周期结构两端界面反射造成的驻波波纹，不代表有禁带。若在某个小区间gap_flag是 1 但传递率还在 -10 dB 附近，先怀疑透射公式取错矩阵元素，再怀疑两侧宿主阻抗不一致还硬用 1/M22。

判断禁带边界还有一个工程习惯：以传递率第一次穿过 -20 dB 的交点为准，不要直接把 trace 过 ±2 的点当成绝对边界。有限 N 的禁带边缘本来就有过渡区，N 越大过渡区越窄，谱上的边缘点向 trace 临界点靠近。所以在 v13 这类包里如果看到能带图和传递率图边界差了几个百分点，优先怀疑 N 设太小。

4. 参数怎么设：把传递率曲线调到想要的位置和深度

跑通只是第一步。实际做隔声或减振设计时，要回答的问题是“禁带能不能压到我的工作频段”，这就要学会调参数。

4.1 三个必调参数：密度比/声速比、厚度比、周期数

周期总厚度 a=d1+d2 是第一个要动的量。禁带中心频率大体上反比于 a：a 缩小一半，整段禁带向高频平移一倍左右。常见做法是先定 a 再定频率，别单独缩放某一层厚度，那样会同时改变周期和占空比，曲线变化不好归因。

厚度比 d1:d2 是第二个要动的量。保持 a 不变，把两层的相位厚度往接近的方向调，禁带宽度会变大；当两层的光学厚度 ωd/c 接近相等时，第一禁带一般做到最宽。这也是为什么很多声子晶体算例都选钢和环氧各一半：不是为了材料便宜，是为了把禁带调宽。

阻抗比 Z1/Z2 是第三个要动的量。密度比和声速比共同决定阻抗对比度，对比度越大，同一厚度比下的禁带越宽。钢/环氧的阻抗比大约在十几倍，出来的第一禁带已经能覆盖接近一个倍频程；想更宽就换橡胶或加截面变化，但阻抗比过大会让状态基矩阵的条件数变得很难看，这时更要坚持用幅值基。

4.2 带阻尼材料怎么写：复声速和复模量的落点

真实结构里树脂、橡胶这类材料必须有阻尼。常见做法是给阻尼层一个复声速或复弹性模量，比如环氧层写成：

eta = 0.02 c[1] = c[1] * (1 + 1j * eta)

在 e^{+jωt} 约定下，c 的虚部取正意味着波数 k 有负虚部，右行波 e^{-jkx} 沿 x 衰减，这是物理正确的阻尼方向。如果算完发现不加阻尼的系统透射率在某些频段算出了大于 1 的值，先查这一行虚部符号，而不是去调周期数。

加了阻尼后，禁带判断就不能再看|trace(Q)| > 2了，因为 trace 变成复数，临界条件模糊掉。更稳的判断是直接看 Q 的特征值模长：禁带内最大特征值模长明显大于 1，通带内约等于 1。这样也顺便能提取衰减常数，留到后面最终验证用。

4.3 用扫频结果反推有效参数：传递率低谷与 Bloch 衰减线对照

一个实用的自检是在扫频循环里顺手记录 Q 的最大特征值模长：

lam_max = np.max(np.abs(np.linalg.eigvals(Q))) alpha = np.log(lam_max) / a # 单位: Np/m

把 alpha 曲线叠在传递率曲线下面，你会发现传递率的深谷和 alpha 的高原严格对齐。若某段 alpha 明明大于 0 而传递率还在 -20 dB 左右晃荡，说明有限结构的边界反射把禁带里的衰减“短路”了；常见原因是两侧宿主介质选得太软，能量从表面漏过去了。这时该改的是宿主介质或入射端边界条件，而不是继续加周期数。

5. 声子晶体传递率计算的 5 个高频踩坑：现象、原因与解决

下面这五条是我在类似压缩包代码里反复见到的问题，每条都按现象、原因、解决三步写，排查起来可以直接照着走。

5.1 频率单位混乱：色散曲线整体平移 6.28 倍

现象是：能带图里第一个禁带的中心频率比自己手算的值大了约 2π 倍，或者反向光谱出现不存在的“鬼禁带”。原因是：扫频循环里用的是 Hz，但传递矩阵里的 k 需要的是波数，波数里的 ω 必须用角频率。直接把 f 塞进 φ 里，等效把声速放大了 6.28 倍。解决：统一在循环入口写omega = 2*np.pi*f；检查单层相位用phi = omega*d/c，不是f*d/c。

5.2 阻抗比悬殊时矩阵条件数爆炸：结果出现 NaN 或透射率大于 1

现象：钢和橡胶这类阻抗差两三数量级的组合，状态基写法算到 N=10 附近开始出现 NaN，或者明明是无源结构透射率却出现几十甚至上百。原因是：物理状态基的层矩阵里含有 Z 和 1/Z 作为离对角项，相乘几次后动态范围直接冲出浮点精度。解决：改用幅值基，也就是前面代码里的interface+propagation写法；如果必须保留状态基，先把状态向量做阻抗归一化，比如换成[√Z·u, σ/√Z]，再去做连乘。

5.3 直接把含宿主界面的 M1 做 N 次幂：透射率曲线出现假波纹

现象：传递率曲线在通带内多出一组 1/N 周期的细波纹，禁带深度反而变浅。原因是：有人图省事把单个胞含左右宿主界面的矩阵整体做幂次，而实际上相邻胞之间没有宿主层，中间只应该重复 Q 的内部界面。解决：按标准写法先搭内部单胞矩阵，再做matrix_power(Q, N-1)，最后只在两端各乘一个宿主界面矩阵。代码里我把这个顺序写死了，照着抄基本不会触发这条。

5.4 取错 M 矩阵元素：禁带位置对但透射深度全错

现象：gap_flag 标出来的禁带位置和传递率低谷大体能对上，但谷底只到 -6 dB 而不是预期的 -40 dB，或者反射率曲线长得像透射率。原因是：透射系数在幅值基下的公式是t = 1/M22，不是t = M22，也不是M11；不少人把反射系数和透射系数拿反了。另一个隐藏原因：两侧宿主介质阻抗不一致，M 已经不满足“右侧无背向波时 M22 的倒数就是透射系数”这个前提。解决：先拿一块均质板做 sanity check，令 A、B 材料参数完全相同，此时透射率应该是|exp(-jφ)|^2；如果对不上，说明 M 的界面向或取行列式有问题。

5.5 阻尼虚部符号反了：无源系统算出了能量增益

现象：加了复声速后，高频段透射率反而涨到 1 以上，尤其接近周期谐振频率时出现尖峰。原因是：时间约定和复模量虚部符号不匹配。e^{+jωt} 下 c 的虚部取正代表衰减，取负代表介质主动供能；e^{-iωt} 下则反过来。解决：先在一个单层均质板里试，令入射和出射介质完全一样，透射率应严格小于等于 1；如果大于 1，把1+1j*eta改成1-1j*eta重试。

6. 进阶：用有限周期透射率曲线反推 Bloch 衰减常数，多一条验证手段

参数扫完、坑也踩过一轮，我习惯再加一道验证：把 Q 矩阵的特征值拉开，画一条衰减常数曲线，用来判断有限周期透射率低谷是“真禁带”还是“两端边界阻抗不匹配的假象”。这条曲线不依赖 N，只依赖单胞材料和几何，所以它比直接看传递率更深一层。

实现上在扫频循环里补三行：

eigs = np.abs(np.linalg.eigvals(Q)) lam_max = np.max(eigs) alpha[i] = np.log(lam_max) / (d[0] + d[1]) # Np/m

通带内lam_max约等于 1，alpha 接近 0；禁带内 lam_max 明显大于 1，alpha 变成正值。把 alpha 大于 0 的区间和gap_flag叠在一起，两者应该重合。重合不上的时候，问题多半出在 Q 的构造顺序，而不是材料参数——调材料是白费功夫。

这个 alpha 还能用来估算周期数对隔声量的贡献。每增加一个胞，禁带中心频率附近的隔离度增量约为20*log10(lam_max)dB。比如用钢/环氧 20mm/20mm，N=8 和 N=16 在禁带中心的传递率差，应该约等于 8 倍的这个值。如果实际算出来差太多，去检查 M 的朝向，而不是怀疑声速数据少了个零。注意这个估算只在离开带边界足够远的频率才准，带边界附近有限结构的透射深谷还受两端界面反射影响，衰减常数反而不如表面阻抗条件重要。

我拿到这类压缩包，第一件事就是把 trace(Q) 的临界边界、alpha 曲线和有限周期传递率的 -20 dB 交点叠在一张图上，三线对不上就先查矩阵朝向，而不是急着调厚度。希望帮到你。

本文还有配套的精品资源，点击获取