高频电磁场仿真-主题099-脑机接口电磁设计

kkchenjj

112人浏览 · 2026-04-05 20:57:30

kkchenjj · 2026-04-05 20:57:30 发布

主题099：脑机接口电磁设计

摘要

脑机接口（Brain-Computer Interface, BCI）是神经工程领域的前沿技术，通过直接读取和解码大脑神经信号，实现人脑与外部设备的直接通信。电磁设计在脑机接口系统中扮演着核心角色，涉及神经信号采集、无线数据传输、电磁安全性等关键技术。本主题系统阐述脑机接口的电磁理论基础，包括神经电信号的产生与传播机制、植入式天线设计、无线能量传输、电磁兼容性分析等内容。通过Python仿真，深入分析神经信号采集系统的性能、植入式天线的辐射特性、组织电磁安全性等关键问题，为脑机接口系统的设计与优化提供理论指导。

关键词

脑机接口；神经信号采集；植入式天线；电磁安全性；无线能量传输；神经工程；生物电磁学；电磁兼容

在这里插入图片描述

1. 引言

1.1 脑机接口技术概述

脑机接口（Brain-Computer Interface, BCI）是一种不依赖于外周神经和肌肉的直接通信通路，能够实现大脑与外部设备之间的信息交换。这一技术的核心在于通过传感器采集大脑的神经活动信号，经过信号处理和模式识别，将其转化为控制指令，驱动外部设备执行相应动作。

脑机接口技术的发展可以追溯到20世纪70年代，Vidal首次提出了"脑机接口"的概念。经过几十年的发展，脑机接口技术已经从实验室研究逐步走向临床应用。2010年代，随着微电子技术、信号处理技术和神经科学的进步，脑机接口技术取得了突破性进展。2020年代，Neuralink、Synchron等公司的介入，使得脑机接口技术进入了快速发展期。

根据信号采集方式的不同，脑机接口可分为以下几类：

非侵入式脑机接口：通过头皮表面采集脑电信号（EEG），具有安全性高、操作简便的优点，但信号空间分辨率较低，易受噪声干扰。

半侵入式脑机接口：将电极植入颅骨下方但不进入脑组织，如皮层脑电图（ECoG），在信号质量和安全性之间取得平衡。

侵入式脑机接口：将微电极阵列直接植入脑组织内部，能够获取高质量的神经信号，但存在手术风险和长期稳定性问题。

1.2 电磁设计在脑机接口中的重要性

电磁设计贯穿脑机接口系统的各个环节：

神经信号采集：微弱神经信号（微伏级）的精确采集需要高灵敏度、低噪声的电磁传感器设计。

无线数据传输：植入式设备需要将采集的神经信号无线传输到外部接收器，这要求高效的植入式天线设计。

无线能量传输：为避免频繁手术更换电池，植入式设备通常采用无线能量传输方式供电。

电磁安全性：植入式设备产生的电磁场必须控制在安全范围内，避免对脑组织造成热损伤或其他不良影响。

电磁兼容性：脑机接口系统需要在复杂的电磁环境中正常工作，避免与其他医疗设备或通信设备产生干扰。

1.3 脑机接口的应用前景

脑机接口技术具有广阔的应用前景：

医疗康复：帮助瘫痪患者控制假肢、轮椅等辅助设备，恢复部分运动功能；治疗帕金森病、癫痫等神经系统疾病。

神经科学研究：深入理解大脑的工作机制，探索意识的本质。

人机交互：开发更自然、更高效的人机交互方式，如意念控制、脑波游戏等。

增强人类能力：未来可能实现记忆增强、直接知识传输等功能。

2. 神经电信号的电磁特性

2.1 神经元的电生理基础

神经元是神经系统的基本功能单元，具有产生和传导电信号的能力。神经元的电活动主要基于离子通道的开关和离子跨膜流动。

静息电位：在静息状态下，神经元膜内外存在约-70mV的电位差，主要由钾离子外流维持。

动作电位：当神经元受到足够强的刺激时，钠离子通道开放，产生快速的膜电位去极化，形成动作电位。动作电位的幅度约为100mV，持续时间约1-2ms。

突触传递：神经元之间通过突触进行信号传递。当动作电位到达突触前膜时，引起神经递质释放，作用于突触后膜，产生兴奋性或抑制性突触后电位。

2.2 神经信号的电场分布

神经活动产生的电场可以在头皮表面被检测到。根据容积导体理论，脑内电流源在头皮表面产生的电位分布可以通过以下方式描述：

泊松方程：在准静态近似下，脑内电势分布满足泊松方程：

$∇⋅(σ∇ϕ)=−Iv\nabla \cdot (\sigma \nabla \phi) = -I_v$

其中， $σ\sigma$ 是电导率张量， $ϕ\phi$ 是电势， $I_v$ 是电流源密度。

多极子展开：对于远离电流源的区域，电势可以表示为多极子展开的形式：

$]\phi(\mathbf{r}) = \frac{1}{4\pi\sigma} \left[ \frac{Q}{r} + \frac{\mathbf{p} \cdot \mathbf{r}}{r^3} + \cdots \right]$

其中， $Q$ 是单极矩， $p\mathbf{p}$ 是偶极矩。

正问题与逆问题：

正问题：已知电流源分布，计算头皮表面电位分布
逆问题：已知头皮表面电位分布，反推电流源位置（病态问题，需要正则化方法）

2.3 神经信号的频谱特性

不同类型的神经信号具有不同的频率范围：

信号类型	频率范围	幅度范围	来源
局部场电位（LFP）	1-100 Hz	10-1000 μV	神经元群体突触活动
动作电位（Spike）	300-3000 Hz	50-500 μV	单个神经元放电
皮层脑电图（ECoG）	0.5-200 Hz	10-1000 μV	皮层表面电位
脑电图（EEG）	0.5-100 Hz	5-100 μV	头皮表面电位

脑电波的典型频段：

δ波：0.5-4 Hz，深度睡眠状态
θ波：4-8 Hz，浅睡眠或放松状态
α波：8-13 Hz，清醒放松状态
β波：13-30 Hz，清醒专注状态
γ波：30-100 Hz，认知活动

3. 神经信号采集系统设计

3.1 微电极阵列设计

微电极阵列是侵入式脑机接口的核心组件，用于直接记录神经元的电活动。

电极材料选择：

贵金属：铂、铱、金等，具有良好的生物相容性和电化学稳定性
导电聚合物：聚吡咯（PPy）、聚苯胺（PANI）等，具有柔性、低阻抗的特点
碳基材料：碳纳米管、石墨烯等，具有优异的电化学性能和生物相容性

电极几何设计：

电极尺寸：典型值为10-50 μm，需要在信号质量和组织损伤之间权衡
电极形状：圆形、方形、圆锥形等，影响电场分布和阻抗特性
电极间距：典型值为100-400 μm，影响空间分辨率和串扰

电极阻抗模型：
电极-电解质界面的阻抗可以用Randles等效电路描述：

$Ze=Rs+11Rct+jωCdlZ_e = R_s + \frac{1}{\frac{1}{R_{ct}} + j\omega C_{dl}}$

其中， $R_s$ 是溶液电阻， $R_{ct}$ 是电荷转移电阻， $C_{dl}$ 是双电层电容。

3.2 信号调理电路

神经信号采集需要高增益、低噪声的信号调理电路。

前置放大器：

输入阻抗：需要极高（>100 MΩ），避免信号衰减
噪声性能：等效输入噪声应低于5 μVrms
共模抑制比（CMRR）：>80 dB，抑制共模干扰

滤波器设计：

高通滤波：截止频率0.5-1 Hz，消除基线漂移
低通滤波：截止频率100-500 Hz（LFP）或3-10 kHz（Spike）
陷波滤波：50/60 Hz，消除工频干扰

模数转换器（ADC）：

分辨率：12-16位，满足动态范围要求
采样率：LFP 1-2 kHz，Spike 20-40 kHz
功耗：需要极低功耗设计（<100 μW/通道）

3.3 噪声分析与抑制

神经信号采集面临多种噪声源的干扰：

热噪声：由电阻和放大器产生，与带宽成正比
$Vthermal=4kTRBV_{thermal} = \sqrt{4kTRB}$

闪烁噪声（1/f噪声）：低频段占主导，与频率成反比

生物噪声：

肌电干扰（EMG）：来自头部肌肉活动
眼动伪迹（EOG）：来自眼球运动
心电干扰（ECG）：来自心脏活动

环境噪声：

工频干扰：50/60 Hz及其谐波
射频干扰：来自无线通信设备

噪声抑制技术：

屏蔽与接地设计
差分放大结构
自适应滤波算法
独立成分分析（ICA）

4. 植入式天线设计

4.1 生物组织电磁特性

人体组织是色散介质，其电磁特性随频率变化。对于脑组织，主要电磁参数包括：

相对介电常数：
$εr(f)=ε∞+Δε1+(2πfτ)1−α+σDCj2πfε0\varepsilon_r(f) = \varepsilon_{\infty} + \frac{\Delta\varepsilon}{1+(2\pi f\tau)^{1-\alpha}} + \frac{\sigma_{DC}}{j2\pi f\varepsilon_0}$

其中， $ε∞\varepsilon_{\infty}$ 是高频极限介电常数， $Δε\Delta\varepsilon$ 是介电增量， $τ\tau$ 是弛豫时间， $α\alpha$ 是Cole-Cole分布参数， $σDC\sigma_{DC}$ 是直流电导率。

电导率：
脑组织的电导率随频率增加而增大，典型值在0.1-1 S/m范围。

穿透深度：
电磁波在生物组织中的穿透深度：
$δ=1πfμσ\delta = \frac{1}{\sqrt{\pi f \mu \sigma}}$

在2.4 GHz频率下，脑组织的穿透深度约为1-2 cm。

4.2 植入式天线类型

平面倒F天线（PIFA）：

结构紧凑，适合植入式应用
通过短路针调节阻抗匹配
典型尺寸：10-20 mm

环形天线：

磁场耦合为主，受组织影响较小
适合近场能量传输
可实现圆极化辐射

螺旋天线：

具有较好的带宽特性
可实现小型化设计
适合多频段应用

贴片天线：

方向性较好
易于与电路集成
需要接地平面

4.3 天线性能优化

阻抗匹配：
植入式天线需要在组织环境中实现良好的阻抗匹配。常用技术包括：

短路针调节
开槽技术
加载电容/电感
匹配网络设计

效率提升：

使用高介电常数基板减小天线尺寸
采用磁电偶极子结构提高辐射效率
优化天线与组织的耦合

带宽扩展：

多谐振技术
耦合谐振器结构
缺陷地结构（DGS）

4.4 天线仿真方法

有限元法（FEM）：
使用ANSYS HFSS、COMSOL等软件，适合复杂几何结构的精确建模。

时域有限差分法（FDTD）：
使用XFdtd、SEMCAD等软件，适合宽带分析和瞬态仿真。

矩量法（MOM）：
使用FEKO、IE3D等软件，适合金属结构的快速分析。

等效电路法：
建立天线的集总参数模型，快速估算性能。

5. 无线能量传输

5.1 无线能量传输原理

植入式脑机接口设备通常采用无线能量传输方式供电，避免频繁手术更换电池。

电磁感应耦合：
基于法拉第电磁感应定律，通过初级线圈和次级线圈的磁耦合实现能量传输。

$Vsecondary=−MdIprimarydtV_{secondary} = -M \frac{dI_{primary}}{dt}$

其中， $M$ 是互感系数。

磁共振耦合：
当发射和接收线圈调谐到相同频率时，实现高效的能量传输。品质因数 $Q$ 越高，传输效率越高。

射频能量传输：
使用微波或毫米波频段，通过天线辐射实现能量传输，适合远距离应用。

5.2 能量传输系统效率

无线能量传输系统的总效率包括：

$ηtotal=ηPA×ηcoil×ηrectifier×ηregulator\eta_{total} = \eta_{PA} \times \eta_{coil} \times \eta_{rectifier} \times \eta_{regulator}$

其中：

$ηPA\eta_{PA}$ ：功率放大器效率（50-80%）
$ηcoil\eta_{coil}$ ：线圈耦合效率（30-70%）
$ηrectifier\eta_{rectifier}$ ：整流效率（60-85%）
$ηregulator\eta_{regulator}$ ：稳压器效率（80-95%）

典型系统的总效率为10-40%。

5.3 安全性考虑

无线能量传输需要严格控制电磁暴露水平：

比吸收率（SAR）：
单位质量组织吸收的电磁功率：
$\frac{\sigma |E|^2}{\rho}$

其中， $σ\sigma$ 是电导率， $E$ 是电场强度， $ρ\rho$ 是组织密度。

IEEE C95.1标准规定：

头部SAR限值：1.6 W/kg（1g平均）
全身SAR限值：0.08 W/kg

温升控制：
组织温升与SAR的关系：
$ΔT=SAR⋅tc\Delta T = \frac{SAR \cdot t}{c}$

其中， $c$ 是比热容， $t$ 是暴露时间。

6. 电磁安全性分析

6.1 生物电磁效应

电磁场与生物组织的相互作用包括：

热效应：
电磁能量被组织吸收转化为热能，导致温度升高。这是高强度电磁场的主要生物效应。

非热效应：
低强度电磁场可能引起的非热生物学效应，包括：

离子跨膜运输改变
酶活性变化
基因表达调控
神经兴奋性改变

刺激效应：
时变电磁场可在神经组织中感应电流，引起神经兴奋或抑制。

6.2 安全标准与限值

国际非电离辐射防护委员会（ICNIRP）：

职业暴露：电场强度 100 V/m，磁场强度 0.5 A/m
公众暴露：电场强度 40 V/m，磁场强度 0.2 A/m

IEEE C95.1标准：

频率范围：3 kHz - 300 GHz
基于SAR和功率密度的限值

医疗设备标准：

IEC 60601-1：医疗电气设备通用安全要求
ISO 14708：有源植入式医疗器械

6.3 植入式设备的电磁兼容性

植入式脑机接口需要在复杂的电磁环境中正常工作：

抗干扰设计：

屏蔽设计：金属屏蔽壳、导电涂层
滤波设计：电源滤波、信号滤波
接地设计：单点接地、浮地设计

电磁敏感性测试：

静电放电（ESD）测试
射频电磁场辐射抗扰度测试
电快速瞬变脉冲群（EFT）测试

与其他医疗设备的兼容性：

MRI兼容性：避免在强磁场中产生热损伤或图像伪影
心脏起搏器兼容性：避免相互干扰

7. Python仿真实现

7.1 仿真概述

本节通过Python实现脑机接口电磁设计的仿真分析，包括：

神经信号采集系统仿真
植入式天线性能分析
无线能量传输仿真
电磁安全性评估

7.2 仿真代码

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
主题099：脑机接口电磁设计 - 仿真程序
=========================================

本程序实现脑机接口电磁设计的仿真分析，包括：
1. 神经信号采集系统性能分析
2. 植入式天线辐射特性仿真
3. 无线能量传输效率分析
4. 电磁安全性评估
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle, FancyArrowPatch, Ellipse
from matplotlib.collections import PatchCollection
import matplotlib.patches as mpatches
import warnings
warnings.filterwarnings('ignore')

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

# 创建输出目录
import os
output_dir = r'd:\文档\500仿真领域\工程仿真\高频电磁场仿真\主题099'
os.makedirs(output_dir, exist_ok=True)


# =============================================================================
# 仿真1：神经信号采集系统分析
# =============================================================================
def simulation_1_neural_signal_acquisition():
    """神经信号采集系统性能仿真"""
    print("[1/5] 运行神经信号采集系统仿真...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 子图1: 神经元动作电位波形
    ax1 = axes[0, 0]
    
    t = np.linspace(0, 5, 1000)  # ms
    
    # 动作电位模型（基于Hodgkin-Huxley方程简化）
    def action_potential(t, t0):
        """生成动作电位波形"""
        mask = (t >= t0) & (t <= t0 + 2)
        v = np.zeros_like(t)
        if np.any(mask):
            t_local = t[mask] - t0
            v[mask] = -70 + 100 * (t_local/0.5) * np.exp(1 - t_local/0.5) * np.exp(-(t_local-0.5)**2/0.3)
        return v
    
    # 生成多个神经元的放电
    np.random.seed(42)
    v_total = -70 * np.ones_like(t)  # 静息电位
    
    for i in range(5):
        t_spike = np.random.uniform(0.5, 4.5)
        v_total += action_potential(t, t_spike)
    
    # 添加噪声
    noise = 2 * np.random.randn(len(t))
    v_noisy = v_total + noise
    
    ax1.plot(t, v_total, 'b-', linewidth=2, label='Clean Signal')
    ax1.plot(t, v_noisy, 'r-', alpha=0.5, linewidth=1, label='With Noise')
    ax1.axhline(y=-70, color='k', linestyle='--', alpha=0.5, label='Resting Potential')
    ax1.axhline(y=-55, color='g', linestyle='--', alpha=0.5, label='Threshold')
    
    ax1.set_xlabel('Time (ms)')
    ax1.set_ylabel('Membrane Potential (mV)')
    ax1.set_title('Neural Action Potential Waveform', fontsize=11)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_ylim(-90, 50)
    
    # 子图2: 脑电波频谱分析
    ax2 = axes[0, 1]
    
    f = np.linspace(0.5, 100, 500)  # Hz
    
    # 模拟不同频段的功率谱
    # δ波: 0.5-4 Hz
    delta = 50 * np.exp(-(f-2)**2/2)
    # θ波: 4-8 Hz
    theta = 30 * np.exp(-(f-6)**2/2)
    # α波: 8-13 Hz
    alpha = 60 * np.exp(-(f-10)**2/3)
    # β波: 13-30 Hz
    beta = 25 * np.exp(-(f-20)**2/20)
    # γ波: 30-100 Hz
    gamma = 15 * np.exp(-(f-50)**2/200)
    
    psd = delta + theta + alpha + beta + gamma + 5
    
    ax2.semilogy(f, psd, 'b-', linewidth=2)
    
    # 标记频段
    ax2.axvspan(0.5, 4, alpha=0.2, color='red', label='δ (0.5-4 Hz)')
    ax2.axvspan(4, 8, alpha=0.2, color='orange', label='θ (4-8 Hz)')
    ax2.axvspan(8, 13, alpha=0.2, color='yellow', label='α (8-13 Hz)')
    ax2.axvspan(13, 30, alpha=0.2, color='green', label='β (13-30 Hz)')
    ax2.axvspan(30, 100, alpha=0.2, color='blue', label='γ (30-100 Hz)')
    
    ax2.set_xlabel('Frequency (Hz)')
    ax2.set_ylabel('Power Spectral Density (μV²/Hz)')
    ax2.set_title('EEG Frequency Bands', fontsize=11)
    ax2.legend(loc='upper right', fontsize=8)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 100)
    
    # 子图3: 电极阻抗特性
    ax3 = axes[0, 2]
    
    f_imp = np.logspace(0, 4, 100)  # 1 Hz - 10 kHz
    
    # Randles等效电路模型
    Rs = 100  # 溶液电阻
    Rct = 10000  # 电荷转移电阻
    Cdl = 1e-6  # 双电层电容
    
    Zc = 1 / (1j * 2 * np.pi * f_imp * Cdl)
    Z_parallel = (Rct * Zc) / (Rct + Zc)
    Z_total = Rs + Z_parallel
    
    Z_mag = np.abs(Z_total)
    Z_phase = np.angle(Z_total, deg=True)
    
    ax3_twin = ax3.twinx()
    
    line1 = ax3.loglog(f_imp, Z_mag, 'b-', linewidth=2, label='|Z|')
    line2 = ax3_twin.semilogx(f_imp, Z_phase, 'r--', linewidth=2, label='Phase')
    
    ax3.set_xlabel('Frequency (Hz)')
    ax3.set_ylabel('Impedance Magnitude (Ω)', color='b')
    ax3_twin.set_ylabel('Phase (deg)', color='r')
    ax3.tick_params(axis='y', labelcolor='b')
    ax3_twin.tick_params(axis='y', labelcolor='r')
    
    ax3.set_title('Electrode Impedance Characteristics', fontsize=11)
    ax3.grid(True, alpha=0.3, which='both')
    
    lines = line1 + line2
    labels = [l.get_label() for l in lines]
    ax3.legend(lines, labels, loc='upper right')
    
    # 子图4: 信号噪声比分析
    ax4 = axes[1, 0]
    
    # 不同噪声水平下的SNR
    signal_power = 100  # μV²
    noise_levels = np.linspace(1, 50, 50)  # μVrms
    
    snr = 20 * np.log10(np.sqrt(signal_power) / noise_levels)
    
    ax4.plot(noise_levels, snr, 'b-', linewidth=2)
    ax4.axhline(y=20, color='r', linestyle='--', label='Good SNR (20 dB)')
    ax4.axhline(y=10, color='orange', linestyle='--', label='Acceptable SNR (10 dB)')
    
    # 标记工作点
    noise_design = 10
    snr_design = 20 * np.log10(np.sqrt(signal_power) / noise_design)
    ax4.plot(noise_design, snr_design, 'go', markersize=10)
    ax4.text(noise_design + 2, snr_design, f'Design Point\n{snr_design:.1f} dB', fontsize=9)
    
    ax4.set_xlabel('Noise Level (μVrms)')
    ax4.set_ylabel('SNR (dB)')
    ax4.set_title('Signal-to-Noise Ratio Analysis', fontsize=11)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    ax4.set_xlim(0, 50)
    ax4.set_ylim(0, 50)
    
    # 子图5: 多通道采集系统
    ax5 = axes[1, 1]
    
    # 绘制电极阵列布局
    n_channels = 16
    grid_size = 4
    
    for i in range(grid_size):
        for j in range(grid_size):
            x = j * 2
            y = i * 2
            circle = Circle((x, y), 0.3, facecolor='gold', edgecolor='black', linewidth=1)
            ax5.add_patch(circle)
            ax5.text(x, y, f'{i*grid_size+j+1}', ha='center', va='center', fontsize=8)
    
    # 绘制连接线
    for i in range(grid_size):
        for j in range(grid_size-1):
            ax5.plot([j*2, (j+1)*2], [i*2, i*2], 'b-', linewidth=1, alpha=0.5)
    for i in range(grid_size-1):
        for j in range(grid_size):
            ax5.plot([j*2, j*2], [i*2, (i+1)*2], 'b-', linewidth=1, alpha=0.5)
    
    ax5.set_xlim(-1, 7)
    ax5.set_ylim(-1, 7)
    ax5.set_aspect('equal')
    ax5.set_xlabel('x (mm)')
    ax5.set_ylabel('y (mm)')
    ax5.set_title('Multi-Channel Electrode Array (4x4)', fontsize=11)
    ax5.grid(True, alpha=0.3)
    
    # 子图6: 信号采集链路效率
    ax6 = axes[1, 2]
    
    stages = ['Electrode', 'Preamp', 'Filter', 'ADC', 'DSP']
    gains = [0, 40, -3, 0, 0]  # dB
    noise_figures = [0, 3, 1, 2, 0.5]  # dB
    
    x_pos = np.arange(len(stages))
    width = 0.35
    
    bars1 = ax6.bar(x_pos - width/2, gains, width, label='Gain (dB)', color='steelblue')
    bars2 = ax6.bar(x_pos + width/2, noise_figures, width, label='Noise Figure (dB)', color='coral')
    
    ax6.set_xticks(x_pos)
    ax6.set_xticklabels(stages)
    ax6.set_ylabel('Value (dB)')
    ax6.set_title('Signal Acquisition Chain Analysis', fontsize=11)
    ax6.legend()
    ax6.grid(True, alpha=0.3, axis='y')
    ax6.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig1_neural_signal_acquisition.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ Figure 1 saved: Neural Signal Acquisition System")


# =============================================================================
# 仿真2：植入式天线设计
# =============================================================================
def simulation_2_implantable_antenna():
    """植入式天线性能仿真"""
    print("[2/5] 运行植入式天线仿真...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 子图1: 天线结构示意图
    ax1 = axes[0, 0]
    
    # 绘制PIFA天线结构
    # 辐射贴片
    patch = Rectangle((1, 1), 8, 6, facecolor='gold', edgecolor='black', linewidth=2)
    ax1.add_patch(patch)
    
    # 短路针
    ax1.plot([2, 2], [1, 0.2], 'k-', linewidth=3)
    ax1.add_patch(Circle((2, 0.2), 0.15, facecolor='silver', edgecolor='black'))
    
    # 馈电点
    ax1.plot([7, 7], [1, 0.2], 'k-', linewidth=3)
    ax1.add_patch(Circle((7, 0.2), 0.15, facecolor='red', edgecolor='black'))
    
    # 接地面
    ground = Rectangle((0, -0.5), 10, 0.5, facecolor='silver', edgecolor='black', linewidth=2)
    ax1.add_patch(ground)
    
    # 基板
    substrate = Rectangle((0, 0), 10, 8, facecolor='lightblue', edgecolor='black', 
                          linewidth=1, alpha=0.3, linestyle='--')
    ax1.add_patch(substrate)
    
    # 标注
    ax1.text(5, 4, 'Patch', ha='center', fontsize=10)
    ax1.text(2, -1, 'Short', ha='center', fontsize=9)
    ax1.text(7, -1, 'Feed', ha='center', fontsize=9, color='red')
    ax1.text(5, 8.5, 'PIFA Antenna Structure', ha='center', fontsize=11, fontweight='bold')
    
    ax1.set_xlim(-1, 11)
    ax1.set_ylim(-2, 10)
    ax1.set_aspect('equal')
    ax1.axis('off')
    
    # 子图2: 生物组织介电特性
    ax2 = axes[0, 1]
    
    f_tissue = np.logspace(6, 10, 100)  # 1 MHz - 10 GHz
    
    # 脑组织的介电常数和电导率（简化模型）
    # 基于Cole-Cole模型
    eps_inf = 4
    eps_s = 80
    tau = 10e-12  # s
    sigma_dc = 0.5  # S/m
    
    eps_tissue = eps_inf + (eps_s - eps_inf) / (1 + (2*np.pi*f_tissue*tau)**2)
    sigma_tissue = sigma_dc + 2*np.pi*f_tissue*8.854e-12*(eps_s - eps_inf)*(2*np.pi*f_tissue*tau)**2/(1+(2*np.pi*f_tissue*tau)**2)
    
    ax2_twin = ax2.twinx()
    
    line1 = ax2.loglog(f_tissue/1e6, eps_tissue, 'b-', linewidth=2, label='εr')
    line2 = ax2_twin.loglog(f_tissue/1e6, sigma_tissue, 'r--', linewidth=2, label='σ')
    
    ax2.axvline(x=402, color='g', linestyle=':', alpha=0.7, label='MICS 402 MHz')
    ax2.axvline(x=2450, color='orange', linestyle=':', alpha=0.7, label='ISM 2.45 GHz')
    
    ax2.set_xlabel('Frequency (MHz)')
    ax2.set_ylabel('Relative Permittivity εr', color='b')
    ax2_twin.set_ylabel('Conductivity σ (S/m)', color='r')
    ax2.tick_params(axis='y', labelcolor='b')
    ax2_twin.tick_params(axis='y', labelcolor='r')
    
    ax2.set_title('Brain Tissue Electromagnetic Properties', fontsize=11)
    ax2.grid(True, alpha=0.3, which='both')
    
    lines = line1 + line2
    labels = [l.get_label() for l in lines]
    ax2.legend(lines, labels, loc='upper right')
    
    # 子图3: 天线回波损耗
    ax3 = axes[0, 2]
    
    f_ant = np.linspace(400, 450, 100) * 1e6  # MICS频段
    
    # 模拟天线在组织中的回波损耗
    f0 = 402e6  # 谐振频率
    Q = 20  # 品质因数
    BW = f0 / Q
    
    S11 = -10 * np.log10(1 + ((f_ant - f0)/BW)**2)
    S11 = np.clip(S11, -30, 0)
    
    ax3.plot(f_ant/1e6, S11, 'b-', linewidth=2)
    ax3.axhline(y=-10, color='r', linestyle='--', label='-10 dB Bandwidth')
    ax3.fill_between(f_ant/1e6, -30, S11, alpha=0.3)
    
    # 标记带宽
    f_lower = f0/1e6 - BW/2e6
    f_upper = f0/1e6 + BW/2e6
    ax3.axvline(x=f_lower, color='g', linestyle=':', alpha=0.5)
    ax3.axvline(x=f_upper, color='g', linestyle=':', alpha=0.5)
    ax3.text(f0/1e6, -25, f'BW = {BW/1e6:.1f} MHz', ha='center', fontsize=9)
    
    ax3.set_xlabel('Frequency (MHz)')
    ax3.set_ylabel('S11 (dB)')
    ax3.set_title('Antenna Return Loss (in Tissue)', fontsize=11)
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_ylim(-30, 0)
    
    # 子图4: 天线辐射方向图
    ax4 = axes[1, 0]
    
    theta = np.linspace(0, 2*np.pi, 360)
    
    # 模拟在组织中的辐射方向图（全向性降低）
    pattern = 1 - 0.3 * np.cos(theta)**2  # 稍微偏离全向
    pattern_dB = 10 * np.log10(pattern)
    
    ax4 = plt.subplot(2, 3, 4, projection='polar')
    ax4.plot(theta, pattern_dB + 30, linewidth=2)
    ax4.fill(theta, pattern_dB + 30, alpha=0.3)
    ax4.set_ylim(0, 30)
    ax4.set_title('Radiation Pattern (E-plane)', fontsize=11, pad=20)
    ax4.set_theta_zero_location('N')
    ax4.set_theta_direction(-1)
    
    # 子图5: 天线增益与效率
    ax5 = axes[1, 1]
    
    frequencies = np.linspace(400, 450, 50)  # MHz
    
    # 模拟增益和效率随频率变化
    gain = -15 - 0.1 * (frequencies - 402)**2  # dBi
    efficiency = 0.1 * np.exp(-(frequencies - 402)**2/100)  # %
    
    ax5_twin = ax5.twinx()
    
    line1 = ax5.plot(frequencies, gain, 'b-', linewidth=2, label='Gain')
    line2 = ax5_twin.plot(frequencies, efficiency * 100, 'r--', linewidth=2, label='Efficiency')
    
    ax5.axvline(x=402, color='g', linestyle=':', alpha=0.7)
    
    ax5.set_xlabel('Frequency (MHz)')
    ax5.set_ylabel('Gain (dBi)', color='b')
    ax5_twin.set_ylabel('Efficiency (%)', color='r')
    ax5.tick_params(axis='y', labelcolor='b')
    ax5_twin.tick_params(axis='y', labelcolor='r')
    
    ax5.set_title('Antenna Gain and Efficiency', fontsize=11)
    ax5.grid(True, alpha=0.3)
    
    lines = line1 + line2
    labels = [l.get_label() for l in lines]
    ax5.legend(lines, labels, loc='lower left')
    
    # 子图6: 不同天线类型对比
    ax6 = axes[1, 2]
    
    antenna_types = ['PIFA', 'Loop', 'Spiral', 'Dipole', 'Patch']
    sizes = [10, 15, 12, 20, 18]  # mm
    gains = [-12, -15, -14, -18, -10]  # dBi
    bandwidths = [15, 8, 25, 5, 10]  # MHz
    
    x = np.arange(len(antenna_types))
    width = 0.25
    
    # 归一化显示
    size_norm = np.array(sizes) / max(sizes) * 100
    gain_norm = (np.array(gains) + 20) / 10 * 100
    bw_norm = np.array(bandwidths) / max(bandwidths) * 100
    
    ax6.bar(x - width, size_norm, width, label='Size (normalized)', color='skyblue')
    ax6.bar(x, gain_norm, width, label='Gain (normalized)', color='lightgreen')
    ax6.bar(x + width, bw_norm, width, label='Bandwidth (normalized)', color='salmon')
    
    ax6.set_xticks(x)
    ax6.set_xticklabels(antenna_types)
    ax6.set_ylabel('Normalized Value (%)')
    ax6.set_title('Antenna Type Comparison', fontsize=11)
    ax6.legend()
    ax6.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig2_implantable_antenna.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ Figure 2 saved: Implantable Antenna Design")


# =============================================================================
# 仿真3：无线能量传输
# =============================================================================
def simulation_3_wireless_power_transfer():
    """无线能量传输仿真"""
    print("[3/5] 运行无线能量传输仿真...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 子图1: 能量传输系统示意图
    ax1 = axes[0, 0]
    
    # 外部发射线圈
    tx_coil = Circle((2, 5), 1.5, fill=False, edgecolor='blue', linewidth=3)
    ax1.add_patch(tx_coil)
    ax1.text(2, 5, 'Tx', ha='center', va='center', fontsize=12, color='blue')
    ax1.text(2, 3, 'External Coil', ha='center', fontsize=10)
    
    # 植入接收线圈
    rx_coil = Circle((8, 5), 1, fill=False, edgecolor='red', linewidth=3)
    ax1.add_patch(rx_coil)
    ax1.text(8, 5, 'Rx', ha='center', va='center', fontsize=12, color='red')
    ax1.text(8, 3, 'Implanted Coil', ha='center', fontsize=10)
    
    # 皮肤/组织层
    ax1.axhspan(4, 6, xmin=0.4, xmax=0.6, alpha=0.3, color='pink')
    ax1.text(5, 6.5, 'Skin/Tissue', ha='center', fontsize=10)
    
    # 磁场线
    for i in range(5):
        y_offset = 4 + i * 0.5
        ax1.annotate('', xy=(7, y_offset), xytext=(3, y_offset),
                    arrowprops=dict(arrowstyle='->', color='green', lw=2, 
                                   connectionstyle='arc3,rad=0.1'))
    
    ax1.set_xlim(0, 10)
    ax1.set_ylim(0, 10)
    ax1.set_aspect('equal')
    ax1.set_title('Wireless Power Transfer System', fontsize=11)
    ax1.axis('off')
    
    # 子图2: 耦合系数与距离关系
    ax2 = axes[0, 1]
    
    distance = np.linspace(5, 50, 100)  # mm
    
    # 耦合系数随距离变化（近似模型）
    k = 0.5 * (1 / (1 + (distance/15)**3))
    
    ax2.plot(distance, k, 'b-', linewidth=2)
    ax2.fill_between(distance, 0, k, alpha=0.3)
    
    # 标记工作点
    d_design = 20
    k_design = 0.5 * (1 / (1 + (d_design/15)**3))
    ax2.plot(d_design, k_design, 'ro', markersize=10)
    ax2.text(d_design + 2, k_design, f'k = {k_design:.3f}', fontsize=9)
    
    ax2.set_xlabel('Distance (mm)')
    ax2.set_ylabel('Coupling Coefficient k')
    ax2.set_title('Coupling vs Distance', fontsize=11)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(5, 50)
    ax2.set_ylim(0, 0.6)
    
    # 子图3: 传输效率分析
    ax3 = axes[0, 2]
    
    k_range = np.linspace(0.01, 0.5, 100)
    Q = 50  # 线圈品质因数
    
    # 最大效率公式
    eta_max = k_range**2 * Q**2 / (1 + np.sqrt(1 + k_range**2 * Q**2))**2
    
    ax3.plot(k_range, eta_max * 100, 'b-', linewidth=2)
    ax3.fill_between(k_range, 0, eta_max * 100, alpha=0.3)
    
    # 标记不同Q值
    for Q_val in [30, 50, 100]:
        eta_q = k_range**2 * Q_val**2 / (1 + np.sqrt(1 + k_range**2 * Q_val**2))**2
        ax3.plot(k_range, eta_q * 100, '--', alpha=0.5, label=f'Q={Q_val}')
    
    ax3.set_xlabel('Coupling Coefficient k')
    ax3.set_ylabel('Maximum Efficiency (%)')
    ax3.set_title('Power Transfer Efficiency', fontsize=11)
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_xlim(0, 0.5)
    ax3.set_ylim(0, 100)
    
    # 子图4: 整流电路效率
    ax4 = axes[1, 0]
    
    P_in = np.linspace(-20, 20, 100)  # dBm
    P_in_linear = 10**((P_in - 30)/10)  # W
    
    # 整流效率模型
    eta_rect = 0.8 * (1 - np.exp(-P_in_linear * 1000)) * (1 - 0.05 * P_in_linear)
    eta_rect = np.clip(eta_rect, 0, 0.85)
    
    ax4.plot(P_in, eta_rect * 100, 'b-', linewidth=2)
    ax4.axhline(y=70, color='r', linestyle='--', label='Target Efficiency 70%')
    ax4.axvline(x=0, color='g', linestyle='--', label='Design Point 0 dBm')
    
    ax4.set_xlabel('Input Power (dBm)')
    ax4.set_ylabel('Rectification Efficiency (%)')
    ax4.set_title('Rectifier Efficiency Curve', fontsize=11)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    ax4.set_xlim(-20, 20)
    ax4.set_ylim(0, 100)
    
    # 子图5: 系统功率预算
    ax5 = axes[1, 1]
    
    stages = ['Tx Amp', 'Coil', 'Rectifier', 'Regulator', 'Load']
    power_levels = [100, 70, 50, 45, 40]  # mW
    losses = [0, 30, 20, 5, 5]  # %
    
    x_pos = np.arange(len(stages))
    
    bars = ax5.bar(x_pos, power_levels, color=['blue', 'cyan', 'green', 'yellow', 'red'],
                   alpha=0.7, edgecolor='black')
    
    # 添加数值标签
    for bar, power in zip(bars, power_levels):
        height = bar.get_height()
        ax5.text(bar.get_x() + bar.get_width()/2., height + 1,
                f'{power} mW', ha='center', va='bottom', fontsize=9)
    
    # 添加效率箭头
    for i in range(len(stages) - 1):
        ax5.annotate('', xy=(i+1, power_levels[i+1]), xytext=(i+0.4, power_levels[i]),
                    arrowprops=dict(arrowstyle='->', color='red', lw=2))
        ax5.text(i+0.7, (power_levels[i] + power_levels[i+1])/2, 
                f'{losses[i+1]}%', fontsize=8, color='red')
    
    ax5.set_xticks(x_pos)
    ax5.set_xticklabels(stages)
    ax5.set_ylabel('Power Level (mW)')
    ax5.set_title('System Power Budget', fontsize=11)
    ax5.grid(True, alpha=0.3, axis='y')
    ax5.set_ylim(0, 120)
    
    # 子图6: 频率选择分析
    ax6 = axes[1, 2]
    
    frequencies = np.array([13.56, 402, 868, 915, 2450]) * 1e6  # Hz
    labels = ['HF\n13.56', 'MICS\n402', 'ISM-EU\n868', 'ISM-US\n915', 'ISM-2.4G\n2450']
    
    # 各频段的特性评分（1-10）
    penetration = [3, 9, 7, 7, 5]  # 穿透能力
    data_rate = [2, 4, 6, 6, 9]  # 数据速率
    efficiency = [8, 6, 7, 7, 5]  # 能量传输效率
    
    x = np.arange(len(labels))
    width = 0.25
    
    ax6.bar(x - width, penetration, width, label='Penetration', color='skyblue')
    ax6.bar(x, data_rate, width, label='Data Rate', color='lightgreen')
    ax6.bar(x + width, efficiency, width, label='Efficiency', color='salmon')
    
    ax6.set_xticks(x)
    ax6.set_xticklabels(labels, fontsize=9)
    ax6.set_ylabel('Score (1-10)')
    ax6.set_title('Frequency Band Comparison', fontsize=11)
    ax6.legend()
    ax6.grid(True, alpha=0.3, axis='y')
    ax6.set_ylim(0, 10)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig3_wireless_power_transfer.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ Figure 3 saved: Wireless Power Transfer")


# =============================================================================
# 仿真4：电磁安全性评估
# =============================================================================
def simulation_4_em_safety():
    """电磁安全性评估仿真"""
    print("[4/5] 运行电磁安全性评估仿真...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 子图1: SAR分布图
    ax1 = axes[0, 0]
    
    # 创建简化的头部模型
    head_radius = 8  # cm
    
    x = np.linspace(-10, 10, 200)
    y = np.linspace(-10, 10, 200)
    X, Y = np.meshgrid(x, y)
    R = np.sqrt(X**2 + Y**2)
    
    # 模拟SAR分布（高斯分布）
    implant_depth = 3  # cm
    sar_max = 1.5  # W/kg
    sar = sar_max * np.exp(-((R - implant_depth)**2) / 2)
    sar[R > head_radius] = 0
    
    levels = [0.1, 0.5, 1.0, 1.5, 2.0]
    cs = ax1.contourf(X, Y, sar, levels=levels, cmap='hot')
    plt.colorbar(cs, ax=ax1, label='SAR (W/kg)')
    
    # 绘制头部轮廓
    head = Circle((0, 0), head_radius, fill=False, edgecolor='blue', linewidth=2)
    ax1.add_patch(head)
    
    # 标记植入位置
    ax1.plot(0, implant_depth, 'b+', markersize=15, mew=2)
    ax1.text(1, implant_depth + 1, 'Implant', fontsize=9)
    
    # 安全限值线
    safe_limit = 1.6
    ax1.contour(X, Y, sar, levels=[safe_limit], colors='green', linewidths=2)
    
    ax1.set_xlabel('x (cm)')
    ax1.set_ylabel('y (cm)')
    ax1.set_title('SAR Distribution (Cross Section)', fontsize=11)
    ax1.set_aspect('equal')
    
    # 子图2: 温度升高分布
    ax2 = axes[0, 1]
    
    # 基于SAR计算温升（简化模型）
    rho = 1000  # kg/m³
    c = 3600  # J/(kg·K)
    t_exposure = 3600  # s (1 hour)
    
    delta_T = sar * t_exposure / (rho * c) * 1000  # °C
    
    levels_T = [0.1, 0.5, 1.0, 2.0, 3.0]
    cs2 = ax2.contourf(X, Y, delta_T, levels=levels_T, cmap='coolwarm')
    plt.colorbar(cs2, ax=ax2, label='ΔT (°C)')
    
    head2 = Circle((0, 0), head_radius, fill=False, edgecolor='blue', linewidth=2)
    ax2.add_patch(head2)
    
    ax2.set_xlabel('x (cm)')
    ax2.set_ylabel('y (cm)')
    ax2.set_title('Temperature Rise Distribution', fontsize=11)
    ax2.set_aspect('equal')
    
    # 子图3: 安全限值对比
    ax3 = axes[0, 2]
    
    standards = ['ICNIRP\nGeneral', 'ICNIRP\nOccupational', 'IEEE\nC95.1', 'FCC\nPart 18']
    sar_limits = [0.08, 0.4, 0.08, 0.4]  # W/kg (全身)
    sar_local = [2.0, 10.0, 1.6, 8.0]  # W/kg (局部)
    
    x = np.arange(len(standards))
    width = 0.35
    
    bars1 = ax3.bar(x - width/2, sar_limits, width, label='Whole Body', color='lightblue')
    bars2 = ax3.bar(x + width/2, sar_local, width, label='Local (1g/10g)', color='coral')
    
    # 标记限值线
    ax3.axhline(y=0.08, color='r', linestyle='--', alpha=0.5, label='Most Stringent')
    
    ax3.set_xticks(x)
    ax3.set_xticklabels(standards, fontsize=9)
    ax3.set_ylabel('SAR Limit (W/kg)')
    ax3.set_title('EM Safety Standards Comparison', fontsize=11)
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='y')
    ax3.set_ylim(0, 12)
    
    # 子图4: 暴露时间与SAR限值
    ax4 = axes[1, 0]
    
    time = np.linspace(1, 3600, 100)  # s
    
    # 短时暴露允许更高的SAR（简化模型）
    sar_short = 0.08 * np.sqrt(3600 / time)
    sar_short = np.clip(sar_short, 0.08, 2.0)
    
    ax4.loglog(time, sar_short, 'b-', linewidth=2, label='Allowable SAR')
    ax4.axhline(y=0.08, color='r', linestyle='--', label='Baseline Limit')
    ax4.fill_between(time, 0.08, sar_short, alpha=0.3)
    
    ax4.set_xlabel('Exposure Time (s)')
    ax4.set_ylabel('SAR Limit (W/kg)')
    ax4.set_title('SAR Limit vs Exposure Time', fontsize=11)
    ax4.legend()
    ax4.grid(True, alpha=0.3, which='both')
    
    # 子图5: 不同组织的电磁特性
    ax5 = axes[1, 1]
    
    tissues = ['Skin', 'Fat', 'Muscle', 'Bone', 'Brain', 'CSF']
    epsilon_r = [40, 5, 55, 15, 50, 70]  # 相对介电常数 @ 402 MHz
    conductivity = [0.6, 0.04, 0.8, 0.1, 0.5, 2.0]  # S/m
    
    x = np.arange(len(tissues))
    width = 0.35
    
    ax5_twin = ax5.twinx()
    
    bars1 = ax5.bar(x - width/2, epsilon_r, width, label='εr', color='skyblue')
    bars2 = ax5_twin.bar(x + width/2, conductivity, width, label='σ', color='coral')
    
    ax5.set_xticks(x)
    ax5.set_xticklabels(tissues, rotation=45, ha='right')
    ax5.set_ylabel('Relative Permittivity εr', color='b')
    ax5_twin.set_ylabel('Conductivity σ (S/m)', color='r')
    ax5.tick_params(axis='y', labelcolor='b')
    ax5_twin.tick_params(axis='y', labelcolor='r')
    
    ax5.set_title('Tissue Properties @ 402 MHz', fontsize=11)
    ax5.grid(True, alpha=0.3, axis='y')
    
    # 子图6: 风险评估矩阵
    ax6 = axes[1, 2]
    
    # 创建风险评估矩阵
    risk_factors = ['SAR Level', 'Exposure Time', 'Frequency', 'Tissue Type', 'Device Power']
    risk_scores = [3, 2, 2, 4, 3]  # 1-5风险等级
    mitigation = [4, 5, 3, 2, 4]  # 1-5缓解措施有效性
    
    colors = ['green', 'yellow', 'orange', 'red', 'darkred']
    
    y_pos = np.arange(len(risk_factors))
    
    # 绘制风险等级
    bars = ax6.barh(y_pos, risk_scores, color=[colors[s-1] for s in risk_scores], alpha=0.7)
    
    # 添加缓解措施标记
    for i, (score, mit) in enumerate(zip(risk_scores, mitigation)):
        ax6.plot(score + 0.1, i, 'bo', markersize=8)
        ax6.plot([score + 0.1, score + 0.1 + mit*0.3], [i, i], 'b-', linewidth=2)
        ax6.plot(score + 0.1 + mit*0.3, i, 'b>', markersize=8)
    
    ax6.set_yticks(y_pos)
    ax6.set_yticklabels(risk_factors)
    ax6.set_xlabel('Risk Level (1-5)')
    ax6.set_title('EM Safety Risk Assessment', fontsize=11)
    ax6.set_xlim(0, 6)
    ax6.axvline(x=3, color='red', linestyle='--', alpha=0.5, label='Acceptable Threshold')
    ax6.legend()
    ax6.grid(True, alpha=0.3, axis='x')
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig4_em_safety.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ Figure 4 saved: EM Safety Assessment")


# =============================================================================
# 仿真5：脑机接口系统综合仿真
# =============================================================================
def simulation_5_bci_system():
    """脑机接口系统综合仿真"""
    print("[5/5] 运行脑机接口系统综合仿真...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 子图1: 系统架构图
    ax1 = axes[0, 0]
    
    # 绘制系统框图
    blocks = [
        {'name': 'Neural\nSignals', 'pos': (1, 5), 'color': 'lightblue'},
        {'name': 'Electrode\nArray', 'pos': (3, 5), 'color': 'lightgreen'},
        {'name': 'Signal\nConditioning', 'pos': (5, 5), 'color': 'lightyellow'},
        {'name': 'ADC', 'pos': (7, 5), 'color': 'lightcoral'},
        {'name': 'DSP', 'pos': (9, 5), 'color': 'plum'},
        {'name': 'Wireless\nTx', 'pos': (11, 5), 'color': 'wheat'},
        {'name': 'External\nReceiver', 'pos': (13, 5), 'color': 'lightgray'},
    ]
    
    for block in blocks:
        rect = Rectangle((block['pos'][0]-0.8, block['pos'][1]-0.5), 1.6, 1,
                         facecolor=block['color'], edgecolor='black', linewidth=2)
        ax1.add_patch(rect)
        ax1.text(block['pos'][0], block['pos'][1], block['name'],
                ha='center', va='center', fontsize=8)
    
    # 绘制连接线
    for i in range(len(blocks) - 1):
        ax1.annotate('', xy=(blocks[i+1]['pos'][0]-0.8, blocks[i+1]['pos'][1]),
                    xytext=(blocks[i]['pos'][0]+0.8, blocks[i]['pos'][1]),
                    arrowprops=dict(arrowstyle='->', color='blue', lw=2))
    
    # 添加电源模块
    power = Rectangle((5, 2.5), 6, 1, facecolor='gold', edgecolor='black', linewidth=2)
    ax1.add_patch(power)
    ax1.text(8, 3, 'Wireless Power Transfer', ha='center', va='center', fontsize=9)
    
    # 电源连接线
    ax1.plot([8, 8], [3.5, 4.5], 'r-', linewidth=2)
    ax1.plot([7, 8], [4.5, 4.5], 'r-', linewidth=2)
    
    ax1.set_xlim(0, 14)
    ax1.set_ylim(1, 7)
    ax1.set_aspect('equal')
    ax1.set_title('BCI System Architecture', fontsize=11)
    ax1.axis('off')
    
    # 子图2: 信号处理链路
    ax2 = axes[0, 1]
    
    # 模拟信号处理流程
    t = np.linspace(0, 1, 1000)  # s
    
    # 原始信号（含噪声）
    signal_raw = np.sin(2*np.pi*10*t) + 0.5*np.sin(2*np.pi*50*t) + 0.3*np.random.randn(len(t))
    
    # 滤波后信号
    from scipy import signal as sig
    b, a = sig.butter(4, [8, 12], btype='band', fs=1000)
    signal_filtered = sig.filtfilt(b, a, signal_raw)
    
    # 特征提取（包络）
    analytic_signal = sig.hilbert(signal_filtered)
    envelope = np.abs(analytic_signal)
    
    ax2.plot(t, signal_raw, 'b-', alpha=0.3, linewidth=1, label='Raw Signal')
    ax2.plot(t, signal_filtered, 'g-', linewidth=2, label='Filtered (α band)')
    ax2.plot(t, envelope, 'r--', linewidth=2, label='Envelope')
    
    ax2.set_xlabel('Time (s)')
    ax2.set_ylabel('Amplitude (μV)')
    ax2.set_title('Signal Processing Chain', fontsize=11)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 1)
    
    # 子图3: 通信链路预算
    ax3 = axes[0, 2]
    
    components = ['Tx Power', 'Antenna Gain', 'Path Loss', 'Rx Gain', 'Rx Power']
    values = [0, -10, -40, 10, -40]  # dBm/dB
    cumulative = np.cumsum([0, -10, -40, 10, 0])  # 累积
    
    x_pos = np.arange(len(components))
    
    ax3.bar(x_pos, values, color=['blue', 'green', 'red', 'green', 'orange'],
           alpha=0.7, edgecolor='black')
    
    # 添加数值标签
    for i, (val, cum) in enumerate(zip(values, cumulative)):
        ax3.text(i, val + (2 if val > 0 else -4), f'{val} dB', ha='center', fontsize=9)
    
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels(components, rotation=15, ha='right')
    ax3.set_ylabel('Power / Gain (dB)')
    ax3.set_title('Communication Link Budget', fontsize=11)
    ax3.grid(True, alpha=0.3, axis='y')
    ax3.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
    ax3.set_ylim(-60, 20)
    
    # 子图4: 系统功耗分析
    ax4 = axes[1, 0]
    
    modules = ['Amplifier', 'ADC', 'DSP', 'RF Tx', 'Clock', 'Other']
    power_consumption = [20, 15, 30, 25, 5, 5]  # μW
    
    colors = plt.cm.Set3(np.linspace(0, 1, len(modules)))
    
    wedges, texts, autotexts = ax4.pie(power_consumption, labels=modules, colors=colors,
                                       autopct='%1.1f%%', shadow=True, startangle=90)
    
    total_power = sum(power_consumption)
    ax4.set_title(f'System Power Consumption\nTotal: {total_power} μW', fontsize=11)
    
    # 子图5: 性能指标雷达图
    ax5 = axes[1, 1]
    
    categories = ['Signal\nQuality', 'Power\nEfficiency', 'Data\nRate', 'Safety', 'Reliability', 'Size']
    values = [8, 7, 6, 9, 7, 8]  # 1-10评分
    
    # 闭合雷达图
    angles = np.linspace(0, 2*np.pi, len(categories), endpoint=False).tolist()
    values_plot = values + values[:1]
    angles += angles[:1]
    
    ax5 = plt.subplot(2, 3, 5, projection='polar')
    ax5.plot(angles, values_plot, 'o-', linewidth=2, color='blue', label='Current Design')
    ax5.fill(angles, values_plot, alpha=0.25, color='blue')
    
    # 添加目标值
    target = [9, 8, 8, 10, 8, 9]
    target_plot = target + target[:1]
    ax5.plot(angles, target_plot, 'o--', linewidth=2, color='red', alpha=0.5, label='Target')
    
    ax5.set_xticks(angles[:-1])
    ax5.set_xticklabels(categories, fontsize=9)
    ax5.set_ylim(0, 10)
    ax5.set_title('System Performance Metrics', fontsize=11, pad=20)
    ax5.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1))
    
    # 子图6: 系统性能随时间变化
    ax6 = axes[1, 2]
    
    time_months = np.linspace(0, 24, 100)  # 24个月
    
    # 模拟性能衰减（生物相容性、电池、信号质量等）
    signal_quality = 90 * np.exp(-time_months/36) + 10  # 信号质量
    power_efficiency = 85 * np.exp(-time_months/30) + 5  # 功率效率
    biocompatibility = 95 - 5 * (1 - np.exp(-time_months/18))  # 生物相容性
    
    ax6.plot(time_months, signal_quality, 'b-', linewidth=2, label='Signal Quality')
    ax6.plot(time_months, power_efficiency, 'r--', linewidth=2, label='Power Efficiency')
    ax6.plot(time_months, biocompatibility, 'g:', linewidth=2, label='Biocompatibility')
    
    ax6.axhline(y=80, color='orange', linestyle='--', alpha=0.5, label='Threshold')
    ax6.axvline(x=12, color='purple', linestyle=':', alpha=0.5)
    ax6.text(12.5, 50, '1 Year', fontsize=9, color='purple')
    
    ax6.set_xlabel('Time (months)')
    ax6.set_ylabel('Performance (%)')
    ax6.set_title('Long-term Performance Degradation', fontsize=11)
    ax6.legend(fontsize=8)
    ax6.grid(True, alpha=0.3)
    ax6.set_xlim(0, 24)
    ax6.set_ylim(0, 100)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig5_bci_system.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ Figure 5 saved: BCI System Integration")


# =============================================================================
# 主程序
# =============================================================================
if __name__ == '__main__':
    print("=" * 70)
    print("主题099：脑机接口电磁设计 - 仿真程序")
    print("=" * 70)
    print()
    
    # 运行所有仿真
    simulation_1_neural_signal_acquisition()
    simulation_2_implantable_antenna()
    simulation_3_wireless_power_transfer()
    simulation_4_em_safety()
    simulation_5_bci_system()
    
    print()
    print("=" * 70)
    print("所有仿真完成！")
    print(f"结果保存在: {output_dir}")
    print("=" * 70)

7.3 仿真结果分析

仿真1：神经信号采集系统

仿真结果展示了神经信号采集系统的关键性能指标：

动作电位波形：模拟了神经元动作电位的产生过程，展示了从静息电位（-70mV）到去极化峰值（约+30mV）再到复极化的完整过程。信号中包含多个神经元的放电活动和随机噪声。

脑电波频谱：展示了EEG信号的主要频段分布，包括δ波（0.5-4Hz）、θ波（4-8Hz）、α波（8-13Hz）、β波（13-30Hz）和γ波（30-100Hz）。α波在清醒放松状态下最为显著。

电极阻抗特性：基于Randles等效电路模型，展示了电极阻抗随频率的变化。在低频段，阻抗主要由电荷转移电阻主导；在高频段，双电层电容起主要作用。

信噪比分析：分析了不同噪声水平下的信号质量。当噪声水平低于10μVrms时，SNR可达到20dB以上，满足神经信号采集的要求。

多通道阵列：展示了4×4电极阵列的布局设计，通道间距为2mm，可实现较高的空间分辨率。

仿真2：植入式天线设计

PIFA天线结构：展示了平面倒F天线的基本结构，包括辐射贴片、短路针和馈电点。这种结构紧凑，适合植入式应用。

组织电磁特性：展示了脑组织的相对介电常数和电导率随频率的变化。在MICS频段（402MHz），介电常数约为50，电导率约为0.5S/m。

回波损耗：天线在组织环境中的回波损耗曲线，-10dB带宽约为20MHz，满足MICS频段通信需求。

辐射方向图：展示了天线在组织中的辐射特性，由于组织的吸收作用，天线呈现准全向辐射模式。

增益与效率：植入式天线在组织中的增益约为-12dBi，效率约为10%，这是典型的植入式天线性能水平。

仿真3：无线能量传输

耦合系数：展示了发射和接收线圈之间的耦合系数随距离的变化。在20mm距离处，耦合系数约为0.15，这是无线能量传输的典型工作点。

传输效率：分析了不同耦合系数和品质因数下的最大传输效率。当Q=50、k=0.15时，理论最大效率约为50%。

整流效率：展示了整流电路效率随输入功率的变化。在0dBm输入功率时，整流效率可达70%以上。

功率预算：展示了无线能量传输系统的功率分配，从发射端100mW到接收端负载40mW，总效率约为40%。

频段选择：对比了不同频段的特性，MICS频段（402MHz）在穿透能力和能量传输效率之间取得了良好平衡。

仿真4：电磁安全性评估

SAR分布：展示了植入设备周围的SAR分布，峰值SAR为1.5W/kg，低于IEEE C95.1标准的1.6W/kg限值。

温度升高：基于SAR计算了组织温升分布，最大温升约为2°C，在安全范围内。

安全标准对比：对比了不同安全标准的SAR限值，ICNIRP和IEEE C95.1是最常用的标准。

暴露时间：展示了SAR限值随暴露时间的变化，短时暴露允许更高的SAR值。

组织特性：对比了不同组织在402MHz频率下的电磁特性，脑脊液（CSF）具有最高的介电常数和电导率。

风险评估：对各项风险因素进行了评估，组织类型和SAR水平是需要重点关注的风险因素。

仿真5：脑机接口系统综合

系统架构：展示了完整的脑机接口系统框图，包括神经信号采集、信号调理、模数转换、数字信号处理、无线传输和外部接收等模块。

信号处理链路：展示了从原始信号到滤波后信号再到特征提取的完整处理流程。

通信链路预算：分析了无线通信链路的功率预算，确保可靠的信号传输。

功耗分析：展示了系统各模块的功耗分布，总功耗约为100μW，满足植入式设备的低功耗要求。

性能指标：通过雷达图展示了系统的综合性能，在信号质量和安全性方面表现良好。

长期性能：模拟了系统性能随时间的衰减，24个月后主要性能指标仍可维持在80%以上。

8. 应用案例

8.1 案例一：运动皮层脑机接口

应用场景：帮助瘫痪患者控制机械臂

系统设计：

电极阵列：96通道Utah阵列，植入运动皮层
采样率：30kHz/通道
无线传输：MICS频段，数据率200kbps
无线供电：感应耦合，传输距离15mm

性能指标：

信号质量：SNR > 15dB
解码准确率：> 90%
系统功耗：< 50mW
连续工作时间：> 24小时

8.2 案例二：视觉皮层刺激器

应用场景：为盲人提供人工视觉

系统设计：

电极阵列：1024通道，植入视觉皮层
刺激模式：双相电流脉冲
无线通信：双向数据传输
安全设计：严格的 charge-balanced 刺激

性能指标：

刺激分辨率：100×100像素等效
刺激频率：单通道最高100Hz
功耗预算：刺激模块< 20mW
安全性：电荷密度< 30μC/cm²

8.3 案例三：深部脑刺激（DBS）

应用场景：治疗帕金森病、癫痫等神经系统疾病

系统设计：

刺激靶点：丘脑底核（STN）或苍白球内侧（GPi）
电极设计：四触点电极，可编程刺激
通信方式：感应遥测
电源：可充电电池+无线充电

性能指标：

刺激参数：幅度0-10V，脉宽60-450μs，频率2-250Hz
电池寿命：> 5年
充电时间：< 2小时
疗效：运动症状改善> 50%

附录：关键参数汇总

A.1 神经信号参数

参数	典型值	单位
静息电位	-70	mV
动作电位幅度	100	mV
动作电位持续时间	1-2	ms
LFP频率范围	1-100	Hz
Spike频率范围	300-3000	Hz
EEG幅度	5-100	μV

A.2 植入式天线参数

参数	典型值	单位
工作频率	402/2450	MHz
增益	-15 ~ -10	dBi
效率	5-15	%
尺寸	10-20	mm
带宽	10-30	MHz

A.3 无线能量传输参数

参数	典型值	单位
传输距离	10-30	mm
耦合系数	0.1-0.3	-
线圈Q值	30-100	-
传输效率	10-40	%
接收功率	10-100	mW

A.4 安全限值参数

参数	限值	单位
头部SAR	1.6	W/kg
全身SAR	0.08	W/kg
局部温升	< 2	°C
MICS频段EIRP	25	μW

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

Web之深入解析Cookie的安全防御与跨域实践

Cookie是网站存储在用户浏览器中的小型文本数据，用于识别用户身份和记录偏好。它通过键值对存储信息，并包含Domain、Path、Expires等属性控制作用范围和有效期。常见应用场景包括用户登录状态保持、购物车存储和A/B测试分组。Cookie可由服务端或客户端生成修改，但受HttpOnly和SameSite等安全属性限制。默认遵循同源策略，但在特定条件下可实现跨域共享。使用时需注意编码、路径