ADC採集噪聲問題及均方根值濾波與Kalman濾波比較
有一陣子筆者在做一個PT100熱電阻的調理電路的時候採用了，使用恆流源的方式測熱電阻的阻值。為了採集方便，將0.3mA的電流接入PT100直接把ADC輸入端接在了PT100的兩端。之後再輸出溫度的時候資料非常亂。查閱資料受到啟發，採用求該訊號的有效值（均方根值）方法進行濾波。
-
ADC取樣、恆流源
ADC使用的是STM32F10x系列的內部12位ADC。恆流源電流值為0.3mA。
-
PT100熱電阻
PT100是在零度下，阻值為100Ω的熱電阻，其溫度特性曲線如圖1
圖1：PT100溫度特性曲線
基本上在-200—600攝氏度之間溫度與電阻是線性關係。故測量PT100電阻值的其中一種辦法就是給PT100接入一個恆流源，然後測其兩端電壓進而求出PT100電阻值。
-
噪聲分析
電路上電（電源採用的是兩節18650串聯，AMS1115-3.3V降壓穩壓）將熱電阻兩端接入示波器（輸入阻抗1MΩ）如圖2。Single一下，然後將介面的所有采樣資料做FFT（快速傅立葉變換）得到圖3、圖4.不難看出噪聲主要來自於1Mhz的頻段和一些高斯白噪聲。
圖2：原始訊號
圖3：FFT
圖4：FFT
-
均方根值去除噪聲
求一段訊號的均方根值其實就是求這段訊號的的有效值（RMS）即訊號平方的均值再開方。ADC取樣頻率為10Khz。每次取100個數據進行計算。計算後資料由串列埠輸出，見圖:6、圖7.（C程式碼見程式1）
圖5：串列埠輸出原始資料
圖6：濾波資料
圖7：給PT100加溫度
-
一階低通濾波器
y(n) = q*x(n) + (1-q)*y(n-1)
中Y（n）為輸出，x（n）為輸入。其中Q取0.5輸入波形見圖8.（C程式碼見程式2）
圖8：加入一階低通濾波串列埠輸出
-
kalman濾波
第一步將資料匯入MATLAB，畫出原始資料曲線。見圖9.。將資料做卡爾曼濾波觀測矩陣。將該資料的t-2到t時刻的均值作為推測矩陣（詳細見附加程式碼）。最後將卡爾曼輸出做滑動平均。卡爾曼輸出模型（藍色）見圖10.（matlab程式碼見程式3）
圖9：MATLAB匯入原始資料
圖10：原始資料與濾波後的資料
-
總結：
噪聲的產生有很大一部分原因應該是在畫電路板的時候沒有把數字地和模擬地分開造成的。Kalman濾波演算法是一種基於時域的濾波演算法。筆者認為該電路中的高斯白噪聲主要來源於電源和數字訊號對模擬訊號的干擾。就本電路而言求訊號的均方根值是比較有效的辦法。筆者學識有限，如有錯誤的地方歡迎指正。
-
程式1
unsigned int disposePT100v()
{
unsigned char i;
double dData[100];
double dSum;
for(i=0;i<100;i++)
{
dData[i] = g_iADC[0];//獲取ADC資料
dSum +=pow(dData[i],2);//做平方和
delay_us(100);
}
dSum /= 100;//求均值
dSum = sqrt(dSum);//開方
return (unsigned int)dSum*3300/4096;//mv電壓值
}-
程式2
unsigned int lowV(unsigned int com)
{
static unsigned int iLastData;
unsigned int iData;
double dPower=0.5;
iData = (com*dPower)+(1-dPower)*iLastData;//計算
iLastData = iData;//jilu
return iData;//返回資料
}-
程式3
clear
clc
ydata=textread('a123.txt','%s');%%匯入資料
cData = hex2dec(ydata);%轉換為10進位制資料
%t-1時刻的輸出作為本次時刻的系統推測數值 最優偏差（初值）g_zer 推測誤差(初值)g_ter 測量值偏差（初值）g_cer
%測量結果(矩陣) ner1 推測結果（矩陣）ner2 (end)
g_zer = 5;%%最優偏差初值
g_ter = 3;%%推測偏差初值
g_cer = 2;%%測量偏差初值
kdata = 0;
for i=1:5792%%迭代5000次
ner1 = cData(i);
if i>3
ner2 = (cData(i-2)+cData(i-1)+cData(i))/3;%%t-2時刻到t時刻的平均值作為t時刻的推測值
else
ner2 = kdata;
end
if g_zer>g_ter
g_er = sqrt((g_zer^2)-(g_ter^2));%計算本次誤差
else
g_er = sqrt((g_ter^2)-(g_zer^2));%計算本次誤差
end
kk = sqrt((g_er^2)/((g_er^2)+(g_cer^2)));%計算卡爾曼增益
kdata = ner1+kk*(ner2-ner1);%更新卡爾曼輸出
g_zer = ((1-kk)*g_ter^2)^0.5;%計算本次最優偏差
cKdata(i) = kdata;
end
%%%%%%%滑動平均%%%%%%%%%%
t = 1:5792;
a = 0;
b = 0;
pKdata = zeros(1,5792);
for i=1:5700
for j=1:50
a = a + cKdata(i+j);
end
b = a/50;%%求平均
a = 0;
pKdata(i)=b;
end
plot(t,cData,'R',t,cKdata,'G',t,pKdata,'B');