研究土壤电弧的电离过程。
测量到的电弧电流波形会有直流分量，这个直流分量不知道是系统产生的还是罗氏线圈的零漂，鉴于罗氏线圈是无法测量直流分量的，那么大概率应该是干扰导致的零漂，在分析电弧的电离区域时，需要首先把直流分量去掉。
去掉直流分量并不难，传统经验是用FFT分析波形，将直流量提取出来即可。但电弧是按照半周波为单位完成一次击穿燃弧到熄弧的过程的，所以按照每个周期为单位对直流分量进行识别处理。
首先需要将每个周期的电流波形提取出来，子程序如下。
%%%%%%%%%子函数——数据处理程序，将数据处理成以周期为单位的矩阵%%%%%%%%%
function [U_matrix,Num,Circles] = DataToCircles(U,f)
N = f / 50;%每周期数据点数
n = length(U);%测量向量内共有多少点
Circles = fix(n/N);%取整，周期数
Num = N * Circles;%共有多少点
U = U(1:Num);%取得合适的周期数据，为编程矩阵做初处理
U = U';%转置
U_matrix = reshape(U,N,[]);
end
1. 利用FFT获取直流偏置分量
提取每周期的电流波形后，对单个周波进行FFT分析，FFT子程序如下。
%%%%%%%%%子函数——FFT谐波程序，获取规定次谐波的幅值和相位%%%%%%%%%%%%%
function [Amp_U,Ph_U] = FFT(U_matrix,Circles,L,harmonic)
Y_U_matrix = fft(U_matrix,L);%快速傅里叶变换
P2_U_matrix = abs(Y_U_matrix/L); % 双侧频谱P2
P1_U_matrix = P2_U_matrix(1:L/2 +1,:); % 单侧频谱P1
P1_U_matrix(2:end -1) = 2 * P1_U_matrix(2:end -1);
num = harmonic + 1;
for i = 1:Circles
PY_U_matrix(:,i) = phase(Y_U_matrix(:,i));%频域弧度幅值
PY_U_matrix(:,i) = PY_U_matrix(:,i) * 180 /pi;%频域相位
end
for i = 1 : Circles
Amp_U(i) = P1_U_matrix(num,i);
Ph_U(i) = PY_U_matrix(num,i);
end
end
通过查看波形，可以看出直流分量大约在2.5A左右，可是FFT的分析结果并不准确，很多周波的直流分量远低于2.5A，原因可能是存在大量的高频谐波，这些高频谐波会影响直流分量的幅值，那么只能采用更暴力的处理方法，找到每周波电压为0时的电流值设定为直流偏置分量。
2. 利用电压信号过零点获取直流偏置分量
同样利用DataToCircles子函数将电压和电流均转换为周期数据，建立一个程序子脚本。
%%%%%%%%%子函数——获取电压为零时的电流值，即直流分量%%%%%%%%%
function [I_DC_matrix] = FindCurrentAsVoltageEqualZero(U_matrix,I_matrix,Circles,f)
N = f / 50;%每周期数据点数
for i = 1 : Circles%按照周期进行分析
location = find(U_matrix(:,i) == 0);%找到每周期中的U=0的位置
k = length(location);%共有多少个U=0的点
for j = 1:k%对每一点的位置找到对应的电流值
I_DC_value(j) = I_matrix(location(j),i);
end
I_DC_matrix(i) = mean(I_DC_value);%取平均值
end
end
3. 确定电离区域的起点
大约按照肉眼观察电离区域的起点中的有效数据，即躲避开罗氏线圈的精工电流区间，大概按照电流最大值的10%，认为是有效测量电流，即用10%~90%的电流部分进行模型参数数据分析，其他部分进行递推处理。
%%%%%%%%%子函数——获取半周波电压、电流电离区信号并去除电流直流分量%%%%%%%%%
function [IonizationZoneStartPoint,IonizationZoneEndPoint,IonizationZoneMaxPoint] = FindHalfCircleCurrentIonizationZone(HalfCircleNum,SingleHalfCircleNum,Id,I_DC_matrix)
k = 1;
for i = 1 : HalfCircleNum %分别分析每一个半周波
Id_cal = Id(k : k + SingleHalfCircleNum(i) - 1);%提取第i个半周波电流信号
[Id_cal_max IonizationZoneMaxPoint(i)] = max(abs(Id_cal));%将第i个半周波电流取绝对值后找到最大值及位置
% IonizationZoneMaxPoint(i) = find(Id_cal == Id_cal_max);%找到最大值的位置
for j = 1 : IonizationZoneMaxPoint(i)
if ((Id_cal(j) - I_DC_matrix(i)) / (Id_cal(IonizationZoneMaxPoint(i)) - I_DC_matrix(i))) > 0.1%上升沿电流超过最大值的10%，则记录电离开始点
IonizationZoneStartPoint(i) = j;
break;
else
IonizationZoneStartPoint(i) = nan;
end
end
for j = IonizationZoneMaxPoint(i) + 1 : SingleHalfCircleNum(i)
if ((Id_cal(j) - I_DC_matrix(i)) / (Id_cal(IonizationZoneMaxPoint(i)) - I_DC_matrix(i))) < 0.1%下降沿电流低于最大值的10%，则记录电离结束点
IonizationZoneEndPoint(i) = j;
break;
else
IonizationZoneEndPoint(i) = nan;
end
end
k = k + SingleHalfCircleNum(i);
end
end
观察第一个半周波。
电压，电流和电阻的曲线，可以看到电压的最大值和电流的最大值时刻并不完全相同。V/R特性曲线如图所示。
用一个更直观的方法画出则是下图这个样子的。
