一种基于FPGA和DAC的实时可重构通用忆阻器仿真方法

一种基于fpga和dac的实时可重构通用忆阻器仿真方法
技术领域
1.本发明属于忆阻器仿真技术领域，更为具体地讲，涉及一种基于fpga和dac的实时可重构通用忆阻器仿真方法。


背景技术：

2.1971年，蔡少棠教授首次根据电路完备性定理提出了忆阻器的理论模型。蔡少棠教授认为忆阻器是除电阻、电容和电感以外的第四个基本二端电路元件，它描述了电荷和磁通量之间的非线性关系。同时，蔡少棠教授也给出了忆阻器的3个基本特征：
①
在双极性周期性电信号的激励下，器件在v-i平面的曲线是一个捏滞回线；
②
当电信号扫描频率增大时，捏滞回线的波瓣面积持续减小；
③
当扫描频率趋近无穷大时，捏滞回线将收缩为一条单值函数。
3.直到2008年，因惠普实验室的stan williams团队首次在基于二氧化钛(tio2)材料的纳米级薄膜中实现了物理忆阻器，忆阻器的特性及其应用研究才得到了大量学者的关注。此时，蔡少棠教授提出的忆阻器不再是一个理论数学模型，而是一个实际存在的器件。目前忆阻器模型已经广泛应用到神经网络、机器学习、混沌理论、保密通信、图形加密、非易失性存储器、滤波电路等领域中。但不幸的是，一方面，忆阻器不同的制备原材料对应着各异的物理机制，使得忆阻特性也各不相同，这就限制了忆阻器的大范围推广应用。例如，研究人员已经在二元氧化物、复杂钙钛矿氧化物、固态电解质材料、非晶碳材料、有机高分子材料等各种材料中都发现捏滞回线这一特征现象，提出各异的物理机理来解释其忆阻特性，如氧空位迁移导致的导电通道形成和断裂、界面势垒调制、活性电极金属化反应导致的金属导电通道的形成和断裂、注入载流子的捕获和释放，以及金属-绝缘体转变机制等。另一方面，因忆阻器生产成本过高，目前，市场上还不能直接购买商用忆阻器。例如，knowm忆阻器被实现为商业组件，但由于其结构复杂且成本高，难以大范围推广。面对忆阻器复杂多变的物理机理和高昂的生产成本，设计一个面向多种物理模型的实时可重构通用忆阻器仿真电路，以用于通过实验探索方法解析忆阻器相关特性仍然是非常急迫的和重要的。
4.一直以来，相关学者对提高忆阻器仿真电路的带宽(频率)进行了大量的研究。模拟、数字模拟混合和专用集成电路等方法被用来构建忆阻器仿真电路并进行验证和带宽测试。
5.经过大量学者深厚的理论研究和专业的硬件调试能力，模拟方法主要使用大量无源器件(如电阻、电容、电感)和有源器件(如运算放大器、跨导运算放大器(ota)、电流反馈运算放大器(cfoa)、差分电流传送器(ddcc)和模拟乘法器)进行面包板或者电路板实验。忆阻器仿真电路的实验带宽已经由早期的500hz发展到当前的1.3mhz。同时，忆阻器仿真电路采用正负电源供电。在模拟方法中，对于不同的忆阻器模型和带宽，往往需要重新设计忆阻器仿真电路并花费大量时间进行电路调试，特别是针对mhz的输入信号，采用面包板的方式已经无法验证实验结果，因为信号传输线的寄生参数、阻抗匹配和信号间的串扰问题必须解决。电路板输出的信号可以采用示波器的李沙育图功能显示捏滞回线。值得注意的是模
拟方法通常存在较大误差，误差的来源主要如下：大多数商用电阻、电容、电感的误差大约分别是5％，20％，20％。运算放大器输入电压失调和输入偏置电流会对结果产生误差。有源器件模拟带宽的非线性特征对不同信号幅度频率曲线的影响较大。根据现有文献报道，目前模拟方法还没有实现带宽10mhz及以上的忆阻器仿真电路。
6.由于电路设计和调试需要非常专业的电路知识，这对于大多数从事基础研究的学者来说非常困难。所以相关学者又提出了数字模拟混合方法。数字模拟混合方法通常采用adc、可编程通用处理器、dac的方式来实现忆阻器仿真电路，其中，adc和dac分别用于输入信号的模拟-数字量化和输出信号的数字-模拟转换，可编程通用处理器用于根据数字模型实时计算忆阻或忆导。和模拟方法相比，数字模拟方法实现更加简单和快速，但数字模拟方法也存在如下限制。可编程通用处理器的计算速度、adc和dac工作频率和模拟带宽限制了数字模拟混合方法，不能仿真高工作频率的忆阻器。而adc、dac和数字示波器的分辨率以及dac输出信号的范围限制了数字模拟混合方法不能对忆阻器模型进行高精度硬件仿真。
7.对于已知的忆阻器模型，相关学者又提出了专用集成电路方法。在文献u.e.ayten,s.minaei,and m.“memristor emulator circuits using single cbta,”aeu-int.j.electron.commun.,vol.82,pp.109
–
118,dec.2017,doi：10.1016/j.aeue.2017.08.008中，使用23个cmos晶体管在
±
0.9v的电源下实现了忆阻器仿真电路。在文献a.yesil,“a new grounded memristor emulator based on mosfet-c,”aeu-int.j.electron.commun.,vol.91,pp.143
–
149,jul.2018,doi：10.1016/j.aeue.2018.05.004.中，介绍了仅由7个mos管和1个接地电容组成的忆阻器仿真电路。电路采用tsmc 0.18μm工艺参数的cadence environment布局，其布局尺寸仅12μm x 38μm，不包括电容器面积。布局后的忆阻器模拟仿真与理论分析非常吻合。显然采用专用集成电路的方法可以大幅缩小忆阻器仿真电路的尺寸并提高实现的精度。但复制的设计流程和高昂的电路制作成本是一个阻止改方法推广的一个主要因素。专用集成电路方法仅适用于一个已知的模型，而对于研究多模型下的忆阻器特征并不合适。
8.总的来讲，首先，大多数现有忆阻器仿真电路只能在相对较低的频率下模拟忆阻器的动态行为，即这些忆阻器仿真电路在低于临界频率时表现出收缩的磁滞回线，并在超过该临界频率时成为线性电阻。这限制了忆阻器仿真电路在高速度高带宽场合的应用，如高频随机信号产生、高速数据传输、高速数据存储。其次，现有忆阻器仿真电路都是面向固定模型设计。从理论模型到实际电路输出结果需要花费研究人员大量的时间和精力，大多数情况下需要具有专业电路背景的工程师参与硬件调试过程。现有的三种方法都无法根据模型实现实时动态可重构，所以硬件实现忆阻器仿真电路是一个耗时和低效率的繁琐过程。最后，忆阻器仿真电路中元器件固有误差、电路传输损耗、测量仪器误差对于忆阻器理论模型的高精度电路实现是非常重要的，但目前，忆阻器仿真电路的实验精度并没有引起相关学者的重视。


技术实现要素：

9.本发明的目的在于克服现有技术的不足，提供一种基于fpga的实时可重构通用忆阻器仿真方法，实现实时可重构以适应不同模型的忆阻器，并能够仿真高工作频率的忆阻器，同时提高忆阻器仿真电路的实验精度。
10.为实现上述发明目的，本发明基于fpga的实时可重构通用忆阻器仿真方法，其特征在于，包括以下步骤：
11.(1)、建立忆阻器的数学模型f(h[n])，判断其是否是关于系统状态变量h[n]的多项式，如果不是，进入步骤(2)；如果是，则确定数学模型f(h[n])关于系统状态变量h[n]的阶次m，并进入步骤(5)；
[0012]
(2)*根据输入信号的零直流分量交流信号的幅度和频率，分别确定最大幅度a
max
和最小频率ω
min
，并进一步确定系统状态变量h[n]的取值区间：
[0013][0014]
(3)、在系统状态变量h[n]的取值区间内，采用麦克劳林公式对数学模型即忆导值或忆阻值f(h[n])进行关于系统状态变量h[n]的m阶多项式拟合，得到拟合函数g(h(n))，并计算最大拟合相对误差εm：
[0015][0016]
设ε0为可接受的最大拟合相对误差，则多项式阶数m的取值应满足εm≤ε0；
[0017]
(4)、根据麦克劳林公式确定m+1个多项式系数ki,i＝0,1,2,...,m；
[0018]
(5)、计算归一化系数α1：
[0019][0020]
其中，r
ch
为当前设置的模拟通道的量程，且模拟通道的量程需要满足约束条件α1≤1，同时α1的取值范围为α1∈[-1,1]；
[0021]
(6)、将h(n)＝a
max
/ω
min
代入：
[0022][0023]
得到f(h[n])的最大值f
max
(h[n])；
[0024]
再将f(h[n])＝f
max
(h[n])和x(n)＝α1代入：
[0025]
y]＝n]＝f(h[n])
·
x[n]
[0026]
得到模型输出y[n]的最大值y
max
；
[0027]
其中，输入信号x[n]为adc量化得到的有符号的数据；
[0028]
(7)、基于fpga进行忆阻器实时仿真即在fpga中进行以下运算：
[0029]
7.1)、对于输入信号x[n]，首先使用fpga内部的定点数转浮点数ip核将其转换为单精度浮点数据f_x[n]，其中的取值范围为f_x[n]∈[-1,1]，然后计算系统状态变量h[n]即n时刻的磁通量或电荷量：
[0030]
[0031]
其中，ts为时间间隔，f_x[j]为输入信号的j个采样点，h[0]为系统状态变量的初始值；
[0032]
7.2)*计算忆导值或忆阻值f(h[n])：
[0033][0034]
7.3)、将单精度浮点数据f_x[n]同时送入fifo中进行延时处理，使得fifo读端口的数据即延时数据f_dly_x[n]与忆导值或忆阻值f(h[n])在时间上对齐，并计算输出信号y[n]：
[0035]
y[n]＝f(h[n])
·
f_dly_x[n]
[0036]
其中，输入信号x[n]为电压信号或电流信号，输出信号y[n]为电流信号或电压信号；
[0037]
7.4)、在fpga中，利用一个除法器ip核，将输出信号y[n]除以最大值y
max
，得到输出信号div2[n]：
[0038][0039]
对延时数据f_dly_x[n]进行延时处理，使得延时后的数据f_dly1_x[n]与div2[n]在时间是对齐的；
[0040]
7.5)、在fpga中，利用两个乘法器，将数据f_dly1_x[n]、div2[n]均与相乘，得到输出：
[0041][0042][0043]
其中，n
dac
是用于模拟波形输出dac的分辨率，且dac能够输出双极性波形；
[0044]
7.6)、在fpga中，利用dac将输出信号multi1和multi2转换为对应的模拟信号x
′
(t)、y
′
(t)；
[0045]
(8)、将模拟信号x
′
(t)、y
′
(t)输入到数字示波器中进行李沙育图形显示，显示忆阻器的捏滞迟滞回线。
[0046]
本发明的目的是这样实现的。
[0047]
本发明基于fpga的实时可重构通用忆阻器仿真方法，利用weierstrass定理和麦克劳林公式，将忆阻器数学模型通过m项多项式进行非线性拟合，其中m与输入信号的幅度和频率及拟合精度有关，这样通过更新多项式的阶次、多项式系数和采样间隔即可简单快速地适配指定的忆阻器模型。在此基础上，基于fpga进行忆阻器实时仿真：计算出系统状态变量、忆导值或忆阻值的基础上计算出输出信号，然后，对输入信号、输出信号进行归一化处理、dac输入处理以及dac数模转换，得到对应的模拟信号，最后送入数字示波器进行李沙育图形显示，以清晰显示忆阻器的捏滞迟滞回线。本发明实时可重构通用忆阻器的仿真方法，通过改变多项式系数即可实时可重构忆阻器，从而实现实时可重构以适应不同模型忆阻器的发明目的。并且本发明可通过多项式系数设置，采用低频信号实验等效验证高频信号的忆阻特性即可以仿真高工作频率的忆阻器，这一突出的优势使得研究忆阻器的特性不再需要特别关心仿真电路的模拟带宽，这对于大多数不具备专业电路调试能力的研究人员
是非常方便的。同时，本发明基于fpga进行重构仿真，实验精度得到了提高。
附图说明
[0048]
图1是二阶忆阻器模型不同系数ki绘制的q-f(q)曲线图，其中，(a)为系数k2≠0,k1≠0,k0≠0,f(q)单极性，(b)为系数k2≠0,k1≠0,k0≠0,f(q)双极性，(c)为系数k2≠0,k1＝0。f(q)单/双极性，(d)为系数k2≠0,k1≠0,k0＝0,f(q)双极性、(e)为系数k2＝0,k1≠0,k0≠0,f(q)双极性，(f)为系数k2＝0,k1≠0,k0＝0，f(q)双极性；
[0049]
图2是ω-l的近似变化曲线图；
[0050]
图3是数值仿真实验1系数矩阵k＝k1时电压电流特性曲线图，其中，(a)为输入电压时域波形图，(b)输出电流时域波形图，(c)为u-i曲线；
[0051]
图4是数值仿真实验1不同激励频率和幅度形成的特性曲线比较图，其中，(a)输入信号幅度为1v，(b)输入信号幅度为2v：
[0052]
图5是数值仿真实验2中，不同系数迟滞曲线对比图，其中，(a)为系数矩阵k＝k2，(b)为系数k0放大10倍，c)为系数k1放大10倍，(d)为系数k2放大10倍；
[0053]
图6是数值仿真实验3中采样率分别为100msa/s和200gsa/s下，输入不同的幅度和频率信号得到的图，其中，(a)为采样率100msa/s、幅度为0.1v的五个频率信号的捏滞迟滞回线图，(b)为采样率100msa/s、幅度为1v的五个频率信号的捏滞迟滞回线图，(c)为采样率100msa/s、幅度为2v的五个频率信号的捏滞迟滞回线图，(d)为采样率200gsa/s、幅度为0.1v的五个频率信号的捏滞迟滞回线图，(e)为采样率200msa/s、幅度为1v的五个频率信号的捏滞迟滞回线图，(f)为采样率200msa/s、幅度为2v的五个频率信号的捏滞迟滞回线图；
[0054]
图7是基于过采样的捏滞迟滞回线图，其中，(a)对应500khz的正弦输入信号；(b)对应1mhz的正弦输入信号，(c)对应2mhz的正弦输入信号；
[0055]
图8是数值仿真实验4中忆阻器的具有有源特性仿真，其中，(a)为系数矩阵k
4.0
下的捏滞迟滞回线图，(b)为系数矩阵k
4.1
下的捏滞迟滞回线图,(c)为系数矩阵k
4.2
下的捏滞迟滞回线图，(d)为系数矩阵k
4.3
下的捏滞迟滞回线图；
[0056]
图9是数值仿真实验5中改变系数k对迟滞回路面积影响对比图，其中，(a)为系数矩阵k
5.0
下的捏滞迟滞回线图，(b)为系数矩阵k
5.1
下的捏滞迟滞回线图,(c)为系数矩阵k
5.2
下的捏滞迟滞回线图，(d)为系数矩阵k
5.3
下的捏滞迟滞回线图；
[0057]
图10是数值仿真实验6中200gsa/s和100msa/s各参数时域波形及对比图，其中，(a)对应200gsa/s，(b)对应100msa/，(c)是两者合并在一起；
[0058]
图11是本发明基于fpga的实时可重构通用忆阻器仿真方法一种具体实施方式流程图；
[0059]
图12是基于fpga的实时可重构通用忆阻器仿真硬件实现部分原理示意图；
[0060]
图13是本发明基于fpga的实时可重构通用忆阻器仿真方法一仿真时序图；
[0061]
图14是k＝k0时，采用matlab软件实现和本发明仿真的结果对比图，其中，(a)为捏滞迟滞回线对比图，(b)为误差局部放大对比图。
具体实施方式
[0062]
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地
理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
[0063]
1、忆阻器通用数学模型
[0064]
自从蔡少棠教授在1971年根据电路完备性理论推导出忆阻器的本构关系，许多科学家都对忆阻器数学模型进行了研究。2008年hp实验室在实验室观察到了忆阻显现，这一发现实现了忆阻器数学模型到物理实现的推进。
[0065]
忆阻器本质是一个非线性动力学系统，一个时不变的x控制的忆阻系统也称为一个广义忆阻，定义为：
[0066][0067]
式中，f(x(t),h)是一个关于个自变量的连续有界函数。其中x(t)，y(t)分别是广义忆阻器系统的输入电压(电流)信号和输出电流(电压)信号，f(h,m,t)表示忆导值(忆阻值)，h是内部n维态变量，为状态矢量的时间导数。若x＝u(t)，则公式(1)称为一个广义电压控制忆阻器。若x＝i(t)，则公式(1)称为一个广义电流控制忆阻器。
[0068]
根据定理1，当f(x(t),h)可用m阶多项式函数序列fm(x(t),h)一致逼近时，公式(1)所表示的广义忆阻器模型可进一步表示为如下形式。
[0069][0070]
其中k＝(k0,k1,...,ki,...,km)是系数矩阵，且k≠0。h(x,t)表示磁通量(电荷量)，h(t0)表示磁通量(电荷量)的初始值。参数α是y(t)的幅度因子。当x(t)为电压(或电流)时，公式(2)则表示电压控制忆阻器(或电流控制忆阻器)。从数学模型上来看，两类忆阻器也仅仅是输入变量x(t)、输出变量y(t)和系统变量h(x,t)物理单位的不同。因此，为了简化表达，本发明都以电压控制忆阻器模型为例进行分析，并固定α＝1。
[0071]
公式(2)是提出忆阻器通用数学模型的连续模式表达式，而连续模型无法直接在通用数字处理电路中进行处理。以均匀时间间隔ts对公式(2)进行量化处理，得到对应的忆阻器离散模型如下：
[0072][0073]
其中，ts是步进时间，n是时刻，也是数字积分窗口长度，为了表达形式一致，ki重新记为k[i]，其他变量与公式(1)中一致。根据黎曼积分(riemann integral)的定义，公式(2)
中的第3个表达式的离散数学模型是公式(3)第3个表达式。
[0074]
黎曼积分的定义为：函数f(x)在闭区间[a,b]有定义，在区间[a,b]上插入n-1个分点，使其分成n个小区间[x
i-1
,xi]，i＝1,2,3,...,n，任取一点ξi∈[x
i-1
,xi]，做和式为：
[0075][0076]
设λ为所有小区间的最大长度，即λ＝max{δx1,δx2,...,δxn}。当分割越来越小的时候，即λ
→
0，n
→
+∞，分割后的每个小条就可以看作矩形，考虑式子：
[0077][0078]
如果上式值存在，且此极限不依赖于区间[a,b]的分法，也不依赖于点ξi的取法，则称此极限为f(x)在区间[a,b]上的定积分，记为：
[0079][0080]
在实际工程应用中，使用模拟数字转换器(adc)在采样时钟f
adc
(1/f
adc
＝ts)的驱动下实现对连续输入信号x(t)的均匀量化。当采样时钟f
adc
频率越高，均匀时间间隔ts也就越小。也就是黎曼积分中的分割步长越来越精细。设fpga系统工作时钟是f
fpga
。当公式(3)在fpga中实现时，如果f
fpga
≥f
adc
，此时在fpga中可实现忆阻器模型的实时计算。当f
fpga
＜f
adc
时，可采用fpga内部存储资源如fifo或外接存储器如ddr3，先对adc量化数据进行缓存，再进行忆阻器模型的硬件实现。根据奈奎斯特采样定理有f
adc
≥2f
in
，所以adc的采样率越高，输入信号的带宽越大，fpga中实时数据处理带宽越大。
[0081]
2、忆阻器通用数学模型的通用性
[0082]
weierstrass逼近定理：设f(x)是[a,b]上的连续函数，则存在多项式函数列{fn(x)}，使得{fn(x)}一致收敛于f(x)。证明如下：
[0083]
由bernstein多项式的定义，在[0,1]上，误差函数满足：
[0084][0085]
由于f(x)∈c[0,1]，因此f(x)在[0,1]上有界且一致连续，设存在正数m,使得在[0,1]上|f(x)|＜m，因此：
[0086][0087]
由f(x)在[0,1]上一致连续可知，对任意的ε＞0，都存在δ(ε)＞0，使得当时，
[0088][0089]
将公式(7)式拆分为两部分：
[0090][0091]
由一致连续性公式(9)，可得：
[0092][0093]
再利用一致有界性公式(8)，当时，
[0094][0095]
因此，当n满足时，由公式(10)-(12)，可得
[0096][0097]
这就证明了weierstrass第一定理。
[0098]
泰勒公式可在局部将复杂函数f(x)逼近近似地表示为n次多项式函数之和，泰勒公式这种化繁为简的功能，使得它成为分析和研究许多数学问题的有力工具。
[0099]
泰勒公式：若函数f(x)在包含x0的某个开区间(a,b)上具有(n+1)阶的导数，那么对于任一x∈(a,b)，有：
[0100][0101]
其中称为n次泰勒多项式，误差项称为n阶泰勒余项，ε为x与x0之间的某个值。
[0102]
如果函数f(x)的n+1阶导数在n(x0)上有界m，从而有：
[0103][0104]
表明rn(x)＝o((x-x0)n)，另外也可证明对固定的x，当n
→
∞时，rn(x)
→
0，即，要想使f(x)与pn(x)误差减小，则可将|x-x0|取小，也可将n取大。在n阶泰勒公式中，如果x0＝0，从而可得：
[0105][0106]
该式称为函数f(x)在x＝0处的n阶泰勒公式，也称作f(x)的n阶麦克劳林(maclaurin)公式，其余项常写为o(xn)或者两种形式，用n+1阶导数表示的余项叫拉格朗日余项，用o(xn)或者o((x-x0)n)表示的余项叫皮亚诺(peano)余项。
[0107]
在所提出的忆阻器模型中，输入信号的幅度a和频率ω是已知的，因此h(t)的取值范围是固定有界的。结合weierstrass第一定理f(h(t))都可以用一多项式进行拟合。而泰勒公式和麦克劳林(maclaurin)公式则给出了m阶多项式的具体实现方法。参数m取值越大，f(h(t))的逼近精度越高，f(h(t))的计算复杂度也就越高。所以在通用处理器如fpga中实现时，需要在精度、速度以及逻辑资源三个方面做出平衡。对于离散忆阻器模型，可以用公式(3)进行逼近。例如当a＜＜ω时，采用麦克劳林(maclaurin)公式对常见函数逼近时系数矩阵k取值情况如表1所示。从表1中可以看出，通过配置不用的系数矩阵k＝{ki}即可实现对常用函数忆阻数学模型的多项式逼近。
[0108][0109][0110]
表1
[0111]
所提出的忆阻器模型对现有文献中的模型也展示出了良好的通用性。例如，2008年，hp团队发现的忆阻器模型如公式(17)所示。
[0112][0113]
其中d为二氧化钛薄膜的厚度，w(t)是掺杂层的厚度。通过外加电压，w(t)也随之
改变，r
on
和r
off
分应忆阻器电阻的最小值(w(t)＝d)和最大值(w(t)＝0)。当系数矩阵k中的元素k0＝r
off
，k1＝(r
on-r
off
)/d，m＝1时，提出的忆阻器模型就是经典的hp忆阻器模型。同理分别配置系数矩阵k和多项式系数m对其他现有模型的通用性如表2所示，忆阻器模型对现有模型具有兼容性。
[0114][0115]
表2
[0116]
3、本发明忆阻器模型对有源或无源忆阻器的通用性
[0117]
在公式(1-2)中f(h(t))的极性与系数矩阵k和h(t)的值有关。为了详细说明各参数对q-f(q)曲线的影响，我们选择m＝2时的电荷控制忆阻器模型进行详细分析，其中q表示磁通量。当m＝2时，f(q)的数学表达式为：
[0118][0119]
系数k1、k2不能同时为0，否则忆阻器会变成线性的时不变电阻器。根据k0,k1和k2取值不同，q-f(q)曲线共有如图1所示的24种图形。从图1中可以看出，曲线q-f(q)可出现在横轴的上方，下方或者与横轴有交点。
[0120]
如果一个忆阻器在任意零直流分量周期性信号作用下，捏滞迟滞回线只在第一和第三象限有分布，我们称其为无源的。否则，称为有源的。有源忆阻器的捏滞回线可以不过零点。因为当捏滞回线落于伏安平面的第二第四象限或第四象限时,x(t)y(t)≤0，有源忆阻器的行为类似于一个能量源。在公式(2)中，等号右边的自变量f(q)和x(t)的极性决定了等号左边因变量y(t)的极性。当输入零直流分量的双极性信号x(t)时，u-i平面的回线有如表3所示的3种情况。第一，如果f(q)正极性，则u-i平面的回线经过第1象限和第3象限，此时提出的忆阻器属于无源忆阻器。第二，如果f(q)负极性，则u-i平面的回线经过第2象限和第
4象限，此时提出的忆阻器属于有源忆阻器。第三，如果f(q)是双极性，则u-i平面的回线经过所有象限，此时提出的忆阻器属于局部有源忆阻器。图1中的曲线(1)、(2)、(11)和曲线(3)、(4)、(14)对应的忆阻器模型分别属于无源忆阻器和有源忆阻器。而其他曲线对应的忆阻器模型都是局部有源忆阻器模型。
[0121][0122][0123]
表3
[0124]
4、本发明忆阻器模型的非线性特征分析
[0125]
为了验证公式(2)中改进型忆阻器数学模型的输入信号带宽，我们使用余弦信号作为输入进行理论带宽分析。设输入信号为x(t)＝acos(ωt)，其中a是输入信号的幅度，ω是输入信号的角频率。假设信号开始激励系统的时刻为t0，同时也假设有x(t0)＝0，h(t0)＝0。h(t)的数学模型如下：
[0126][0127]
其中，有因为所以我们可以很容易地计算出
[0128][0129]
将公式(19)带入公式(2)中，忆阻器h(t)的表达式如下：
[0130][0131]
公式(20)带入公式(2)的第一个表达式，我们就可以得到忆阻器的输入输出信号
表达式如下：
[0132][0133]
在公式(22)中，由关于x(t)的线性项y1(t)和非线性项y2(t)组成。其中y1(t)＝k0x(t)，
[0134]
当k0≠0时，在公式(22)中，随时输入信号角频率ω的减小，那么y2(t)的值将逐渐增大，此时y(t)的值主要由非线性项y2(t)决定。相反地，随时输入信号频率ω的增大，那么y2(t)的值将逐渐减小到0。此时y(t)的值主要取决于线性项y1(t)的大小，且y1(t)可近似理解为一个标准线性器件。比如，当输入为电压信号时，y(t)就是一个标准导纳模型。
[0135]
考虑实际应用的场景，产生更高频率的信号要比产生更高信号幅度要容易的多。比如常见信号源可以产生频率ω高达ghz的正弦信号，而幅度a通常才数十v。因此，在公式(19)中有0≤|h(t)|＜1。函数(|h(t)|)m(m≥1)是关于幂指数m的减函数，所以改变系数k0的绝对值，f(h(t))的变化量也就最大。同时，当n≥1时，也有max{(|h(t)|)m}≤|h(t)|。在系数k1到km共m个系数中，调整系数k1可使f(h(t))的变化量最大。换句话说，为了快速的增加或者减小y2(t)的非线性特征，在实际工程应用中，我们应当调整系数k0或者k1的变化量。而变化量相同的其他系数则不容易使捏滞曲线形状或面积发生变化。
[0136]
根据棣莫弗定理(de moivres theorem)有如下降幂公式成立，当m是偶数时，有：
[0137][0138]
当m是奇数时，有：
[0139][0140]
再结合三角函数的积化和差公式：
[0141][0142]
根据公式(23)-公式(25)，在公式(22)中有如下表达式成立。当m是偶数时，有：
[0143][0144]
当m是奇数时，有：
[0145]
[0146]
根据公式(26)和公式(27)，当输入为单频信号x(t)＝acos(ωt)时，输出信号y(t)是输入信号基波和各次谐波之和。其中基波幅度为ak0/ω，单个谐波的幅度同时受多项式系数ki(i＝1,...,m)，m及a/ω的影响。
[0147]
设则有：
[0148][0149]
ω-l的近似变化曲线如图2所示，从图中我们可以发现当信号角频率ω从0开始增加到无穷大时，l的取值由大变小，公式(2)的数学模型也从非线性模型变化为线性模型。换句话说，就是提出的忆阻器模型随着角频率ω的增加逐渐由非线性忆阻(忆导)模型变为线性电阻模型。假设当比值l小于等于阈值δ时，即0＜l≤δ，改进后的忆阻器已成为电阻器件。δ的取值则可以根据实际应用情况来选择以满足系统的精度要求。
[0150]
在图2中，我们把l＝δ对应的输入信号角频率ω
max
称为忆阻器最大理论带宽。因此为了实现忆阻器的迟滞回路特性，参数设置必须满足公式(29)中的约束条件。
[0151][0152]
公式(29)对提高改进型忆阻器模型的记忆特点包括如下几点：首先，当k0绝对值的减小，ki(i＞0)绝对值的增大，多项式的次数m的增加，输入正弦信号幅度a的增大或者角频率ω的减小，公式(2)则更容易满足忆阻器的非线性特性。其次，当输入正弦信号角频率ω已知时，为了使公式(2)的数学模型表现出忆阻器特征，我们可以增大输入信号的幅度a，增加ki(i＞0)绝对值或者减小的k0绝对值。
[0153]
4、提出的忆阻器模型对宽带信号的通用性
[0154]
根据数字采样原理，有：
[0155][0156]
其中fi和f
sap
分别是信号频率和采样频率。当相同幅度的低频信号f
il
和高频信号f
ih
分别使用采样率为f
sl
和f
sh
的adc进行量化，且满足约束条件：
[0157][0158]
则有低采样率量化序列x
l
[n]和高采样率量化序列xh[n]满足如下关系：
[0159]
x
l
[n]＝xh[n]。
ꢀꢀꢀ
(32)
[0160]
设两个情况下对应的系数矩阵分别为k
l
和kh。如果系数矩阵k
l
和kh满足公式(33)的约束条件(i∈[0,m])，则低频采样量化和高频采样量化形成的迟滞回路曲线完全一样。证明如下：
[0161]kil
＝ri·kih
ꢀꢀꢀ
(33)
[0162]
证明：设r
·fil
＝f
ih
，则有r-1
·
t
sl
＝t
sh
。其中t
sl
和t
sh
分别是f
sl
和f
sh
对应的采样周期。再根据公式(3)有：
[0163]hl
[n]＝r
·hh
[n]
ꢀꢀꢀ
(34)
[0164]
其中h
l
[n]和hh[n]分别是f
sl
和f
sh
对应的磁通量。进一步展开公式(3)有：
[0165]fl
(h,m,n)＝k
0l
+k
1l
·hl
[n]+...+k
il
·
(h
l
[n])i+k
ml
·
(h
l
[n])mꢀꢀꢀ
(35)
[0166]
和
[0167]fh
(h,m,n)＝k
0h
+k
1h
·hh
[n]+...+k
ih
·
(hh[n])i+k
mh
·
(hh[n])mꢀꢀꢀ
(36)
[0168]
将公式(33)-(34)带入公式(35)：
[0169][0170]
公式(36)和公式(37)等号右边相等，因此有：
[0171]fl
(h
l
[n])＝fh(hh[n])
ꢀꢀꢀ
(38)
[0172]
进一步地，根据公式(3)有：
[0173]yl
[n]＝yh[n]
ꢀꢀꢀ
(4)
[0174]
所以根据公式(32)和公式(39)有：低频采样量化和高频采样量化形成的迟滞回路曲线完全一样。证毕！
[0175]
根据上述结论，设osr＝20，则有输入信号f
ih
＝10ghz和采样率f
sh
＝200gsps形成的迟滞回路曲线，与输入信号f
il
＝5mhz和采样率f
sl
＝100msps形成的迟滞回路曲线完全一样。通过配置提出忆阻器模型的多项式系数k，将高频高采样率的数学模型搬移到低频低采样率进行处理和硬件验证，对忆阻器模型的工程化推广和应用具有非常重要的意义。原因主要如下：一是ghz模拟信号属于微波段高频信号，此时电路中的寄生参数对信号电路的模拟带宽和传输影响非常大，模拟电路的设计和调试难度也非常高，除非具有资深背景的专业公司可以完成硬件平台的实现，大多数从事理论研究的科学家们是无法在实验室进行硬件验证的。二是现在还无商用单片adc的采样率可以达到数百ghz。三是高采样率带来的数据带宽急速上升，后端数据处理fpga无法实现实时运算。比如8bit 200gsps的采样率实时数据带宽高达200gsa/s，现有商用fpga无法实现实时运算。
[0176]
5、改进忆阻器模型的数值仿真
[0177]
首先，matlab被用于对提出的忆阻器通用数学模型进行数值仿真，以证明所提出的忆阻器通用数学模型正确性。我们分别设计了6个数值仿真实验，每个实验的参数和实验结果分别如下所示。
[0178]
数值仿真实验1
[0179]
我们设置多项式系数矩阵k1＝(-0.4,4e6,-1500,-1.25e-6,-3.125e-9)，m＝5，ts＝10ns。考虑余弦激励电压信号幅度a＝1v和频率fi＝2mhz时，得到的时域波形曲线如图3所示。图3(a)、(b)分别显示了输入电压x(t)和输出电流y(t)的时域波形。图3(c)描述了改进忆阻器在u-i平面上的特性曲线，显示了一个典型的捏滞迟滞回线。
[0180]
再次改变余弦信号的幅度(1v和2v)和频率(1mhz，2mhz，2.5mhz，5mhz和10mhz)，当
采样率固定为100msa/s时，仿真结果如图4所示。图4(a)和图4(b)分别显示了当输入信号幅度为1v和2v时，不同频率的特性曲线的比较。很明显，随着激励频率的增加，磁滞回线的面积单调减小，磁滞回线的斜率变化也逐渐减小，越来越像直线。但u-i平面的曲线总是固定经过坐标原点且在原点处收紧的磁滞回线。另外，从图4中还可以看出：首先，当输入信号幅度不变时，迟滞回线的面积随着频率的增加而减小。其次，当输入信号频率不变时，捏滞迟滞回线的面积随着幅度的增加而增加。因此同时增大(减小)输入信号幅度，或减小(增大)输入信号频率，都可以快速增大(减小)捏滞迟滞回线的面积。综述所述，公式(2)中描述的忆阻器模型可以呈现用于识别忆阻器的基本特征。
[0181]
数值仿真实验2
[0182]
根据公式(29)，公式(2)中的系数ki与提出忆阻器模型的迟滞回线和工作频率带宽有关。首先设置正弦激励电压a＝1v和fi＝5mhz，当系数矩阵k2＝(1,5e6,5e6,5e6,7,100,500,0,-10,0,1)时,绘制的u-i平面的迟滞曲线如图5(a)所示。然后将分别将k0、k1和k2分别扩大10倍，绘制的迟滞曲线分别如图5(b)-(d)所示。从图5中可以看到，增加k0磁滞回线将缩小为线性函数，此时忆阻器的中的线性部分占比较大。而增加k1磁滞回线的面积将增加，此时忆阻器的中的非线性部分占比较大。而增加k2迟滞回路的面积几乎不变。这与第二节的理论计算是一致的。
[0183]
数值仿真实验3
[0184]
分别设置系数矩阵为：
[0185]k3.0
＝(-0.4,4*10^6,-1500,-1.26*10^-6,-20,0,0,0,0,0,0)
[0186]k3.1
＝(-0.40,8.0*10^9,-6*10^9,-10^4,-5*10^4,508,-1.25*10^5,-6*10^7,10^4,0,-200)
[0187]
采样率为100msa/s和200ggsa/s，输入信号的幅度0.1v，1v，2v共3组，当采样频率分别是100msa/s和200gsa/s时，仿真信号的频率分别是1mhz，2mhz，2.5mhz，50mhz，10mhz和2ghz，4ghz，5ghz，10ghz，20ghz。仿真结果如图6所示。从仿真结果中可以看出：当输入信号幅度不变时，捏滞迟滞回线的面积随着频率的增加而减小。当输入信号频率不变时，捏滞迟滞回线的面积随着幅度的增加而增加。因此同时增大(减小)输入信号幅度，或减小(增大)输入信号频率，都可以快速增大(减小)捏滞迟滞回线的面积。
[0188]
保持系数矩阵k＝k
2.0
不变，频率分别为500khz，1mhz，2mhz的正弦输入信号在10msa/s，20msa/s，50msa/s下得到的捏滞迟滞回线如图7所示。捏滞迟滞回线实际上由多个直线段拟合的方式实现。分段线性函数数量越多，曲线也就越光滑。比如当过采样率为10(20)时，一个信号周期可量化10(20)个离散数据。得到的捏滞迟滞回线由9(19)个直线段组成。9(19)个直线段也可以由9(19)个分段线性函数的组合。与离散数学模型和连续数学模型方法相比较，本发明提出的基于奈奎斯特采样定理的分段线性函数拟合方法实现更加简单，仅需设置适当的过采样率即可实现捏滞迟滞回线的绘制。图6中的曲线也可以证明上述结论。
[0189]
数值仿真实验4
[0190]
本实验目的在于，设置适当的参数使得提出的忆阻模型能够表现出丰富的局部有源特点。激励正弦电压信号的幅度1v，频率10hz。矩阵k的参数分别设置为k
4.0
＝(0.9,-25,50)、k
4.1
＝(0.9,-35,50)、k
4.2
＝(-0.9,25,50)和k
4.3
＝(-0.9,35,50)。依次绘制的迟滞回线
如图8(a)-图8(d)所示。从图8(a)-(b)(图8(c)-(d))可以看到，增加系数k1的绝对值，忆阻器由无源(有源)忆阻器切换为局部有源忆阻器。同时，迟滞回线的面积也增加。对比图8(a)和图8(c)改变系数k0的极性，可实现无源忆阻器到有源忆阻器的切换。综上所述，改变公式(3)中的系数矩阵k的极性和绝对值可以轻松快速的丰富忆阻器的局部有源特性。证明了第3节理论分析的正确性。
[0191]
数值仿真实验5
[0192]
仿真参数如下：m＝10，ts＝1us，输入信号的幅度和频率分别为0.1v和20khz，分别设置为：
[0193]k5.0
＝(10,10^5,10^5,10^5,7,100,500,0,-10^10,0,0)
[0194]k5.1
＝(10,10^6,10^5,10^5,7,100,500,0,-10^10,0,0)
[0195]k5.2
＝(10,10^5,10^6,10^5,7,100,500,0,-10^10,0,0)
[0196]k5.3
＝(10,10^5,10^5,10^6,7,100,500,0,-10^10,0,0)。
[0197]
和系数k
5.0
相比(如图9(a)所示)，在系数k
5.1
，k
5.2
和k
5.3
中，分别将k1，k2和k3增加10倍得到的捏滞迟滞回线如图9(b)、(c)、(d)所示。从图9(a)-(d)中，我们可以清晰地看见增大系数k1的值可以增大捏滞迟滞回线的面积。而改变相同数量级的系数k2和k3对捏滞迟滞回线的面积几乎不产生影响，这也证明了前面理论分析的正确性。
[0198]
数值仿真实验6
[0199]
仿真参数如下：m＝10。
[0200]
高采样率参数配置如下：f
sh
＝200gsa/s，t
sh
＝5ps，f
ih
＝20ghz，a＝1v，k
5.0
＝(-0.4,8.0*10^9,-6*10^9,-10^4,-5*10^4,508,-1.25*10^5,-6*10^7,10^4,0,-200)。
[0201]
低采样率参数配置如下：f
sl
＝100msa/s，t
sl
＝10ns，f
il
＝1mhz，a＝1v，k
5.1
＝(-0.4,8,-6*10^-9,-10^-23,-4*10^-5,5.08*10^-43,-1.25*10^-49,-6*10^-56,10^-68,0,-2*10^-88)。
[0202]
仿真得到的磁通量、忆导量、捏滞迟滞回线、输入电压和输出电流时域时域波形如图10所示。从图10中我们可以看到，所有仿真数据都相同，除了磁通量。因此，当系数矩阵k满足公式(33)和过采样率满足公式(31)的比例关系时，高频信号形成的捏滞迟滞回线与低频信号形成的捏滞迟滞回线完全一样。
[0203]
通过上述数值仿真实验，我们可以发现本发明提出的忆阻器通用数学模型可以显示出极好的收缩磁滞回线。值得注意的是，改进的忆阻器通过改变不同的系数k即使在更高的激励频率下也可以显示出收缩磁滞回线和出现在坐标区域的任意位置。
[0204]
图11是本发明基于fpga的实时可重构通用忆阻器仿真方法一种具体实施方式流程图。
[0205]
在本实施例中，如图11所示，本发明基于fpga的实时可重构通用忆阻器仿真方法方包括以下步骤：
[0206]
步骤s1：建立忆阻器的数学模型f(h[n])并判断是否是多项式
[0207]
建立忆阻器的数学模型f(h[n])，判断其是否是关于系统状态变量h[n]的多项式，如果不是，进入步骤s2；如果是，则确定数学模型f(h[n])关于系统状态变量h[n]的阶次m，并进入步骤s5。
[0208]
步骤s2：确定系统状态变量h[n]的取值区间
[0209]
根据输入信号的零直流分量交流信号的幅度和频率，分别确定最大幅度a
max
和最小频率ω
min
，并进一步确定系统状态变量h[n]的取值区间：
[0210][0211]
步骤s3：确定多项式阶数m
[0212]
在系统状态变量h[n]的取值区间内，采用麦克劳林公式对数学模型即忆导值或忆阻值f(h[n])进行关于系统状态变量h[n]的m阶多项式拟合，得到拟合函数g(h(n))，并计算最大拟合相对误差εm：
[0213][0214]
设ε0为可接受的最大拟合相对误差，则多项式阶数m的取值应满足εm≤ε0。
[0215]
步骤s4：确定多项式系数ki[0216]
根据麦克劳林公式确定m+1个多项式系数ki,i＝0,1,2,...,m。多项式系数ki可以采用matlab中的taylor()函数快速计算。
[0217]
步骤s5：计算归一化系数α1[0218][0219]
其中，r
ch
为当前模拟通道的量程，且模拟通道的量程需要满足约束条件α1≤1，同时α1的取值范围为α1∈[-1,1]。
[0220]
步骤s6：计算模型输出最大值y
max
[0221]
将h(n)＝a
max
/ω
min
代入：
[0222][0223]
得到f(h[n])的最大值f
max
(h[n])；
[0224]
再将f(h[n])＝f
max
(h[n])和x(n)＝α1代入：
[0225]
y[n]＝f(h[n])
·
x[n]
[0226]
得到模型输出y[n]的最大值y
max
，其中，输入信号x[n]为adc量化得到的有符号的数据。
[0227]
步骤s7：基于fpga进行忆阻器实时仿真即在fpga中进行以下运算：
[0228]
步骤s7.1：计算系统状态变量h[n]
[0229]
如图12所示，adc对模拟的输入信号x(t)量化得到有符号的数据即输入信号x[n]。对于输入信号x[n]，首先使用fpga内部的定点数转浮点数ip核(fp ip)将其转换为单精度浮点数据f_x[n]，其中的取值范围为f_x[n]∈[-1,1]，然后计算系统状态变量h[n]即n时刻的磁通量或电荷量：
[0230][0231]
其中，ts为时间间隔，f_x[j]为输入信号的j个采样点，h[0]为系统状态变量的初始值。
[0232]
步骤s7.2：计算忆导值或忆阻值f(h[n])：
[0233][0234]
步骤s7.3：计算输出信号y[n]
[0235]
如图12所示，将单精度浮点数据f_x[n]同时送入fifo中进行延时处理，使得fifo读端口的数据即延时数据f_dly_x[n]与忆导值或忆阻值f(h[n])在时间上对齐，并计算输出信号y[n]：
[0236]
y[n]＝f(h[n])
·
f_dly_x[n]
[0237]
其中，输入信号x[n]为电压信号或电流信号，输出信号y[n]为电流信号或电压信号。
[0238]
步骤s7.4：计算输出信号div2[n]
[0239]
如图12所示，在fpga中，利用一个除法器ip核，将输出信号y[n]除以最大值y
max
，得到输出信号div2[n]：
[0240][0241]
对延时数据f_dly_x[n]进行延时处理，使得延时后的数据f_dly1_x[n]与div2[n]在时间是对齐的。
[0242]
步骤s7.5：计算用于adc输入的输出信号multi1、multi2
[0243]
如图12所示，在fpga中，利用两个乘法器，将数据f_dly1_x[n]、div2[n]均与相乘，得到输出：
[0244][0245][0246]
其中，n
dac
是用于模拟波形输出dac的分辨率，且dac能够输出双极性波形。
[0247]
步骤s7.6：数模转换输出模拟信号x
′
(t)、y
′
(t)
[0248]
在fpga中，利用dac将输出信号multi1和multi2转换为对应的模拟信号x
′
(t)、y
′
(t)。在本实施例中，输出信号multi1和multi2经过浮点-定点ip转换为定点数据后，再转换为对应的模拟信号x
′
(t)、y
′
(t)。
[0249]
步骤s8：李沙育图形显示
[0250]
将模拟信号x
′
(t)、y
′
(t)输入到数字示波器中进行李沙育图形显示，显示忆阻器的捏滞迟滞回线。
[0251]
功能仿真
[0252]
为了充分评估基于fpga实现所提忆阻器模型即本发明的性能，我们进行了功能仿真实验，并与现有方法进行了对比。
[0253]
verilog硬件编程语言、软件开发环境vivado 2020.1和功能仿真软件xsim被用于所提出忆阻器模型的硬件功能仿真。
[0254]
首先，使用matlab软件生成14位有符号共4096点的数字离散正弦信号。同时在fpga实现了dds原理的波形发送模块。通过频率控制字kf(波形rom读地址的步进值)即可输出频率可调的正弦波。
[0255]
系数矩阵为：
[0256]k1.0
＝(0.4,-3.0*10^6,4.5*10^8,-10^2,1.5*10^6,50.8,-125,6,100,10,-200)。
[0257]
采样间隔ts＝10ns。fpga仿真输出数据为双精度浮点格式以为了方便硬件和软件仿真结果进行对不分析。硬件仿真时序图如图13所示，图13中clk，rst_n和ce分别是系统时钟100mhz，低电平有效的复位信号和模块使能信号。数据dfx和dfy分别表示忆阻器模块的输入电压信号和输出电流信号。dfx_valid(dyx_valid)为高表示输出数据dfx(dfy)。图中dfx_valid和dyx_valid一直为高，这也证明了所提出的忆阻器模块能够在fpga中实现实时运算，且实时运算速度为12.8gbps。
[0258]
fpga硬件仿真输出的双精度浮点数据导入matlab中，与采用matlab软件仿真得到的捏滞迟滞回线如图14(a)所示。从图中可以看到，一个迟滞回路曲线由40个采样率组成，也就是过采样的倍数是40。两种实现方法的误差局部放大如图14(b)所示。fpga仅需400ns可完成一个周期的数据计算。
[0259]
接下来，我们在fpga中线性了一个10次非线性多项式对现有文献提出的6种忆阻器模型进行快速可重构实验，6种忆阻器模型的输入参数和输入双极性正弦信号的频率和幅度(固定1v)如表4所示。
[0260]
我们在fpga中对现有忆阻器模型进行仿真实验。每个模型均输入双极性1v正弦信号。在fpga中固定m＝10，按照如表4所示配置多项式系数、采样率、输入信号频率，得到的输出信号幅度和fpga仿真时间如表4所示。从表4中我们可以发现，提出的实时可重构忆阻器模型输出幅度与输入参数有关。
[0261][0262]
表4
[0263]
在表4中，总结了对不同忆阻器模型和不同输入信号频率运算一个周期fpga硬件实现需要的时间(忽略了固定延时390ns)。在matlab中则模拟了fpga中数据处理流程，具体包括定点数转单精度浮点数，使用单精度浮点数计算磁通量和忆导。与fpga处理相同数据量所需要的时间如表4所示。从表4，我们可以看出，基于本发明实时可重构通用忆阻器仿真电路，不仅可以仿真各种现有忆阻器模型，而且相对于matlab软件仿真，其计算时间大幅度减小。
[0264]
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。