基于C#的一维光谱信号处理方法研究与应用
2019-08-26周超方哲宋春苗胡学强
周超 方哲 宋春苗 胡学强
摘要：光谱信号的处理是光谱系统的关键部分，包括插值、光滑、寻峰、提取等步骤，插值是根据已知条件逼近还原真实信号，光滑可以消除统计涨落的影响，寻峰是信号标定的先决条件，信号提取可以得到每一个有用信号单独的数学表征。这些手段综合起来应用，可以提高光谱系统的检出限、精密度、稳定性等主要指标。本文针对上述步骤和方法展开讨论，用C#编程语言实现其在一维光谱系统中的应用，取得了良好的效果。
关键词：一维光谱信号;信号处理方法;C#
中图分类号：TH741;TH842 文献标识码：A 文章编号：1007-9416（2019）05-0043-03
1 一维光谱信号处理方法
1.1 插值
插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。对于常见的高斯信号和高斯叠加多项式的信号，一般需要在谱图横轴平均分布的至少7-9个点，才能展现信号的主要特征，在光谱系统中，由于分辨率的问题，我们所得到的谱图的谱峰往往没有足够的数据来表征。这时候就可以使用插值的方法，补足缺失的点。时域上，f（x）是定义在[a，b]上的已知映射关系，x1，x2，x3...xn为区间内n个互不相同的点，G为给定的某一函数类。若G上有函数g（x）满足：g（xi）=f（xi），i=1，2，...n。则称g（x）为f（x）关于节点x1，x2，x3...xn在G上的插值函数。频域插值是将时域的一维光谱信号数据首先转换成频域信号，对频域信号进行低频插值，一般为补0，然后将频域信号再转换成时域信号。频域插值方法的关键代码如下：
ifftResult = FreqAnalyzer.FFT（data， true）.ToList（）;//傅里叶逆变换
for （int i = 0; i < data.Count（）; i++）
{
ifftResult[i].Real = ifftResult[i].Real / data.Count（）;
ifftResult[i].Image = ifftResult[i].Image / data.Count（）;
}
ListShift（ifftResult）;//序列前后两段换位
fftExpand = ZeroFilling（ifftResult）;//补0
fftResult = FreqAnalyzer.FFT（fftExpand.ToArray（）， false）.ToList（）;//傅里叶变换
List
1.2 光滑
在实际光谱系统中，当信号较弱的时候，信号的统计涨落比较大，往往不是高斯或类高斯分布的期望值，有用信号被淹没在统计涨落之中，从而影响整个系统的精密度。通过光滑处理方法可以解决上述问题。常用的光滑方法有算术平均法、多项式最小二乘法、离散函数褶积变换法、傅里叶变换法等。光滑的目的是通过数学方法，消除信号中统计涨落的影响，同时保留原有信号的特征信息。其方法的过程可以概括为以当前点为中心，利用其左右n个点的测量数据，通过算法修正当前点的值。比如，算术平均法就是用当前点及其左右2n个点的算术平均数的值来替换当前点的值。傅里叶变换法则是通过频域低通滤波，将统计涨落看作是高频的噪声从原始信号中消除。频域滤波方法的关键代码如下：
fftResult = FreqAnalyzer.FFT（data，false）;//傅里叶变换
ListShift（fftResult）;//序列前后两段换位
for （int i = 0; i < data.Length; i++）
{
fftResult[i].Real = fftResult[i].Real * Hd[i];
fftResult[i].Image = fftResult[i].Image * Hd[i];
}
fftResult = FreqAnalyzer.FFT（fftResult， true）;//傅里叶逆变换
data = FreqAnalyzer.Cmp2Mdl（fftResult）;//各项再除以项数得到结果数组
1.3 寻峰
一维光谱信号的寻峰是信号标定的先决条件。信号的形状可用信号峰位、峰强、半高宽加以描述，必要时还应考虑非对称因子[1]。首先判定信号峰是否存在，如果存在就确定峰的位置和边界，然后才能根据信号峰的特征进行定性分析。常用的尋峰方法有算术比较法、导数法、协方差法、对称零面积法等。有信号重叠现象时，算术比较法不能使用。信号较弱时，协方差法的效果会受到较大影响。导数法寻峰会导致一定程度的峰位偏移。对各种情况综合比较，对称零面积法的准确度较好。通过窗函数与一维光谱信号进行褶积变换，当变换谱与它的标准偏差之比出现正极值，且大于固定的阈值，即可判定为峰，然后通过探测器响应函数或数学拟合方法确定峰的边界。对称零面积法的关键代码如下：
int m = （int）（hFWHM + 1）;//hFWHM半峰宽
if （i - m > 0）
{
W = 2 * m + 1; //窗宽
double[] G = new double[W];
double[] C = new double[W];
for （int j = 0; j < W; j++）{G[j] = Math.Cos（Math.PI * （j-m） / 2 / hFWHM）;}
double d = G.Sum（） / W;
for （int j = 0; j < W; j++）{C[j] = G[j] - d;}
for （int j = 0; j < W; j++）
{
y[i] += C[j] * netSpectrum[i - m + j];
dy[i] += C[j] * C[j] * netSpectrum[i - m + j];
}
dy[i] = Math.Sqrt（dy[i]）;
SSi[i] = y[i] / dy[i];//SSi[i]和阈值比较定峰
}
1.4 提取
光谱信号重叠是光谱自身物理特性和光谱系统分辨率不足的结果。受制于光路设计和探测器的响应曲线，光谱系统不能把两个临近的一维光谱信号完全分开时，就需要用数学方法进行提取。一维光谱信号一般可以认为是某种特定分布信号的叠加，如常见的高斯信号和高斯加多项式的信号。提取的目的是通过数学方法，找到用来表征每个峰的最优函数，计算机处理时，可以转化为求矩阵的最优解[2][3]。常用的方法有最小二乘法、剥谱法、小波变换法等。如果提取后出现残差较大的情况，需要再次对所得到的函数参数进行调优。最小二乘法提取方法的关键代码如下：
CalCoefArray（）;//计算系数矩阵
CalConstArray（）;//计算常数矩阵
lEquations.Init（coefArray， constArray）;//coefArray系数矩阵，constArray常数矩阵
lEquations.GetRootsetGauss（resultMatrix）;//求方程组的解
for （int i = 0; i < peakNumber; i++）
{
List
for （int j = 0; j < spectrum.Count（）; j++）
{
singlePeak.Add（resultMatrix[i， 0] * Math.Exp（-Math.Pow（（positionList[j] - peakList[i].position） / peakList[i].peakWidth， 2）））;
}
peakList[i].singlePeakSpectrum = singlePeak;
}
2 實际应用
2.1 单个一维光谱信号的一般处理过程
为了方便叙述信号处理过程，本文设计了一组复杂的一维光谱信号，一共9个峰，32个数据点。由于数据量少，难以进行后续的数学运算，首先对原始信号进行8倍的插值处理，如图1所示。
我们利用上述方法对插值后的结果进行平滑、寻峰、提取，并通过对提取结果与原始信号的残差分析优化提取过程参数，即可得到每个峰最终的数学表达，通过函数曲线表示，如图2所示。
2.2 消除信号统计涨落的效果分析
图3是实际光谱系统在同样条件下采集的10组信号，我们对235-243范围内的数据放大后可以看到信号统计涨落的影响，我们从每组信号中选取10个相同位置的有用信号，利用上文提到的插值和光滑方法进行数据处理，结果如表1所示，可以看到不同组之间相同位置的信号的相对标准偏差均有减小，提高了精密度指标。
2.3 信号提取的效果分析
图4是实际光谱系统采集的1组信号，可以看到信号Sig.1和Sig.2之间有小部分的信号重叠，信号Sig.3的一大部分被信号Sig.4覆盖，如果不进行信号提取，就无法对信号进行下一步的分析和利用。应用上文提到的寻峰和提取方法，可以得到每个信号的数学表达，并通过函数曲线表示在图上。提取信号与实际信号的相对偏差如表2所示，相对大的信号提取效果好，相对小的信号提取效果略差，但是总体上能够控制在3%以内，可以满足一般的信号处理要求。
3 结语
本文介绍了基于C#的一维光谱信号处理方法与应用，介绍了插值、光滑、寻峰、提取的方法，给出了优选方法的相关代码，并在数据逼近还原、消除统计涨落、信号寻找、信号提取等步骤成功应用。
参考文献
[1] 吉昂，陶光仪，卓尚军，等.X射线荧光光谱分析[M].北京：科学出版社，2003：251-262.
[2] 王婷婷.发射光谱仪中光谱干扰校正方法的应用和研[D].杭州：杭州电机科技大学，2011.
[3] 肖筱南.现代数值计算方法[M].北京：北京大学出版社，2003：186-195.
