主题049：电磁兼容与干扰分析

目录

1. 电磁兼容基础理论
2. 电磁干扰源与耦合机制
3. 传导干扰与辐射干扰
4. 屏蔽技术
5. 滤波与接地技术
6. EMC测试与标准
7. Python仿真实现

---

一、电磁兼容基础理论

1.1 电磁兼容的基本概念

电磁兼容（Electromagnetic Compatibility, EMC）是指设备或系统在其电磁环境中能正常工作，且不对该环境中任何事物构成不能承受的电磁骚扰的能力。电磁兼容包含两个方面的要求：

1. 电磁干扰（EMI）：设备或系统在正常运行过程中产生的电磁能量，对其他设备或系统造成的不良影响
2. 电磁敏感性（EMS）：设备或系统抵抗外部电磁干扰的能力，也称为电磁抗扰度

电磁兼容的三要素模型：

$$干扰源\stackrel{耦合路径}{\to }敏感设备$$

要实现电磁兼容，必须对这三个要素进行控制：

• 抑制干扰源：降低干扰源的发射电平
• 切断耦合路径：采用屏蔽、滤波、接地等措施
• 提高敏感设备的抗扰度：增强设备的电磁抗干扰能力

1.2 电磁兼容的重要性

随着电子设备的广泛应用和集成度的提高，电磁兼容问题日益突出：

1. 设备密集度增加：现代电子系统中设备数量大幅增加，相互干扰风险增大
2. 工作频率提高：高速数字电路和射频电路产生更丰富的谐波分量
3. 灵敏度提高：低电压、低功耗器件对电磁干扰更加敏感
4. 法规要求：各国都制定了严格的EMC法规标准，产品必须通过认证

EMC问题的影响：

• 设备性能下降或功能失效
• 数据错误或丢失
• 系统崩溃或损坏
• 安全隐患（如医疗设备、航空电子）

1.3 频谱管理与EMC标准

电磁频谱是有限的资源，需要进行有效管理。国际电信联盟（ITU）负责全球频谱分配，各国设有相应的管理机构（如FCC在美国、CE在欧洲、SRRC在中国）。

主要EMC标准体系：

1. CISPR标准（国际无线电干扰特别委员会）：

   • CISPR 11：工业、科学和医疗（ISM）设备
   • CISPR 22/32：信息技术设备
   • CISPR 25：车辆、船舶和内燃机

2. IEC标准（国际电工委员会）：

   • IEC 61000系列：电磁兼容性基础标准
   • IEC 61000-4-x：抗扰度测试方法
   • IEC 61000-3-x：发射限值

3. 军用标准：

   • MIL-STD-461/462：美国军用EMC标准
   • GJB 151/152：中国军用EMC标准

1.4 EMC设计流程

电磁兼容设计应贯穿产品开发的整个生命周期：

1. 需求分析阶段：

• 确定产品的EMC标准要求
• 识别潜在的干扰源和敏感电路
• 制定EMC设计目标和预算

2. 设计阶段：

• 电路设计：选择合适的器件、降低di/dt和dv/dt
• PCB设计：合理的层叠结构、布线策略、去耦设计
• 结构设计：屏蔽设计、接地设计、电缆布局

3. 仿真分析阶段：

• 信号完整性（SI）分析
• 电源完整性（PI）分析
• 电磁干扰（EMI）仿真

4. 测试验证阶段：

• 预兼容测试
• 正式EMC测试
• 问题诊断与整改

---

二、电磁干扰源与耦合机制

2.1 自然干扰源

1. 雷电：

• 电流幅值：数千至数十万安培
• 上升时间：微秒级
• 频谱范围：DC至数百MHz
• 影响：直接击中、感应耦合、地电位升高

2. 静电放电（ESD）：

• 电压幅值：数千至数万伏
• 上升时间：亚纳秒级
• 频谱范围：DC至GHz
• 影响：设备损坏、数据错误、系统复位

3. 太阳活动和宇宙噪声：

• 太阳耀斑、射电爆发
• 银河背景噪声
• 主要影响卫星通信和射电天文

2.2 人为干扰源

1. 电力系统干扰：

• 开关操作：断路器、接触器、继电器
• 负载切换：电机启动、电容器投切
• 谐波：整流器、变频器、电弧炉
• 电压暂降/中断

2. 数字电路干扰：

• 时钟谐波：高速时钟的各次谐波
• 开关噪声：逻辑电平跳变产生的宽带噪声
• 同步开关噪声（SSN）：多路信号同时切换

3. 射频设备：

• 广播发射机：AM/FM/TV
• 移动通信：手机基站、对讲机
• 雷达系统
• 工业、科学和医疗（ISM）设备

4. 汽车电子：

• 点火系统：火花塞放电
• 电机驱动：雨刮器、电动车窗
• 燃油喷射系统

2.3 耦合机制

电磁干扰从干扰源传输到敏感设备主要通过四种耦合方式：

1. 传导耦合（Conducted Coupling）：

干扰通过导线、电缆、PCB走线等导体传输。

• 共模传导：干扰电流在相线和中线上同向流动
• 差模传导：干扰电流在相线和中线上反向流动

2. 电场耦合（Electric Field Coupling）：

通过电容耦合，干扰源的电场在敏感电路上感应出电压。

耦合电压：

$$V_{coupled}=j\omega C_{mutual}{Z}_{victim}{V}_{source}$$

其中$C_{mutual}$为互电容，$Z_{victim}$为受害电路阻抗。

减小电场耦合的方法：

• 增加间距
• 使用屏蔽层
• 降低源电压或受害电路阻抗

3. 磁场耦合（Magnetic Field Coupling）：

通过互感耦合，干扰源的磁场在敏感电路上感应出电压。

感应电压：

$$V_{induced}=j\omega M{I}_{source}$$

其中$M$为互感，$I_{source}$为干扰源电流。

减小磁场耦合的方法：

• 增加间距
• 使用双绞线或同轴电缆
• 减小环路面积
• 使用高磁导率屏蔽材料

4. 辐射耦合（Radiated Coupling）：

干扰源以电磁波形式辐射能量，被敏感设备接收。

近场/远场分界距离：

$$r=\frac{\lambda }{2\pi }\approx \frac{47.7}{f(MHz)}\text{ (m)}$$

在近场区：

• 高阻抗源（如偶极子天线）主要产生电场
• 低阻抗源（如环形天线）主要产生磁场

在远场区：电场和磁场以固定比例（377Ω）存在。

---

三、传导干扰与辐射干扰

3.1 传导干扰分析

共模干扰与差模干扰：

共模干扰是指同时出现在相线和中线（或信号线和地线）上、相对于地的同相干扰电压。差模干扰是指出现在相线和中线之间、与工作信号串联的干扰电压。

传导干扰的频谱特性：

• 低频段（<1MHz）：主要是差模干扰
• 高频段（>1MHz）：共模干扰占主导

传导干扰的测量：

使用线路阻抗稳定网络（LISN）进行传导发射测量：

$$V_{measured}=V_{CM}+V_{DM}$$

通过LISN可以将共模和差模干扰分离，分别进行测量和分析。

3.2 辐射干扰分析

辐射发射的来源：

1. 差模辐射：由信号环路产生

   辐射电场：

   $$E=\frac{1.316\times 10^{-14}{f}^{2}AI}{r}$$

   其中$A$为环路面积，$I$为电流，$f$为频率，$r$为距离

2. 共模辐射：由共模电流在电缆上产生

   辐射电场：

   $$E=\frac{6.28\times 10^{-7}fL{I}_{CM}}{r}$$

   其中$L$为电缆长度，$I_{CM}$为共模电流

减小辐射发射的方法：

• 减小信号环路面积
• 控制共模电流
• 使用屏蔽电缆
• 合理布局，避免天线效应

3.3 传输线理论在EMC中的应用

传输线模型：

对于高频信号，PCB走线和电缆需要用传输线理论分析：

特性阻抗：

$$Z_0=\sqrt{\frac{L}{C}}$$

传播速度：

$$v=\frac{1}{\sqrt{LC}}=\frac{c}{\sqrt{{\epsilon }_r}}$$

反射与匹配：

当传输线终端阻抗不匹配时，会产生反射：

$$\Gamma =\frac{Z_L-Z_0}{Z_L+Z_0}$$

反射会导致信号完整性问题和额外的辐射发射。

串扰分析：

相邻传输线之间的串扰：

• 近端串扰（NEXT）：

  $$NEXT=\frac{V_{near}}{V_{source}}=\frac{1}{4}\left(\frac{C_m}{C}+\frac{L_m}{L}\right)$$

• 远端串扰（FEXT）：

  $$FEXT=\frac{V_{far}}{V_{source}}=\frac{l}{2t_r}\left(\frac{C_m}{C}-\frac{L_m}{L}\right)$$

其中$C_m$、$L_m$为互容和互感，$l$为耦合长度，$t_r$为上升时间。

---

四、屏蔽技术

4.1 屏蔽原理

屏蔽是利用导电或导磁材料制成的屏蔽体，将电磁能量限制在特定空间内或阻止外部电磁能量进入。

屏蔽效能（Shielding Effectiveness, SE）：

$$SE=20{\log }_{10}\left(\frac{E_{without}}{E_{with}}\right)\text{ (dB)}$$

或

$$SE=20{\log }_{10}\left(\frac{H_{without}}{H_{with}}\right)\text{ (dB)}$$

屏蔽效能通常由三部分组成：

$$SE=A+R+B$$

其中：

• $A$：吸收损耗
• $R$：反射损耗
• $B$：多次反射修正项（通常可忽略）

4.2 吸收损耗

电磁波在屏蔽材料中传播时的衰减：

$$A=8.686\frac{t}{\delta }=8.686t\sqrt{\pi f\mu \sigma }\text{ (dB)}$$

其中：

• $t$：屏蔽材料厚度
• $\delta$：趋肤深度
• $f$：频率
• $\mu$：磁导率
• $\sigma$：电导率

趋肤深度：

$$\delta =\frac{1}{\sqrt{\pi f\mu \sigma }}$$

对于铜（$\sigma =5.8\times 10^7$ S/m）：

频率	趋肤深度
1 kHz	2.1 mm
1 MHz	66 μm
100 MHz	6.6 μm
1 GHz	2.1 μm

4.3 反射损耗

电磁波在屏蔽材料表面反射造成的损耗：

对于电场（平面波）：

$$R_E=168+10{\log }_{10}\left(\frac{{\sigma }_r}{{\mu }_r{f}^3{r}^2}\right)\text{ (dB)}$$

对于磁场（平面波）：

$$R_H=168+10{\log }_{10}\left(\frac{f{r}^2{\sigma }_r}{{\mu }_r}\right)\text{ (dB)}$$

其中：

• ${\sigma }_r$：相对电导率
• ${\mu }_r$：相对磁导率
• $r$：源到屏蔽体的距离

重要结论：

• 低频磁场最难屏蔽（反射损耗小、趋肤深度大）
• 高磁导率材料（如坡莫合金、铁氧体）对磁场屏蔽更有效

4.4 屏蔽材料选择

常用屏蔽材料：

材料	相对电导率	相对磁导率	特点
铜	1.0	1	电导率高，加工性好
铝	0.61	1	重量轻，成本低
钢	0.17	1000	磁屏蔽效果好
坡莫合金	0.04	80000	高磁导率，低频磁屏蔽
铁氧体	0.01	1500	高频磁屏蔽，成本低

屏蔽设计要点：

1. 缝隙和孔洞：

   • 缝隙长度应小于λ/20（最高工作频率）
   • 使用导电衬垫、指形簧片、导电胶
   • 螺钉间距应小于λ/20

2. 通风孔：

   • 使用波导通风板
   • 截止波导长度：$l>3d$
   • 截止频率：$f_c=\frac{17.5}{d}$（GHz，d单位为cm）

3. 电缆穿透：

   • 使用屏蔽电缆
   • 360°端接
   • 使用馈通滤波器

---

五、滤波与接地技术

5.1 EMI滤波器设计

EMI滤波器用于抑制传导干扰，通常采用低通滤波器结构。

基本滤波电路：

1. 电容滤波：并联在干扰源和地之间，提供高频旁路
2. 电感滤波：串联在干扰路径中，提供高频阻抗
3. LC滤波：组合使用电感和电容
4. π型滤波：两个电容中间一个电感
5. T型滤波：两个电感中间一个电容

滤波器设计要点：

• 截止频率：$f_c=\frac{1}{2\pi \sqrt{LC}}$
• 阻抗匹配：源阻抗和负载阻抗决定滤波器拓扑
• 高频性能：注意元件的寄生参数

EMI滤波器与信号滤波器的区别：

特性	EMI滤波器	信号滤波器
目的	抑制干扰	保留有用信号
频率范围	宽频带	特定频段
插入损耗	越大越好	根据需求设计
功率容量	大	小

5.2 接地技术

接地是EMC设计中最重要但最容易被忽视的技术。

接地的目的：

1. 安全接地：防止触电事故
2. 信号参考：为电路提供零电位参考
3. EMI控制：提供干扰电流的返回路径
4. 屏蔽接地：保证屏蔽效能

接地方式：

1. 单点接地：

   • 适用于低频电路（<1MHz）
   • 避免地环路
   • 串联单点接地和并联单点接地

2. 多点接地：

   • 适用于高频电路（>10MHz）
   • 减小接地阻抗
   • 注意地环路问题

3. 混合接地：

   • 电容耦合单点接地
   • 低频单点、高频多点

地阻抗：

$$Z_g=R_g+j\omega L_g$$

对于典型PCB地平面：

• 直流电阻：约0.1-1 mΩ/□
• 电感：约0.5-1 nH/mm

在100MHz时，1cm长的地线电感（约5nH）产生的阻抗：

$$X_L=2\pi fL=2\pi \times 10^8\times 5\times 10^{-9}=3.14\text{ Ω}$$

5.3 去耦设计

去耦电容用于为集成电路提供本地电荷存储，抑制电源噪声。

去耦电容的选择：

1. 大容量电容（10-100μF）：

   • 电解电容或钽电容
   • 提供低频能量存储
   • 每个电源输入端一个

2. 中容量电容（0.1-1μF）：

   • 陶瓷电容
   • 每个IC一个
   • 尽可能靠近IC电源引脚

3. 小容量电容（10-1000pF）：

   • 高频陶瓷电容
   • 抑制高频噪声
   • 用于高速IC

去耦电容的放置：

• 电容到IC的引线长度应最小化
• 使用多个不同容值的电容并联
• 注意电容的自谐振频率

去耦效果分析：

去耦电容与电源平面形成的LC谐振：

$$f_{res}=\frac{1}{2\pi \sqrt{L_{esr}C}}$$

其中$L_{esr}$包括电容ESL和走线电感。

---

六、EMC测试与标准

6.1 EMC测试项目

发射测试（Emission Tests）：

1. 传导发射（CE）：

   • 频率范围：150kHz - 30MHz
   • 使用LISN测量电源线上的干扰电压
   • 分为共模和差模测量

2. 辐射发射（RE）：

   • 频率范围：30MHz - 1GHz（或更高）
   • 在电波暗室中使用天线测量
   • 测量距离：3m、10m或30m

抗扰度测试（Immunity Tests）：

1. 静电放电（ESD）：

   • 接触放电：±4kV至±8kV
   • 空气放电：±8kV至±15kV

2. 射频电磁场辐射抗扰度（RS）：

   • 频率范围：80MHz - 6GHz
   • 场强：3V/m至30V/m

3. 电快速瞬变脉冲群（EFT）：

   • 脉冲幅度：±0.5kV至±4kV
   • 重复频率：5kHz或100kHz

4. 浪涌（Surge）：

   • 开路电压：±0.5kV至±4kV
   • 短路电流：相应匹配

5. 传导抗扰度（CS）：

   • 频率范围：150kHz - 80MHz
   • 注入电流：1mA至100mA

6. 电压暂降和短时中断（Dips）：

   • 暂降：0%至70%额定电压
   • 中断：0%额定电压

6.2 测试设施

电波暗室（Anechoic Chamber）：

• 屏蔽效能：>100dB
• 吸波材料：锥形或楔形吸波材料
• 归一化场地衰减（NSA）：±4dB

开阔试验场（OATS）：

• 无障碍的平坦场地
• 金属接地平面
• 环境电平低于限值6dB

TEM/GTEM小室：

• 用于小设备的辐射测试
• 频率范围：DC至数GHz
• 场均匀性好

6.3 预兼容测试

在产品开发阶段进行预兼容测试，可以及早发现和解决EMC问题。

预兼容测试设备：

• 频谱分析仪
• 近场探头（电场和磁场）
• LISN
• 天线（双锥、对数周期、喇叭）

近场探测技术：

使用近场探头定位干扰源：

• 电场探头：定位高电压点
• 磁场探头：定位高电流环路
• 扫描PCB找出热点区域

---

七、Python仿真实现

本节提供电磁兼容与干扰分析的Python仿真代码，包括屏蔽效能计算、滤波器设计、串扰分析和辐射发射估算。

仿真代码详见同目录下的 run_simulation.py 文件。

"""
主题049：电磁兼容与干扰分析 - Python仿真
============================================

本程序实现以下仿真内容：
1. 屏蔽效能分析
2. EMI滤波器设计
3. 传输线串扰分析
4. 辐射发射估算
5. 接地阻抗分析

作者：电磁兼容仿真专家
日期：2026-03-06
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle, FancyArrowPatch
from matplotlib.collections import LineCollection
import matplotlib.patches as mpatches

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 使用Agg后端，不弹出窗口
plt.switch_backend('Agg')

print("=" * 70)
print("主题049：电磁兼容与干扰分析 - Python仿真")
print("=" * 70)

# ==============================================================================
# 1. 屏蔽效能分析
# ==============================================================================
print("\n[1/5] 屏蔽效能分析...")

def skin_depth(f, sigma, mu_r=1):
    """
    计算趋肤深度
    
    参数:
        f: 频率 (Hz)
        sigma: 电导率 (S/m)
        mu_r: 相对磁导率
    
    返回:
        趋肤深度 (m)
    """
    mu_0 = 4 * np.pi * 1e-7  # 真空磁导率
    mu = mu_r * mu_0
    return 1 / np.sqrt(np.pi * f * mu * sigma)

def absorption_loss(f, t, sigma, mu_r=1):
    """
    计算吸收损耗
    
    参数:
        f: 频率 (Hz)
        t: 屏蔽材料厚度 (m)
        sigma: 电导率 (S/m)
        mu_r: 相对磁导率
    
    返回:
        吸收损耗 (dB)
    """
    delta = skin_depth(f, sigma, mu_r)
    return 8.686 * t / delta

def reflection_loss_electric(f, sigma_r, mu_r, r=1):
    """
    计算电场反射损耗
    
    参数:
        f: 频率 (Hz)
        sigma_r: 相对电导率
        mu_r: 相对磁导率
        r: 源到屏蔽体距离 (m)
    
    返回:
        反射损耗 (dB)
    """
    return 168 + 10 * np.log10(sigma_r / (mu_r * f**3 * r**2))

def reflection_loss_magnetic(f, sigma_r, mu_r, r=1):
    """
    计算磁场反射损耗
    
    参数:
        f: 频率 (Hz)
        sigma_r: 相对电导率
        mu_r: 相对磁导率
        r: 源到屏蔽体距离 (m)
    
    返回:
        反射损耗 (dB)
    """
    return 168 + 10 * np.log10(f * r**2 * sigma_r / mu_r)

# 材料参数
materials = {
    '铜': {'sigma': 5.8e7, 'mu_r': 1, 'sigma_r': 1.0, 'color': '#B87333'},
    '铝': {'sigma': 3.5e7, 'mu_r': 1, 'sigma_r': 0.61, 'color': '#C0C0C0'},
    '钢': {'sigma': 1.0e7, 'mu_r': 1000, 'sigma_r': 0.17, 'color': '#708090'},
    '坡莫合金': {'sigma': 2.0e6, 'mu_r': 80000, 'sigma_r': 0.04, 'color': '#4A4A4A'}
}

# 频率范围
f = np.logspace(3, 10, 500)  # 1 kHz 到 10 GHz
t = 1e-3  # 1 mm 厚度

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (a) 趋肤深度
ax1 = axes[0, 0]
for name, props in materials.items():
    delta = skin_depth(f, props['sigma'], props['mu_r'])
    ax1.loglog(f/1e6, delta*1e6, linewidth=2, label=name, color=props['color'])

ax1.set_xlabel('频率 (MHz)', fontsize=11)
ax1.set_ylabel('趋肤深度 (μm)', fontsize=11)
ax1.set_title('(a) 不同材料的趋肤深度', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right', fontsize=9)
ax1.grid(True, which='both', linestyle='--', alpha=0.5)
ax1.set_xlim([1e-3, 1e4])

# (b) 吸收损耗
ax2 = axes[0, 1]
for name, props in materials.items():
    A = absorption_loss(f, t, props['sigma'], props['mu_r'])
    ax2.semilogx(f/1e6, A, linewidth=2, label=name, color=props['color'])

ax2.set_xlabel('频率 (MHz)', fontsize=11)
ax2.set_ylabel('吸收损耗 (dB)', fontsize=11)
ax2.set_title('(b) 吸收损耗 (厚度=1mm)', fontsize=12, fontweight='bold')
ax2.legend(loc='lower right', fontsize=9)
ax2.grid(True, which='both', linestyle='--', alpha=0.5)
ax2.set_xlim([1e-3, 1e4])

# (c) 电场反射损耗
ax3 = axes[1, 0]
r = 1  # 1m 距离
for name, props in materials.items():
    R_E = reflection_loss_electric(f, props['sigma_r'], props['mu_r'], r)
    ax3.semilogx(f/1e6, R_E, linewidth=2, label=name, color=props['color'])

ax3.set_xlabel('频率 (MHz)', fontsize=11)
ax3.set_ylabel('电场反射损耗 (dB)', fontsize=11)
ax3.set_title('(c) 电场反射损耗', fontsize=12, fontweight='bold')
ax3.legend(loc='upper right', fontsize=9)
ax3.grid(True, which='both', linestyle='--', alpha=0.5)
ax3.set_xlim([1e-3, 1e4])

# (d) 总屏蔽效能
ax4 = axes[1, 1]
for name, props in materials.items():
    A = absorption_loss(f, t, props['sigma'], props['mu_r'])
    R_E = reflection_loss_electric(f, props['sigma_r'], props['mu_r'], r)
    SE = A + R_E
    ax4.semilogx(f/1e6, SE, linewidth=2, label=name, color=props['color'])

ax4.axhline(y=60, color='red', linestyle='--', linewidth=1.5, label='60dB目标')
ax4.set_xlabel('频率 (MHz)', fontsize=11)
ax4.set_ylabel('总屏蔽效能 (dB)', fontsize=11)
ax4.set_title('(d) 总屏蔽效能', fontsize=12, fontweight='bold')
ax4.legend(loc='upper right', fontsize=9)
ax4.grid(True, which='both', linestyle='--', alpha=0.5)
ax4.set_xlim([1e-3, 1e4])

plt.tight_layout()
plt.savefig('shielding_effectiveness.png', dpi=150, bbox_inches='tight')
print("  ✓ shielding_effectiveness.png 已保存")
plt.close()

# ==============================================================================
# 2. EMI滤波器设计
# ==============================================================================
print("\n[2/5] EMI滤波器设计...")

def lc_filter_response(f, L, C, R_source=50, R_load=50):
    """
    LC低通滤波器频率响应
    
    参数:
        f: 频率数组 (Hz)
        L: 电感值 (H)
        C: 电容值 (F)
        R_source: 源阻抗 (Ω)
        R_load: 负载阻抗 (Ω)
    
    返回:
        插入损耗 (dB)
    """
    omega = 2 * np.pi * f
    
    # LC滤波器传递函数
    Z_L = 1j * omega * L
    Z_C = 1 / (1j * omega * C)
    
    # 分压计算
    Z_parallel = (Z_C * R_load) / (Z_C + R_load)
    H = Z_parallel / (R_source + Z_L + Z_parallel)
    
    # 插入损耗
    H_no_filter = R_load / (R_source + R_load)
    IL = 20 * np.log10(np.abs(H_no_filter / H))
    
    return IL

def pi_filter_response(f, L, C1, C2, R_source=50, R_load=50):
    """
    π型滤波器频率响应
    
    参数:
        f: 频率数组 (Hz)
        L: 电感值 (H)
        C1, C2: 电容值 (F)
        R_source: 源阻抗 (Ω)
        R_load: 负载阻抗 (Ω)
    
    返回:
        插入损耗 (dB)
    """
    omega = 2 * np.pi * f
    
    Z_L = 1j * omega * L
    Z_C1 = 1 / (1j * omega * C1)
    Z_C2 = 1 / (1j * omega * C2)
    
    # 等效阻抗
    Z_in = 1 / (1/Z_C1 + 1/(Z_L + 1/(1/Z_C2 + 1/R_load)))
    
    H = (1/(1/Z_C2 + 1/R_load)) / (Z_L + 1/(1/Z_C2 + 1/R_load))
    H = H * Z_in / (R_source + Z_in)
    
    H_no_filter = R_load / (R_source + R_load)
    IL = 20 * np.log10(np.abs(H_no_filter / H))
    
    return IL

# 频率范围
f_filter = np.logspace(3, 9, 500)  # 1 kHz 到 1 GHz

# 滤波器参数
L = 100e-6  # 100 μH
C = 1e-6    # 1 μF
C1 = C2 = 0.47e-6  # 0.47 μF

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (a) LC滤波器电路图
ax1 = axes[0, 0]
ax1.set_xlim(0, 10)
ax1.set_ylim(0, 10)
ax1.axis('off')
ax1.set_title('(a) LC低通滤波器电路', fontsize=12, fontweight='bold')

# 绘制电路
# 源
ax1.add_patch(FancyBboxPatch((0.5, 4), 1.5, 2, boxstyle="round,pad=0.1", 
                              facecolor='#FFE4B5', edgecolor='black', linewidth=2))
ax1.text(1.25, 5, 'V_s', ha='center', va='center', fontsize=11, fontweight='bold')

# 电感
ax1.plot([2, 3], [5, 5], 'k-', linewidth=2)
for i in range(5):
    x = 3 + i * 0.3
    ax1.plot([x, x+0.15], [5, 5.3], 'k-', linewidth=2)
    ax1.plot([x+0.15, x+0.3], [5.3, 5], 'k-', linewidth=2)
ax1.plot([4.5, 5], [5, 5], 'k-', linewidth=2)
ax1.text(3.75, 5.8, 'L', ha='center', fontsize=11, fontweight='bold')

# 电容
ax1.plot([5, 5], [5, 3.5], 'k-', linewidth=2)
ax1.plot([4.5, 5.5], [3.5, 3.5], 'k-', linewidth=2)
ax1.plot([4.5, 5.5], [3.2, 3.2], 'k-', linewidth=2)
ax1.plot([5, 5], [3.2, 2], 'k-', linewidth=2)
ax1.text(6, 3.35, 'C', ha='left', fontsize=11, fontweight='bold')

# 负载
ax1.plot([5, 7], [5, 5], 'k-', linewidth=2)
ax1.plot([7, 7], [5, 3.5], 'k-', linewidth=2)
ax1.add_patch(FancyBboxPatch((6.5, 2.5), 1, 1, boxstyle="round,pad=0.05",
                              facecolor='#E6E6FA', edgecolor='black', linewidth=2))
ax1.text(7, 3, 'R_L', ha='center', va='center', fontsize=10)
ax1.plot([7, 7], [2.5, 2], 'k-', linewidth=2)

# 地线
ax1.plot([2, 7], [2, 2], 'k-', linewidth=2)
ax1.plot([1.5, 7.5], [2, 2], 'k-', linewidth=1)

# (b) 不同电感值的LC滤波器
ax2 = axes[0, 1]
L_values = [10e-6, 47e-6, 100e-6, 220e-6]  # 10, 47, 100, 220 μH
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']

for L_val, color in zip(L_values, colors):
    IL = lc_filter_response(f_filter, L_val, C)
    ax2.semilogx(f_filter/1e6, IL, linewidth=2, 
                 label=f'L={L_val*1e6:.0f}μH', color=color)

ax2.axhline(y=40, color='red', linestyle='--', linewidth=1.5, alpha=0.7)
ax2.set_xlabel('频率 (MHz)', fontsize=11)
ax2.set_ylabel('插入损耗 (dB)', fontsize=11)
ax2.set_title('(b) LC滤波器插入损耗 (C=1μF)', fontsize=12, fontweight='bold')
ax2.legend(loc='lower right', fontsize=9)
ax2.grid(True, which='both', linestyle='--', alpha=0.5)
ax2.set_xlim([0.001, 1000])

# (c) π型滤波器
ax3 = axes[1, 0]
ax3.set_xlim(0, 10)
ax3.set_ylim(0, 10)
ax3.axis('off')
ax3.set_title('(c) π型滤波器电路', fontsize=12, fontweight='bold')

# 绘制π型滤波器电路
# 源
ax3.add_patch(FancyBboxPatch((0.5, 4), 1.5, 2, boxstyle="round,pad=0.1",
                              facecolor='#FFE4B5', edgecolor='black', linewidth=2))
ax3.text(1.25, 5, 'V_s', ha='center', va='center', fontsize=11, fontweight='bold')

# C1
ax3.plot([2, 3], [5, 5], 'k-', linewidth=2)
ax3.plot([3, 3], [5, 3.5], 'k-', linewidth=2)
ax3.plot([2.5, 3.5], [3.5, 3.5], 'k-', linewidth=2)
ax3.plot([2.5, 3.5], [3.2, 3.2], 'k-', linewidth=2)
ax3.plot([3, 3], [3.2, 2], 'k-', linewidth=2)
ax3.text(4, 3.35, 'C1', ha='left', fontsize=10, fontweight='bold')

# L
ax3.plot([3, 4], [5, 5], 'k-', linewidth=2)
for i in range(5):
    x = 4 + i * 0.3
    ax3.plot([x, x+0.15], [5, 5.3], 'k-', linewidth=2)
    ax3.plot([x+0.15, x+0.3], [5.3, 5], 'k-', linewidth=2)
ax3.plot([5.5, 6], [5, 5], 'k-', linewidth=2)
ax3.text(4.75, 5.8, 'L', ha='center', fontsize=11, fontweight='bold')

# C2
ax3.plot([6, 6], [5, 3.5], 'k-', linewidth=2)
ax3.plot([5.5, 6.5], [3.5, 3.5], 'k-', linewidth=2)
ax3.plot([5.5, 6.5], [3.2, 3.2], 'k-', linewidth=2)
ax3.plot([6, 6], [3.2, 2], 'k-', linewidth=2)
ax3.text(7, 3.35, 'C2', ha='left', fontsize=10, fontweight='bold')

# 负载
ax3.plot([6, 8], [5, 5], 'k-', linewidth=2)
ax3.plot([8, 8], [5, 3.5], 'k-', linewidth=2)
ax3.add_patch(FancyBboxPatch((7.5, 2.5), 1, 1, boxstyle="round,pad=0.05",
                              facecolor='#E6E6FA', edgecolor='black', linewidth=2))
ax3.text(8, 3, 'R_L', ha='center', va='center', fontsize=10)
ax3.plot([8, 8], [2.5, 2], 'k-', linewidth=2)

# 地线
ax3.plot([2, 8], [2, 2], 'k-', linewidth=2)
ax3.plot([1.5, 8.5], [2, 2], 'k-', linewidth=1)

# (d) π型滤波器性能对比
ax4 = axes[1, 1]
IL_lc = lc_filter_response(f_filter, L, C)
IL_pi = pi_filter_response(f_filter, L, C1, C2)

ax4.semilogx(f_filter/1e6, IL_lc, linewidth=2.5, label='LC滤波器', color='#1f77b4')
ax4.semilogx(f_filter/1e6, IL_pi, linewidth=2.5, label='π型滤波器', color='#ff7f0e')

ax4.axhline(y=60, color='red', linestyle='--', linewidth=1.5, alpha=0.7, label='60dB目标')
ax4.set_xlabel('频率 (MHz)', fontsize=11)
ax4.set_ylabel('插入损耗 (dB)', fontsize=11)
ax4.set_title('(d) LC与π型滤波器对比', fontsize=12, fontweight='bold')
ax4.legend(loc='lower right', fontsize=10)
ax4.grid(True, which='both', linestyle='--', alpha=0.5)
ax4.set_xlim([0.001, 1000])

plt.tight_layout()
plt.savefig('filter_design.png', dpi=150, bbox_inches='tight')
print("  ✓ filter_design.png 已保存")
plt.close()

# ==============================================================================
# 3. 传输线串扰分析
# ==============================================================================
print("\n[3/5] 传输线串扰分析...")

def crosstalk_near_end(Cm, C, Lm, L, l, tr):
    """
    计算近端串扰
    
    参数:
        Cm: 互电容 (F/m)
        C: 自电容 (F/m)
        Lm: 互感 (H/m)
        L: 自感 (H/m)
        l: 耦合长度 (m)
        tr: 上升时间 (s)
    
    返回:
        近端串扰系数
    """
    return 0.25 * (Cm/C + Lm/L)

def crosstalk_far_end(Cm, C, Lm, L, l, tr, v):
    """
    计算远端串扰
    
    参数:
        Cm: 互电容 (F/m)
        C: 自电容 (F/m)
        Lm: 互感 (H/m)
        L: 自感 (H/m)
        l: 耦合长度 (m)
        tr: 上升时间 (s)
        v: 传播速度 (m/s)
    
    返回:
        远端串扰系数
    """
    td = l / v  # 传输延迟
    return (td / tr) * (Cm/C - Lm/L)

# 传输线参数
L = 400e-9    # 400 nH/m
C = 40e-12    # 40 pF/m
Lm_values = np.linspace(0, 200e-9, 100)  # 0 到 200 nH/m
Cm_values = np.linspace(0, 20e-12, 100)  # 0 到 20 pF/m

# 耦合长度和上升时间
l = 0.3  # 30 cm
tr = 1e-9  # 1 ns
v = 1.5e8  # 传播速度 (m/s), 约0.5c

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (a) 串扰耦合模型
ax1 = axes[0, 0]
ax1.set_xlim(0, 10)
ax1.set_ylim(0, 10)
ax1.axis('off')
ax1.set_title('(a) 传输线串扰耦合模型', fontsize=12, fontweight='bold')

# 绘制两条传输线
# 攻击线
ax1.plot([1, 9], [6, 6], 'b-', linewidth=4, label='攻击线')
ax1.plot([1, 9], [5.5, 5.5], 'b-', linewidth=4)
ax1.fill_between([1, 9], [5.5, 5.5], [6, 6], alpha=0.3, color='blue')

# 受害线
ax1.plot([1, 9], [3.5, 3.5], 'r-', linewidth=4, label='受害线')
ax1.plot([1, 9], [3, 3], 'r-', linewidth=4)
ax1.fill_between([1, 9], [3, 3], [3.5, 3.5], alpha=0.3, color='red')

# 耦合标注
for i in range(5):
    x = 2 + i * 1.5
    # 电容耦合
    ax1.annotate('', xy=(x, 5.5), xytext=(x, 3.5),
                arrowprops=dict(arrowstyle='<->', color='green', lw=1.5))
    # 电感耦合
    circle1 = plt.Circle((x, 5.75), 0.15, fill=False, color='purple', linewidth=1.5)
    circle2 = plt.Circle((x, 3.25), 0.15, fill=False, color='purple', linewidth=1.5)
    ax1.add_patch(circle1)
    ax1.add_patch(circle2)

ax1.text(5, 4.5, 'Cm (互容)', ha='center', fontsize=10, color='green', fontweight='bold')
ax1.text(5, 7, 'Lm (互感)', ha='center', fontsize=10, color='purple', fontweight='bold')

# 近端和远端标注
ax1.text(0.5, 5.75, '近端', ha='right', fontsize=10, fontweight='bold')
ax1.text(9.5, 5.75, '远端', ha='left', fontsize=10, fontweight='bold')
ax1.text(0.5, 3.25, '近端', ha='right', fontsize=10, fontweight='bold')
ax1.text(9.5, 3.25, '远端', ha='left', fontsize=10, fontweight='bold')

ax1.legend(loc='upper center', ncol=2, fontsize=10)

# (b) 近端串扰 vs 互感/互容
ax2 = axes[0, 1]
Cm_fixed = 10e-12  # 10 pF/m
NEXT_Lm = [crosstalk_near_end(Cm_fixed, C, Lm, L, l, tr) for Lm in Lm_values]

Lm_fixed = 100e-9  # 100 nH/m
NEXT_Cm = [crosstalk_near_end(Cm, C, Lm_fixed, L, l, tr) for Cm in Cm_values]

ax2.plot(Lm_values*1e9, NEXT_Lm, 'b-', linewidth=2.5, label=f'Cm={Cm_fixed*1e12:.0f}pF/m')
ax2_twin = ax2.twinx()
ax2_twin.plot(Cm_values*1e12, NEXT_Cm, 'r--', linewidth=2.5, label=f'Lm={Lm_fixed*1e9:.0f}nH/m')

ax2.set_xlabel('互感 Lm (nH/m)', fontsize=11, color='blue')
ax2.set_ylabel('近端串扰', fontsize=11, color='blue')
ax2_twin.set_xlabel('互容 Cm (pF/m)', fontsize=11, color='red')
ax2_twin.set_ylabel('近端串扰', fontsize=11, color='red')
ax2.set_title('(b) 近端串扰 vs 耦合参数', fontsize=12, fontweight='bold')
ax2.tick_params(axis='x', colors='blue')
ax2.tick_params(axis='y', colors='blue')
ax2_twin.tick_params(axis='x', colors='red')
ax2_twin.tick_params(axis='y', colors='red')
ax2.grid(True, linestyle='--', alpha=0.5)

# (c) 远端串扰 vs 耦合长度
ax3 = axes[1, 0]
lengths = np.linspace(0.01, 0.5, 100)  # 1 cm 到 50 cm
Cm = 10e-12
Lm = 100e-9

FEXT = [crosstalk_far_end(Cm, C, Lm, L, l_val, tr, v) for l_val in lengths]

ax3.plot(lengths*100, FEXT, 'g-', linewidth=2.5)
ax3.axhline(y=0, color='black', linestyle='-', linewidth=1)
ax3.set_xlabel('耦合长度 (cm)', fontsize=11)
ax3.set_ylabel('远端串扰', fontsize=11)
ax3.set_title('(c) 远端串扰 vs 耦合长度', fontsize=12, fontweight='bold')
ax3.grid(True, linestyle='--', alpha=0.5)

# 添加注释
ax3.annotate('Cm/C > Lm/L\n正串扰', xy=(30, 0.02), fontsize=10, 
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
ax3.annotate('Cm/C < Lm/L\n负串扰', xy=(30, -0.02), fontsize=10,
            bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))

# (d) 串扰随频率变化
ax4 = axes[1, 1]
f_crosstalk = np.logspace(6, 9, 200)  # 1 MHz 到 1 GHz
tr_values = 1 / (3 * f_crosstalk)  # 上升时间与频率相关

NEXT_freq = []
FEXT_freq = []
for tr_val in tr_values:
    NEXT_freq.append(crosstalk_near_end(Cm, C, Lm, L, l, tr_val))
    FEXT_freq.append(crosstalk_far_end(Cm, C, Lm, L, l, tr_val, v))

ax4.semilogx(f_crosstalk/1e6, NEXT_freq, 'b-', linewidth=2.5, label='近端串扰')
ax4.semilogx(f_crosstalk/1e6, FEXT_freq, 'r--', linewidth=2.5, label='远端串扰')

ax4.set_xlabel('频率 (MHz)', fontsize=11)
ax4.set_ylabel('串扰系数', fontsize=11)
ax4.set_title('(d) 串扰 vs 频率', fontsize=12, fontweight='bold')
ax4.legend(loc='upper right', fontsize=10)
ax4.grid(True, which='both', linestyle='--', alpha=0.5)
ax4.set_xlim([1, 1000])

plt.tight_layout()
plt.savefig('crosstalk_analysis.png', dpi=150, bbox_inches='tight')
print("  ✓ crosstalk_analysis.png 已保存")
plt.close()

# ==============================================================================
# 4. 辐射发射估算
# ==============================================================================
print("\n[4/5] 辐射发射估算...")

def differential_mode_radiation(f, A, I, r):
    """
    差模辐射电场计算
    
    参数:
        f: 频率 (Hz)
        A: 环路面积 (m²)
        I: 电流 (A)
        r: 距离 (m)
    
    返回:
        电场强度 (V/m)
    """
    return 1.316e-14 * f**2 * A * I / r

def common_mode_radiation(f, L, I_cm, r):
    """
    共模辐射电场计算
    
    参数:
        f: 频率 (Hz)
        L: 电缆长度 (m)
        I_cm: 共模电流 (A)
        r: 距离 (m)
    
    返回:
        电场强度 (V/m)
    """
    return 6.28e-7 * f * L * I_cm / r

# 频率范围
f_rad = np.logspace(6, 9, 500)  # 1 MHz 到 1 GHz

# 参数
A = 10e-4  # 10 cm² 环路面积
I = 0.01   # 10 mA 电流
L_cable = 1.0  # 1 m 电缆长度
I_cm = 1e-6  # 1 μA 共模电流
r = 3  # 3 m 测量距离

# CISPR 22 Class B 限值 (3m距离)
f_limit = np.array([30, 88, 230, 1000]) * 1e6  # MHz to Hz
E_limit = np.array([30, 30, 37, 37])  # dBμV/m

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (a) 差模辐射 vs 环路面积
ax1 = axes[0, 0]
A_values = [1e-4, 5e-4, 10e-4, 20e-4]  # 1, 5, 10, 20 cm²
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']

for A_val, color in zip(A_values, colors):
    E_dm = differential_mode_radiation(f_rad, A_val, I, r)
    E_dm_db = 20 * np.log10(E_dm * 1e6)  # Convert to dBμV/m
    ax1.semilogx(f_rad/1e6, E_dm_db, linewidth=2, 
                 label=f'A={A_val*1e4:.0f}cm²', color=color)

# 绘制限值
ax1.semilogx(f_limit/1e6, E_limit, 'r--', linewidth=2.5, label='CISPR 22 Class B')
ax1.fill_between(f_limit/1e6, E_limit, 80, alpha=0.2, color='red')

ax1.set_xlabel('频率 (MHz)', fontsize=11)
ax1.set_ylabel('电场强度 (dBμV/m)', fontsize=11)
ax1.set_title('(a) 差模辐射 vs 环路面积', fontsize=12, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9)
ax1.grid(True, which='both', linestyle='--', alpha=0.5)
ax1.set_xlim([1, 1000])
ax1.set_ylim([0, 80])

# (b) 共模辐射 vs 电缆长度
ax2 = axes[0, 1]
L_values = [0.3, 0.5, 1.0, 2.0]  # 0.3, 0.5, 1, 2 m

for L_val, color in zip(L_values, colors):
    E_cm = common_mode_radiation(f_rad, L_val, I_cm, r)
    E_cm_db = 20 * np.log10(E_cm * 1e6)
    ax2.semilogx(f_rad/1e6, E_cm_db, linewidth=2,
                 label=f'L={L_val:.1f}m', color=color)

ax2.semilogx(f_limit/1e6, E_limit, 'r--', linewidth=2.5, label='CISPR 22 Class B')
ax2.fill_between(f_limit/1e6, E_limit, 80, alpha=0.2, color='red')

ax2.set_xlabel('频率 (MHz)', fontsize=11)
ax2.set_ylabel('电场强度 (dBμV/m)', fontsize=11)
ax2.set_title('(b) 共模辐射 vs 电缆长度', fontsize=12, fontweight='bold')
ax2.legend(loc='upper left', fontsize=9)
ax2.grid(True, which='both', linestyle='--', alpha=0.5)
ax2.set_xlim([1, 1000])
ax2.set_ylim([0, 80])

# (c) 差模 vs 共模辐射对比
ax3 = axes[1, 0]
E_dm = differential_mode_radiation(f_rad, A, I, r)
E_cm = common_mode_radiation(f_rad, L_cable, I_cm, r)

E_dm_db = 20 * np.log10(E_dm * 1e6)
E_cm_db = 20 * np.log10(E_cm * 1e6)

ax3.semilogx(f_rad/1e6, E_dm_db, 'b-', linewidth=2.5, label=f'差模 (A={A*1e4:.0f}cm²)')
ax3.semilogx(f_rad/1e6, E_cm_db, 'r--', linewidth=2.5, label=f'共模 (L={L_cable}m)')
ax3.semilogx(f_limit/1e6, E_limit, 'g:', linewidth=2.5, label='CISPR 22 Class B')

ax3.set_xlabel('频率 (MHz)', fontsize=11)
ax3.set_ylabel('电场强度 (dBμV/m)', fontsize=11)
ax3.set_title('(c) 差模与共模辐射对比', fontsize=12, fontweight='bold')
ax3.legend(loc='upper left', fontsize=10)
ax3.grid(True, which='both', linestyle='--', alpha=0.5)
ax3.set_xlim([1, 1000])
ax3.set_ylim([0, 80])

# (d) 辐射源示意图
ax4 = axes[1, 1]
ax4.set_xlim(0, 10)
ax4.set_ylim(0, 10)
ax4.axis('off')
ax4.set_title('(d) 差模与共模辐射机理', fontsize=12, fontweight='bold')

# 差模辐射 - 信号环路
ax4.add_patch(Rectangle((1, 6), 2, 1.5, fill=True, facecolor='#FFE4B5', 
                         edgecolor='black', linewidth=2))
ax4.plot([1, 3], [6.75, 6.75], 'b-', linewidth=3)
ax4.annotate('', xy=(3, 6.75), xytext=(1, 6.75),
            arrowprops=dict(arrowstyle='->', color='blue', lw=2))
ax4.text(2, 7.5, '信号环路', ha='center', fontsize=11, fontweight='bold')
ax4.text(2, 5.5, '差模辐射', ha='center', fontsize=10, color='blue')

# 绘制辐射波
for angle in np.linspace(0, 2*np.pi, 8, endpoint=False):
    x_start = 2 + 1.5 * np.cos(angle)
    y_start = 6.75 + 1.5 * np.sin(angle)
    x_end = 2 + 2.5 * np.cos(angle)
    y_end = 6.75 + 2.5 * np.sin(angle)
    ax4.annotate('', xy=(x_end, y_end), xytext=(x_start, y_start),
                arrowprops=dict(arrowstyle='->', color='blue', lw=1.5, 
                               connectionstyle='arc3,rad=0.1'))

# 共模辐射 - 电缆天线
ax4.plot([6, 8], [6.75, 6.75], 'r-', linewidth=4)
ax4.plot([6, 6], [6.75, 5.5], 'k-', linewidth=2)
ax4.plot([8, 8], [6.75, 5.5], 'k-', linewidth=2)
ax4.plot([6, 8], [5.5, 5.5], 'k-', linewidth=2)
ax4.text(7, 7.5, '电缆天线', ha='center', fontsize=11, fontweight='bold')
ax4.text(7, 5, '共模辐射', ha='center', fontsize=10, color='red')

# 绘制辐射波
for angle in np.linspace(0, 2*np.pi, 12, endpoint=False):
    x_start = 7 + 1.5 * np.cos(angle)
    y_start = 6.75 + 1.5 * np.sin(angle)
    x_end = 7 + 2.5 * np.cos(angle)
    y_end = 6.75 + 2.5 * np.sin(angle)
    ax4.annotate('', xy=(x_end, y_end), xytext=(x_start, y_start),
                arrowprops=dict(arrowstyle='->', color='red', lw=1.5,
                               connectionstyle='arc3,rad=0.1'))

# 添加说明
ax4.text(5, 2, '差模：电流在信号线中流动，形成小环路\n共模：电流在电缆和地之间流动，形成大天线',
        ha='center', fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.savefig('radiated_emission.png', dpi=150, bbox_inches='tight')
print("  ✓ radiated_emission.png 已保存")
plt.close()

# ==============================================================================
# 5. 接地阻抗分析
# ==============================================================================
print("\n[5/5] 接地阻抗分析...")

def ground_impedance(f, R_dc, L, w=1e-3, sigma=5.8e7):
    """
    计算地平面阻抗
    
    参数:
        f: 频率 (Hz)
        R_dc: 直流电阻 (Ω)
        L: 电感 (H)
        w: 导体宽度 (m)
        sigma: 电导率 (S/m)
    
    返回:
        阻抗 (Ω)
    """
    mu_0 = 4 * np.pi * 1e-7
    delta = skin_depth(f, sigma)
    
    # 高频电阻（考虑趋肤效应）
    if delta < w/2:
        R_ac = R_dc * w / (2 * delta)
    else:
        R_ac = R_dc
    
    Z = R_ac + 1j * 2 * np.pi * f * L
    return np.abs(Z)

# 频率范围
f_gnd = np.logspace(3, 9, 500)  # 1 kHz 到 1 GHz

# 地线参数
lengths_gnd = [0.01, 0.02, 0.05, 0.1]  # 1, 2, 5, 10 cm
R_dc_per_m = 0.01  # 10 mΩ/m
L_per_m = 1e-6     # 1 μH/m

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (a) 地线阻抗 vs 频率
ax1 = axes[0, 0]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']

for l_gnd, color in zip(lengths_gnd, colors):
    R_dc = R_dc_per_m * l_gnd
    L = L_per_m * l_gnd
    Z_gnd = [ground_impedance(f, R_dc, L) for f in f_gnd]
    ax1.loglog(f_gnd/1e6, Z_gnd, linewidth=2, 
               label=f'L={l_gnd*100:.0f}cm', color=color)

ax1.axhline(y=1, color='red', linestyle='--', linewidth=1.5, alpha=0.7, label='1Ω')
ax1.set_xlabel('频率 (MHz)', fontsize=11)
ax1.set_ylabel('阻抗 (Ω)', fontsize=11)
ax1.set_title('(a) 地线阻抗 vs 频率', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right', fontsize=9)
ax1.grid(True, which='both', linestyle='--', alpha=0.5)
ax1.set_xlim([0.001, 1000])

# (b) 单点接地 vs 多点接地
ax2 = axes[0, 1]
ax2.set_xlim(0, 10)
ax2.set_ylim(0, 10)
ax2.axis('off')
ax2.set_title('(b) 接地方式对比', fontsize=12, fontweight='bold')

# 单点接地
ax2.text(2.5, 9, '单点接地', ha='center', fontsize=12, fontweight='bold')
# 电路1
ax2.add_patch(FancyBboxPatch((1, 7), 1.5, 1, boxstyle="round,pad=0.05",
                              facecolor='#FFE4B5', edgecolor='black', linewidth=1.5))
ax2.text(1.75, 7.5, '电路1', ha='center', va='center', fontsize=9)
ax2.plot([1.75, 1.75], [7, 5.5], 'k-', linewidth=2)

# 电路2
ax2.add_patch(FancyBboxPatch((2.5, 7), 1.5, 1, boxstyle="round,pad=0.05",
                              facecolor='#E6E6FA', edgecolor='black', linewidth=1.5))
ax2.text(3.25, 7.5, '电路2', ha='center', va='center', fontsize=9)
ax2.plot([3.25, 3.25], [7, 5.5], 'k-', linewidth=2)

# 电路3
ax2.add_patch(FancyBboxPatch((4, 7), 1.5, 1, boxstyle="round,pad=0.05",
                              facecolor='#98FB98', edgecolor='black', linewidth=1.5))
ax2.text(4.75, 7.5, '电路3', ha='center', va='center', fontsize=9)
ax2.plot([4.75, 4.75], [7, 5.5], 'k-', linewidth=2)

# 公共接地点
ax2.plot([1.75, 4.75], [5.5, 5.5], 'k-', linewidth=2)
ax2.plot([3.25, 3.25], [5.5, 5], 'k-', linewidth=2)
ax2.plot([3, 3.5], [5, 5], 'k-', linewidth=3)
ax2.plot([3.1, 3.4], [4.8, 4.8], 'k-', linewidth=2)
ax2.plot([3.2, 3.3], [4.6, 4.6], 'k-', linewidth=1.5)
ax2.text(3.25, 4.3, '地', ha='center', fontsize=9)

# 多点接地
ax2.text(7.5, 9, '多点接地', ha='center', fontsize=12, fontweight='bold')
# 电路1
ax2.add_patch(FancyBboxPatch((6.5, 7), 1, 0.8, boxstyle="round,pad=0.05",
                              facecolor='#FFE4B5', edgecolor='black', linewidth=1.5))
ax2.text(7, 7.4, '电路1', ha='center', va='center', fontsize=8)
ax2.plot([7, 7], [7, 5.5], 'k-', linewidth=2)
ax2.plot([6.8, 7.2], [5.5, 5.5], 'k-', linewidth=3)
ax2.plot([6.9, 7.1], [5.3, 5.3], 'k-', linewidth=2)

# 电路2
ax2.add_patch(FancyBboxPatch((7.5, 7), 1, 0.8, boxstyle="round,pad=0.05",
                              facecolor='#E6E6FA', edgecolor='black', linewidth=1.5))
ax2.text(8, 7.4, '电路2', ha='center', va='center', fontsize=8)
ax2.plot([8, 8], [7, 5.5], 'k-', linewidth=2)
ax2.plot([7.8, 8.2], [5.5, 5.5], 'k-', linewidth=3)
ax2.plot([7.9, 8.1], [5.3, 5.3], 'k-', linewidth=2)

# 电路3
ax2.add_patch(FancyBboxPatch((8.5, 7), 1, 0.8, boxstyle="round,pad=0.05",
                              facecolor='#98FB98', edgecolor='black', linewidth=1.5))
ax2.text(9, 7.4, '电路3', ha='center', va='center', fontsize=8)
ax2.plot([9, 9], [7, 5.5], 'k-', linewidth=2)
ax2.plot([8.8, 9.2], [5.5, 5.5], 'k-', linewidth=3)
ax2.plot([8.9, 9.1], [5.3, 5.3], 'k-', linewidth=2)

# 地平面对象
ax2.add_patch(Rectangle((6.5, 4.8), 3, 0.3, facecolor='gray', alpha=0.5))
ax2.text(8, 4.5, '地平面', ha='center', fontsize=9)

# 说明
ax2.text(2.5, 3, '优点：无地环路\n缺点：高频阻抗大', ha='center', fontsize=9,
        bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))
ax2.text(7.5, 3, '优点：高频阻抗小\n缺点：可能形成地环路', ha='center', fontsize=9,
        bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.5))

# (c) 地平面阻抗分布
ax3 = axes[1, 0]
# 模拟地平面上的阻抗分布
x_gnd = np.linspace(0, 0.1, 100)  # 10 cm 地平面
y_gnd = np.linspace(0, 0.1, 100)
X_gnd, Y_gnd = np.meshgrid(x_gnd, y_gnd)

# 电流注入点在中心
I_inject_x, I_inject_y = 0.05, 0.05

# 计算阻抗分布（简化模型）
R_sheet = 0.001  # 1 mΩ/□
Z_distribution = R_sheet * np.sqrt((X_gnd - I_inject_x)**2 + (Y_gnd - I_inject_y)**2) / 0.01

im = ax3.contourf(X_gnd*100, Y_gnd*100, Z_distribution*1000, levels=20, cmap='hot')
ax3.plot(I_inject_x*100, I_inject_y*100, 'b+', markersize=15, mew=2, label='电流注入点')
ax3.set_xlabel('x (cm)', fontsize=11)
ax3.set_ylabel('y (cm)', fontsize=11)
ax3.set_title('(c) 地平面阻抗分布', fontsize=12, fontweight='bold')
ax3.legend(loc='upper right', fontsize=9)
plt.colorbar(im, ax=ax3, label='阻抗 (mΩ)')

# (d) 去耦电容效果
ax4 = axes[1, 1]
# 去耦电容的阻抗特性
C_decap = [0.1e-6, 1e-6, 10e-6]  # 0.1, 1, 10 μF
ESL = 2e-9  # 2 nH 等效串联电感
ESR = 0.1   # 0.1 Ω 等效串联电阻

for C_val, color in zip(C_decap, colors[:3]):
    Z_cap = np.abs(ESR + 1j*2*np.pi*f_gnd*ESL + 1/(1j*2*np.pi*f_gnd*C_val))
    ax4.loglog(f_gnd/1e6, Z_cap, linewidth=2, 
               label=f'C={C_val*1e6:.1f}μF', color=color)

# 目标阻抗
ax4.axhline(y=0.1, color='red', linestyle='--', linewidth=1.5, alpha=0.7, label='目标阻抗')
ax4.set_xlabel('频率 (MHz)', fontsize=11)
ax4.set_ylabel('阻抗 (Ω)', fontsize=11)
ax4.set_title('(d) 去耦电容阻抗特性', fontsize=12, fontweight='bold')
ax4.legend(loc='upper right', fontsize=9)
ax4.grid(True, which='both', linestyle='--', alpha=0.5)
ax4.set_xlim([0.001, 1000])

plt.tight_layout()
plt.savefig('grounding_analysis.png', dpi=150, bbox_inches='tight')
print("  ✓ grounding_analysis.png 已保存")
plt.close()

# ==============================================================================
# 完成
# ==============================================================================
print("\n" + "=" * 70)
print("所有仿真完成！")
print("=" * 70)

print("\n生成的仿真图：")
print("  1. shielding_effectiveness.png - 屏蔽效能分析")
print("  2. filter_design.png - EMI滤波器设计")
print("  3. crosstalk_analysis.png - 传输线串扰分析")
print("  4. radiated_emission.png - 辐射发射估算")
print("  5. grounding_analysis.png - 接地阻抗分析")

print("\n所有文件保存在: D:\\文档\\500仿真领域\\工程仿真\\高频电磁场仿真\\主题049")

7.1 仿真内容概述

1. 屏蔽效能分析：计算不同材料的屏蔽效能随频率的变化
2. 滤波器设计：EMI滤波器的插入损耗分析
3. 串扰仿真：传输线之间的近端和远端串扰
4. 辐射发射估算：差模和共模辐射发射计算
5. 接地阻抗分析：地平面阻抗随频率的变化

7.2 运行仿真

python run_simulation.py

仿真将生成以下图片：

• shielding_effectiveness.png - 屏蔽效能分析
• filter_design.png - 滤波器设计
• crosstalk_analysis.png - 串扰分析
• radiated_emission.png - 辐射发射估算
• grounding_analysis.png - 接地分析

7.3 生成GIF动图

python create_gif.py

将生成以下GIF动图：

• shielding_materials.gif - 屏蔽材料对比
• filter_response.gif - 滤波器频率响应
• crosstalk_variation.gif - 串扰变化
• radiation_pattern.gif - 辐射方向图
• impedance_frequency.gif - 阻抗频率特性