磁滞回线是分析磁性薄膜性能的重要依据之一。无论通过何种方式采集信号，获取的数据点都是离散的。因此，如果要相对精确地分析磁场强度和磁化强度的关系，必须对这些离散点做拟合处理。

本文不再赘述磁滞回线产生的原理以及测量装置的搭建等。简单画了个图，用到的是纵向克尔效应。

一句话概括就是用线偏振光（激光经过起偏器后成为线偏振光）照射样品，同时对样品施加一个正弦变化的磁场，收集透射光或者反射光（透射用法拉第效应，反射用克尔效应）。经过一系列处理后，我们将得到一个矩形信号，下面都用克尔信号代替。

磁滞回线便是以磁场信号为横坐标，克尔信号为纵坐标绘制而成。如下图所示。

以上内容随便翻下有关磁学的书应该都能找到。

言归正传，下面是关于磁滞回线的代码及相应的解释说明。

原始数据：实验采用高速数据采集卡，磁场信号和克尔信号分别被数据采集卡的两个通道同时记录。数据被保存至txt文件中。

txt文件中包含720000行2列的原始数据，每个通道记录720000个点。第一列记录了磁场信号，第二列记录了克尔信号，两列数据之间用空格间隔。

import numpy as np  

filename = r'C:\Users\Dell\Desktop\测量原始.txt'  # ’’内为原始数据的路径

data = np.genfromtxt(filename,delimiter='')  #两列数据之间用空格间隔，故delimiter=''，若两列数据之间用逗号间隔，则用delimiter='，'。用np.genfromtxt()函数的原因是因为数据量过大时，直接用列表速度会很慢。

a0 = data[:,0] #data[:,0]表示从头取到尾，0为第一列数据，即磁场信号

a1 = data[:,1] #1为克尔信号

#如果将720000个点直接用曲线连接，由于数据点过多且存在异常值，绘制的曲线相当粗糙，且在峰值处过于密集。因此将每100个数据平均为1个点，如果数据点在平均值加减3倍标准差范围外，则判定为异常值，用平均值代替该值。在下列替代异常值并平顺光滑数据的处理中，与x相关的数据都为磁场数据数据，与y相关的都为克尔信号数据。

#磁场信号处理：用平均值代替异常值，每100个点取平均

x=[]

for i in range(0, len(a0), 100):

    x1=a0[i:i+100]

    x2=np.mean(x1)  #平均数

    x3=np.std(x1)   #标准差

    for i in range(i,i+100,1):

        if math.fabs(a0[i]-x2)>3*x3:

            a0[i]=x2   #用平均值代替异常值

        x2=np.mean(x1)

    x.append(x2)

#上图中，蓝色曲线为原始磁场信号，橙色为平均后的磁场信号。平均后的磁场信号相对比较光滑，但在峰值处起伏仍然明显。

#克尔信号处理

y=[]

for i in range(0, len(a1), 100):

    y1=a1[i:i+100]

    y2=np.mean(y1)  #平均数

    y3=np.std(y1)  #标准差

    for i in range(i,i+100,1):

        if math.fabs(a1[i]-x2)>3*x3:

            a1[i]=y2  #用平均值代替异常值

            y2=np.mean(y1)

    y.append(y2）

#拟合正弦磁场函数。A, k, theta,c分别为正弦函数的振幅，波矢、初始相位、直流分量c四个参数。人们已尝试对其中三个参数进行拟合，发现拟合过程是线性闭合的，且绝对收敛。而四参数正弦曲线拟合由于是非线性迭代拟合过程，如果初始值与真实值差别较大，拟合结果中的频率参数可能无法收敛于实际的频率或收敛到局部最优点而不是全局最小值。如果使用optimize.basinhopping（）全局优化函数则可以找到正弦函数的四个参数。其前两个参数为目标函数和初始值，niter为全局优化算法的迭代次数，迭代次数越多，越有可能找到全局最优解，但同时也增加了运行时间。我们需要根据实验得到的数据，选取恰当的迭代次数。为了保守起见，代码中niter取值为50。因为optimize.basinhopping（）是全局优化函数，所以初始值的选择对优化过程并不怎么重要。代码中将四个参数的初始值都赋值为1（参考：张若愚.Python科学计算（第2版）.北京:清华大学出版社.2016:128）。method可采用Nelder-Mead、Powell、CG、BFGS、Newton-CG、L-BFGS-B等，它们的区别在于搜索路径。

def func(x,p):

    A, k, theta,c=p

    return A*np.sin(k*x+theta)+c

def error(p,y,x):

    return np.sum((y-func(x,p))**2)

time0= np.arange(0, len(x), 1) #建立跟磁场信号长度一致的数组，初始值为0，终值为len(x)-1,即7199(python中，取不到终值len(x)，即7200），间隔为1,共计7200个数值

result0 = optimize.basinhopping(error,[1, 1 ,1 ,1],niter=50,minimizer_kwargs=

                    {"method":"Powell","args":(x,time0)})

f=func(time0,result0.x)#以time0为横坐标，拟合后的磁场信号为纵坐标，建立一个函数

f0=f(max(f)+min(f))/2#减去直流分量，即将磁场信号沿y坐标平移至关于x轴对称。

#由图可以看出，磁场信号（未减去直流分量）的拟合相当成功，而且弥补了原始磁场信号由于跟不上在峰值处较平的缺陷。

#在尝试对克尔信号做相似的拟合处理时发现，拟合后的克尔信号的峰值与实际有所偏差。而且每一次拟合，峰值并不相同，增大niter值拟合效果也并未改善，同时还增加了运算时间。

#拟合克尔信号代码如下：

time1= np.arange(0, len(y), 1)

def func1(x,p):#构建方波克尔函数

    A, f, phi,c=p

    y = np.where(np.mod(x-phi,1/f)<1/(2*f), A+c, 0)

    y = np.where(np.mod(x-phi,1/f)>1/(2*f), -A+c, y)

    return y

def error1(p,y,x):

    return np.sum((y-func1(x,p))**2)

result = optimize.basinhopping(error1,[1, 1 ,1 ,1],niter=50,minimizer_kwargs={"method":"Powell","args":(y,time1)})

f1=func1(time0,result.x)

#对于克尔信号的拟合不理想的原因，个人认为可能是因为构造的矩形函数有问题，也有可能是因为克尔信号的漂移背景对拟合造成了干扰，恳请指正。

#采用样条插值，得到克尔信号的漂移背景，将其从原始数据中分离出去。由于上面的数据周期太小，且漂移背景起伏不明显，所以下面的原始数据增加了周期数。

#采用样条插值，得到克尔信号的漂移背景，将其从原始数据中分离出去。

maximum =[]#正弦极大值

for i in range(1,len(x)-1):

    if float(f0[i-1])<float(f0[i]) and float(f0[i+1])<float(f0[i]):

        maximum .append(i)

minimum=[]#正弦极小值

for i in range(1,len(x)-1):

    if float(f0[i-1])>float(f0[i]) and float(f0[i+1])>float(f0[i]):

        minimum.append(i)

t=maximum+minimum#正弦周期分段

t0=sorted(t)#时间从小到大

t1=[]

y1=[]

for i in range(0,len(minimum)-1): 

    T=int((minimum[i+1]+minimum[i])/2)

    Y=sum(y[minimum[i]:minimum[i+1]])/len(y[minimum[i]:minimum[i+1]])

    t1.append(T)

    y1.append(Y)

#添加首尾数据

t1.insert(0,0)

t1.append(len(y)-1)

y1.insert(0,sum(y[minimum[0]:minimum[1]])/len(y[minimum[0]:minimum[1]]))

y1.append(Y)

#三阶插值

xnew = np.arange(0, len(y), 1) 

func1 = interp1d(t1, y1, kind='cubic')

ynew1 = func1(xnew)

Y=y-ynew1#除去漂移后的克尔效应值,中心在y=0

突然记起来作业还没写，先撤了