基于C#的一維光譜信號(hào)處理方法研究與應(yīng)用

周超　方哲　宋春苗　胡學(xué)強(qiáng)

摘要：光譜信號(hào)的處理是光譜系統(tǒng)的關(guān)鍵部分，包括插值、光滑、尋峰、提取等步驟，插值是根據(jù)已知條件逼近還原真實(shí)信號(hào)，光滑可以消除統(tǒng)計(jì)漲落的影響，尋峰是信號(hào)標(biāo)定的先決條件，信號(hào)提取可以得到每一個(gè)有用信號(hào)單獨(dú)的數(shù)學(xué)表征。這些手段綜合起來(lái)應(yīng)用，可以提高光譜系統(tǒng)的檢出限、精密度、穩(wěn)定性等主要指標(biāo)。本文針對(duì)上述步驟和方法展開(kāi)討論，用C#編程語(yǔ)言實(shí)現(xiàn)其在一維光譜系統(tǒng)中的應(yīng)用，取得了良好的效果。

關(guān)鍵詞：一維光譜信號(hào);信號(hào)處理方法;C#

中圖分類(lèi)號(hào)：TH741;TH842 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1007-9416（2019）05-0043-03

1 一維光譜信號(hào)處理方法

1.1 插值

插值是離散函數(shù)逼近的重要方法，利用它可通過(guò)函數(shù)在有限個(gè)點(diǎn)處的取值狀況，估算出函數(shù)在其他點(diǎn)處的近似值。對(duì)于常見(jiàn)的高斯信號(hào)和高斯疊加多項(xiàng)式的信號(hào)，一般需要在譜圖橫軸平均分布的至少7-9個(gè)點(diǎn)，才能展現(xiàn)信號(hào)的主要特征，在光譜系統(tǒng)中，由于分辨率的問(wèn)題，我們所得到的譜圖的譜峰往往沒(méi)有足夠的數(shù)據(jù)來(lái)表征。這時(shí)候就可以使用插值的方法，補(bǔ)足缺失的點(diǎn)。時(shí)域上，f（x）是定義在[a，b]上的已知映射關(guān)系，x1，x2，x3...xn為區(qū)間內(nèi)n個(gè)互不相同的點(diǎn)，G為給定的某一函數(shù)類(lèi)。若G上有函數(shù)g（x）滿足：g（xi）=f（xi），i=1，2，...n。則稱(chēng)g（x）為f（x）關(guān)于節(jié)點(diǎn)x1，x2，x3...xn在G上的插值函數(shù)。頻域插值是將時(shí)域的一維光譜信號(hào)數(shù)據(jù)首先轉(zhuǎn)換成頻域信號(hào)，對(duì)頻域信號(hào)進(jìn)行低頻插值，一般為補(bǔ)0，然后將頻域信號(hào)再轉(zhuǎn)換成時(shí)域信號(hào)。頻域插值方法的關(guān)鍵代碼如下：

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）;//補(bǔ)0

fftResult = FreqAnalyzer.FFT（fftExpand.ToArray（）， false）.ToList（）;//傅里葉變換

List result = FreqAnalyzer.Cmp2Mdl（fftResult.ToArray（））.ToList（）; //插值結(jié)果

1.2 光滑

在實(shí)際光譜系統(tǒng)中，當(dāng)信號(hào)較弱的時(shí)候，信號(hào)的統(tǒng)計(jì)漲落比較大，往往不是高斯或類(lèi)高斯分布的期望值，有用信號(hào)被淹沒(méi)在統(tǒng)計(jì)漲落之中，從而影響整個(gè)系統(tǒng)的精密度。通過(guò)光滑處理方法可以解決上述問(wèn)題。常用的光滑方法有算術(shù)平均法、多項(xiàng)式最小二乘法、離散函數(shù)褶積變換法、傅里葉變換法等。光滑的目的是通過(guò)數(shù)學(xué)方法，消除信號(hào)中統(tǒng)計(jì)漲落的影響，同時(shí)保留原有信號(hào)的特征信息。其方法的過(guò)程可以概括為以當(dāng)前點(diǎn)為中心，利用其左右n個(gè)點(diǎn)的測(cè)量數(shù)據(jù)，通過(guò)算法修正當(dāng)前點(diǎn)的值。比如，算術(shù)平均法就是用當(dāng)前點(diǎn)及其左右2n個(gè)點(diǎn)的算術(shù)平均數(shù)的值來(lái)替換當(dāng)前點(diǎn)的值。傅里葉變換法則是通過(guò)頻域低通濾波，將統(tǒng)計(jì)漲落看作是高頻的噪聲從原始信號(hào)中消除。頻域?yàn)V波方法的關(guān)鍵代碼如下：

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）;//各項(xiàng)再除以項(xiàng)數(shù)得到結(jié)果數(shù)組

1.3 尋峰

一維光譜信號(hào)的尋峰是信號(hào)標(biāo)定的先決條件。信號(hào)的形狀可用信號(hào)峰位、峰強(qiáng)、半高寬加以描述，必要時(shí)還應(yīng)考慮非對(duì)稱(chēng)因子[1]。首先判定信號(hào)峰是否存在，如果存在就確定峰的位置和邊界，然后才能根據(jù)信號(hào)峰的特征進(jìn)行定性分析。常用的尋峰方法有算術(shù)比較法、導(dǎo)數(shù)法、協(xié)方差法、對(duì)稱(chēng)零面積法等。有信號(hào)重疊現(xiàn)象時(shí)，算術(shù)比較法不能使用。信號(hào)較弱時(shí)，協(xié)方差法的效果會(huì)受到較大影響。導(dǎo)數(shù)法尋峰會(huì)導(dǎo)致一定程度的峰位偏移。對(duì)各種情況綜合比較，對(duì)稱(chēng)零面積法的準(zhǔn)確度較好。通過(guò)窗函數(shù)與一維光譜信號(hào)進(jìn)行褶積變換，當(dāng)變換譜與它的標(biāo)準(zhǔn)偏差之比出現(xiàn)正極值，且大于固定的閾值，即可判定為峰，然后通過(guò)探測(cè)器響應(yīng)函數(shù)或數(shù)學(xué)擬合方法確定峰的邊界。對(duì)稱(chēng)零面積法的關(guān)鍵代碼如下：

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 提取

光譜信號(hào)重疊是光譜自身物理特性和光譜系統(tǒng)分辨率不足的結(jié)果。受制于光路設(shè)計(jì)和探測(cè)器的響應(yīng)曲線，光譜系統(tǒng)不能把兩個(gè)臨近的一維光譜信號(hào)完全分開(kāi)時(shí)，就需要用數(shù)學(xué)方法進(jìn)行提取。一維光譜信號(hào)一般可以認(rèn)為是某種特定分布信號(hào)的疊加，如常見(jiàn)的高斯信號(hào)和高斯加多項(xiàng)式的信號(hào)。提取的目的是通過(guò)數(shù)學(xué)方法，找到用來(lái)表征每個(gè)峰的最優(yōu)函數(shù)，計(jì)算機(jī)處理時(shí)，可以轉(zhuǎn)化為求矩陣的最優(yōu)解[2][3]。常用的方法有最小二乘法、剝譜法、小波變換法等。如果提取后出現(xiàn)殘差較大的情況，需要再次對(duì)所得到的函數(shù)參數(shù)進(jìn)行調(diào)優(yōu)。最小二乘法提取方法的關(guān)鍵代碼如下：

CalCoefArray（）;//計(jì)算系數(shù)矩陣

CalConstArray（）;//計(jì)算常數(shù)矩陣

lEquations.Init（coefArray， constArray）;//coefArray系數(shù)矩陣，constArray常數(shù)矩陣

lEquations.GetRootsetGauss（resultMatrix）;//求方程組的解

for （int i = 0; i < peakNumber; i++）

{

List singlePeak = new 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 實(shí)際應(yīng)用

2.1 單個(gè)一維光譜信號(hào)的一般處理過(guò)程

為了方便敘述信號(hào)處理過(guò)程，本文設(shè)計(jì)了一組復(fù)雜的一維光譜信號(hào)，一共9個(gè)峰，32個(gè)數(shù)據(jù)點(diǎn)。由于數(shù)據(jù)量少，難以進(jìn)行后續(xù)的數(shù)學(xué)運(yùn)算，首先對(duì)原始信號(hào)進(jìn)行8倍的插值處理，如圖1所示。

我們利用上述方法對(duì)插值后的結(jié)果進(jìn)行平滑、尋峰、提取，并通過(guò)對(duì)提取結(jié)果與原始信號(hào)的殘差分析優(yōu)化提取過(guò)程參數(shù)，即可得到每個(gè)峰最終的數(shù)學(xué)表達(dá)，通過(guò)函數(shù)曲線表示，如圖2所示。

2.2 消除信號(hào)統(tǒng)計(jì)漲落的效果分析

圖3是實(shí)際光譜系統(tǒng)在同樣條件下采集的10組信號(hào)，我們對(duì)235-243范圍內(nèi)的數(shù)據(jù)放大后可以看到信號(hào)統(tǒng)計(jì)漲落的影響，我們從每組信號(hào)中選取10個(gè)相同位置的有用信號(hào)，利用上文提到的插值和光滑方法進(jìn)行數(shù)據(jù)處理，結(jié)果如表1所示，可以看到不同組之間相同位置的信號(hào)的相對(duì)標(biāo)準(zhǔn)偏差均有減小，提高了精密度指標(biāo)。

2.3 信號(hào)提取的效果分析

圖4是實(shí)際光譜系統(tǒng)采集的1組信號(hào)，可以看到信號(hào)Sig.1和Sig.2之間有小部分的信號(hào)重疊，信號(hào)Sig.3的一大部分被信號(hào)Sig.4覆蓋，如果不進(jìn)行信號(hào)提取，就無(wú)法對(duì)信號(hào)進(jìn)行下一步的分析和利用。應(yīng)用上文提到的尋峰和提取方法，可以得到每個(gè)信號(hào)的數(shù)學(xué)表達(dá)，并通過(guò)函數(shù)曲線表示在圖上。提取信號(hào)與實(shí)際信號(hào)的相對(duì)偏差如表2所示，相對(duì)大的信號(hào)提取效果好，相對(duì)小的信號(hào)提取效果略差，但是總體上能夠控制在3%以內(nèi)，可以滿足一般的信號(hào)處理要求。

3 結(jié)語(yǔ)

本文介紹了基于C#的一維光譜信號(hào)處理方法與應(yīng)用，介紹了插值、光滑、尋峰、提取的方法，給出了優(yōu)選方法的相關(guān)代碼，并在數(shù)據(jù)逼近還原、消除統(tǒng)計(jì)漲落、信號(hào)尋找、信號(hào)提取等步驟成功應(yīng)用。

參考文獻(xiàn)

[1] 吉昂，陶光儀，卓尚軍，等.X射線熒光光譜分析[M].北京：科學(xué)出版社，2003：251-262.

[2] 王婷婷.發(fā)射光譜儀中光譜干擾校正方法的應(yīng)用和研[D].杭州：杭州電機(jī)科技大學(xué)，2011.

[3] 肖筱南.現(xiàn)代數(shù)值計(jì)算方法[M].北京：北京大學(xué)出版社，2003：186-195.