滑動窗口濾波c語言實現,滑動窗口濾波器

本文目錄一覽：

• 1、　隨機雜訊的濾波
• 2、　實現方法
• 3、AvaReader如何採用濾波,過濾掉雜訊
• 4、平滑濾波的濾波方法
• 5、Gabor濾波器
• 6、plc模擬量輸入濾波程序和方案？

　隨機雜訊的濾波

在觀測的結果中，不可避免地會有隨機雜訊。產生雜訊的原因有三方面，第一是儀器本身讀數不穩定引起的儀器雜訊；第二是觀測者肉眼的分辨力不夠所引起的觀測誤差，這種誤差具有隨機性質；第三是儀器所處的環境引起的干擾，這裡所指的環境包括自然因素和人文因素，干擾則可分為隨機干擾和非隨機干擾，後者往往具有局部場的特點，而隨機干擾則是這裡所指的隨機雜訊。

不失一般性，假定觀測結果是單變數x的函數，x可以是位置變數，也可以是時間變數。設觀測值用t1（x）表示，區域場用t0（x）表示，雜訊用n（x）表示，則有

資訊理論、系統論與地質找礦工作

假定n（x）具有分段均值為零的性質，即

資訊理論、系統論與地質找礦工作

式中x0為計算平均值的點的坐標，△x為採樣點點距，（2n+1）△x為計算平均值的窗口。

選擇削弱雜訊濾波方法的原則是：依據區域場的性質，選擇方法，使濾波所引起的區域場的畸變儘可能地小，而雜訊的削弱程度則儘可能地大。因此，在（7—2）式成立條件下，根據區域場的特點，常用來削弱雜訊的濾波方法有以下五種：

第一種方法，滑動窗口平均法。這種方法適用於區域場在滑動窗口內呈線性變化的情況，即

資訊理論、系統論與地質找礦工作

將坐標原點取在窗口中心，則有xi=i△x及

資訊理論、系統論與地質找礦工作

對（7—1）式在（-n,n）窗口範圍內取平均，並令平均值為A，則有

資訊理論、系統論與地質找礦工作

資訊理論、系統論與地質找礦工作

根據（7—2）及（7—4）式，上式化簡為：

A=a=t0（0）

即平均值即為在窗口中心的區域異常值。圖7—2是一個例子。圖中區域場用下式表示：

資訊理論、系統論與地質找礦工作

式中x=i△x；△x=1;i=1,2,3，…，100。

圖7—2　疊加正負等幅雜訊的信號

1—信號；2—疊加雜訊後的信號；3—雜訊

隨機雜訊為正負等幅，振幅為10，加在8到90號點上，濾波窗口取11△x。由於用3單數點求平均值，故（7—2）式不成立，而有雜訊的平均值為A：

資訊理論、系統論與地質找礦工作

如果對平均的結果再取一次平均，設其均值為A1，則有

資訊理論、系統論與地質找礦工作

圖7—3是對圖7—2中的t（x）兩次取平均的結果。從圖看出，消除隨機雜訊的效果是很好的。但仔細看圖，還可以發現平均的結果在實際值的上下有小的波動。為了圖面清晰，濾波結果每隔5個點取1個點在圖上表示，而區域場則用曲線表示。

第二種方法，自動趨勢面擬合法。這種方法適用區域場可用變數x的低價多項式近似時的情況。例如對（7—5）式，令y=x-50，則有

t0（y）=260000/（10000+y2）=260000（10000+y2）-1

由於-50≤y≤50，故y2/10000的最大值為2500/10000=0.25。因此，在（7—5）式中x的取值範圍內，令x1=y2/1000,x1＜1，故可將上式展開為

t0（y）=260000（10000+y2）-1=26（1+y2/10000）-1

資訊理論、系統論與地質找礦工作

已知

的最大值為0.25×0.25=0.0625，故上述展開式可取到x1的1次項，即y的平方項。

圖7—3　圖7—2曲線2兩次滑動窗口平均的結果（窗口大小為11）

1—信號；2—平均法的結果

圖7—4圖　7—2曲線2用自動趨勢面擬合的結果

1—信號；2—自動趨勢面擬合的結果

圖7—5　圖7—2曲線2用高頻截斷法濾波結果

1—信號；2—高頻截斷濾波結果

因此，對圖7—2中的t（x），也可用自動趨勢面擬合法直接求出區域場。圖7—4是用這種方法對t（x）的擬合結果，計算程序自動確定擬合多項式的階次為2。從圖看出，計算的結果是令人滿意的。

第三種方法，頻率域的高頻截斷法。由於區域異常以低頻成分為主，而隨機雜訊則以高頻為主，故高頻截斷能很好地消除信號中的雜訊。圖7—5是用這種方法對圖7—2中t（x）的濾波結果。從圖可以看出，濾波結果是好的。

高頻截斷所用的濾波運算元或頻率域中的響應函數L（ω）為

資訊理論、系統論與地質找礦工作

式中△ω＝2π/[dx×（2N+1）];dx＝取樣點距；N＝採樣點數的1/2（採樣點數應為偶數）；k為濾波器參數，應根據t（x）的頻譜特點選定。此次濾波選k=0.5N，由於採樣點數為100，即N=50，故有k=0.5N=25。

第四種方法，斯特拉霍夫濾波器[32]。這種濾波器應用的條件是：雜訊n（x）應是各態遍歷的、平均值為零的隨機函數；

，h為區域場源的埋深，t0（ω）為T0（x）的傅里葉變換。圖7—6是用這個方法對圖7—2中t1（x）的濾波結果。濾波所用公式是：

資訊理論、系統論與地質找礦工作

式中C0=0.4452,C±1=0.0944,C±2=0.0932,C±3=0.0911.

上述係數適合於△x/h=0.1的情況。由於此處區域場源的埋深h=100，即△x/h=0.01，故濾波效果不如前述三種方法的濾波效果好。

圖7—6　圖7—2曲線2用斯特拉霍夫濾波器的濾波結果

1—信號；2—斯特拉霍夫濾波器的濾波結果

第五種方法，向上延拓。從數學上講，向上延拓高度h相當於對t1（x）的頻譜T1（ω）乘一個因子exp（-hw）。因此，以高頻成分為主的雜訊衰減比以低頻成分為主的區域異常衰減快得多。圖7—7是用這種方法對圖7—2中t1（x）的濾波結果，上延高度h=4△x。從圖看出，濾波結果是好的。這是因為區域場源埋深為100△x，故上延4△x高度將使區域極大值下降3.8%，影響不大。對雜訊而言，從頻譜圖7—8可以看出，t1（x）中雜訊的譜主要集中於ω＞π/2△x部分，而當ω＝π/2△x時，exp（-4△xπ/2△x）=exp（-2π/△x）＝exp（-2π）＝1.87×10-3。即雜訊的衰減將大於99.8%。

圖7—7　圖7—2曲線2向上延拓的結果（延拓高度為4倍點距）

1—信號；2—向上延拓高度4△x的結果（△x為採樣點距）

圖7—8　功率譜曲線圖

因此，用高頻截斷法或向上延拓法削弱信號中的雜訊時，最好先作頻譜分析，計算出t1（x）的頻譜圖，了解雜訊譜分布的主要區間，據此確定（7—6）式中的k值或向上延拓高度。

　實現方法

4.2.2.1　M周期的譜分析

1.M周期譜分析步驟

M周期的識別與劃分用譜分析來實現，具體步驟包括：

①對測井曲線濾波，目的是了解測井曲線的頻譜組成和消除各種干擾。

②消除測井曲線的整體漂移，減小長周期信號的影響。

③求取頻譜，可用付氏變換直接求取頻譜，也可用Burg最大熵譜法求取頻譜。

④求取M周期，劃分沉積旋迴。根據已知的經驗，我們可以確定出塔北各時期的最大沉積速率。已知M周期分別為413ka、125ka、95ka和21ka，這樣，就可確定各級M周期所對應的最大地層厚度，從而可確定分析時窗的大小（見圖4—6），並可保證選取的時窗大於一個M周期。

圖4—6　分析時窗選擇圖解

具體作法是：在橫軸上找到已知的沉積速率，投映到待研究的周期的對應的直線上，作縱軸的垂線，讀出其厚度值，即為所需窗長。所選擇的窗長不應小於該厚度。假定沉積速率為每千年500mm，若要識別413ka的周期，對應的厚度是206.5m，窗長至少要達到這一厚度，而要識別125ka的周期則需要62.5m的窗長，而濾波分析時滑動窗口的步長的選擇則很靈活，根據研究的詳細程度選2m、5m、10m、20m等。

如果測井曲線中某個沉積旋迴正好全部處於分析時窗內，則功率譜上就會顯示幾個不同波長的譜峰，雖然由於沉積速率不同，譜峰的波長是會變化的，但譜峰的頻率的倒數（即其所代表的厚度）的比值，接近於M周期的比，即410:123:41。如果分析時窗跨越兩個沉積旋迴，則譜峰的數目與波長就會失去規律性。據此，我們可以識別沉積旋迴和劃分旋迴邊界。

測井曲線的譜分析，除了用於確認M周期外，還可用於沉積速率分析，計算沉積速率。

2.M周期的確認

識別M周期可用譜峰比值法或譜峰強度。

①譜峰比值：在圖4—7中諧波的波長比在5241.7m深度處，第一、第三、第五個譜的頻率比為0.23∶0.68∶1.8，其倒數比（厚度比）接近410∶123∶41，即接近M周期的比值。表明在以此深度點為中心的分析窗口內，包含著一個完整的M周期。

圖4—7　Milankovitch周期的譜分析

圖4—7中深度5293.7m處的第一、第二、第五個大的譜峰的比為0.1:0.3:0.9，其倒數也接近410∶123∶41。

圖4—7中深度5189.7m處的第一、第三、第五個譜大的譜峰的比值為0.19∶0.58:1.79，倒數的比值接近M周期410:123:41。

功率譜上極大譜峰的比接近M周期的比是M周期存在著強有力的證據。如有古地磁資料、同位素測年資料、微體古生物定年資料，亦可作為另外的證據。

②譜峰的強度：圖4—7中對應譜峰最強的是410ka的M周期，這與大的周期對沉積的控制作用有關，小周期疊加在大周期的變化作用之上，它們的比例分別對應其周期長度。

3.沉積速率變化的譜嶂顯示

沉積速率的漸變：在譜峰上顯示為峰值向右（沉積速率變小）或向左（沉積速度變大）發生小的移動。如圖4—7中5245.7m至5247.7m的第一個最大峰由0.35移到0.19（周/m），並且下一深度還在向左偏移，說明沉積速度向上減小（向下增大）。據沙11井的綜合柱狀圖知道，該井段為淺灰色粉砂質泥岩夾綠灰色泥岩，在51.2m的時窗內頂部偏泥而底部偏砂（在自然伽馬曲線上的顯示更明顯），反應了向上沉積速度減小。

沉積速率突變：5177.7m至5179.7m上下頻譜特徵的突變使我們了解到沉積速率的突變。該深度上下的譜特徵有明顯的差別，必然是沉積速率變化的結果.在岩性柱狀圖上，頂部為綠灰色粉砂質泥岩與泥質粉砂岩互層，而其下有一層綠灰色長石岩屑細砂岩。反應了向上岩性突變，沉積速率突變。

5221.7m上下的譜峰突變也同樣反應了沉積速率突變（見圖4—7）。

4.計算沉積速率

雖然在幾萬年的時間裡沉積速率不可能保持不變，對應的垂向周期的波長也將有變化，但是由於M周期的時間比恆定，這些周期在窗口內的比也保持恆定，窗長內存在M周期則說明此段地層是在一個完整周期內沉積下來的，這樣就不難求得沉積速率。各代表點的沉積速率標在圖4—7中。計算公式如下：

V（H）＝1/（F*T）

式中，V（H）表示H深度的沉積速率、F表示由頻譜圖上讀得的周期為T的諧波的頻率（單位：周/m，相當於該周期波波長的倒數）。

以下幾例分析的結果說明，利用M周期求取沉積速率是有實際意義的。

在以5221.7m為中心的長度為51.2m時窗內，頂部為淺灰白色細砂岩，中間為大段的淺綠灰色泥質粉砂岩夾一層岩屑細砂岩，說明該段為一較快速的沉積段，計算得沉積速率為48.78mm/ka。

在以5229.7m為中心的長度為51.2m的時窗內，是大段的綠灰色泥岩與淺灰色泥質粉砂岩，泥岩發育很好，表示了一種慢速的沉積過程，計算得沉積速率為7.31mm/ka。從5229.7m到5221.7m岩性變粗，由綠灰色泥岩過渡到粉砂質泥岩夾細砂岩，是向上水深變淺、粒度變粗的過程，沉積速率也由小變大。

在以5241.7m為中心的長度為51.2m時窗內，為綠灰色泥岩與粉砂質泥岩，沉積仍為一慢速過程，計算得沉積速率為10.6mm/ka。

在以5247.7m為中心的長度為51.2m時窗內，仍為綠灰色泥岩與粉砂質泥岩，5250m上下自然伽馬對應兩個降低段，是粉砂質泥岩含砂較多的表現，在該時窗及上下時窗範圍內含粉砂較多，計算得沉積速率為16.3mm/ka，大於上下時窗內的沉積速率。

以5255.7m為中心的長度為51.2m時窗內，為泥岩與粉砂質泥岩，沉積較慢，計算得的沉積速率為13.6mm/ka。

以5267.7m為中心的長度為51.2m時窗內，為大段的泥岩與粉砂質泥岩，計算得沉積速率為16.6mm/ka。

Milankovitch周期是地層中普遍存在的，並且是在測井解析度意義下可體現的周期信號，是利用測井分析地層周期性、進行層序劃分與可容納空間分析的基礎。

4.2.2.2　Fischer圖解應用條件修正

Fischer圖是根據潮坪淺水環境中形成的碳酸鹽岩層序來求取海平面變化的一種圖解。其理論依據是：在潮坪淺水環境中，碳酸鹽生長速率很高，在任何時候，由相對海平面上升所增加的可容納空間，隨時都能被新產生的碳酸鹽物質所充填。因此，在相同周期的海平面變化過程中所形成的碳酸鹽岩旋迴的厚度，間接地反映了該周期相對海平面升降的幅度。這一基本假設前提在硅質碎屑環境中則不滿足，在深水中，往往沉積地層厚度小，淺水環境中沉積地層厚度大，這就需要作水深校正，使該旋迴的厚度經校正後能反映可容納空間的大小。

根據統計分析，自然伽馬、自然電位和視電阻率等曲線與岩石平均粒度有較好的對應關係（如圖4—8,4—9，），黃智輝等人在不同地區的研究也證實了這種規律，如圖4—10所示。對應的公式如下：

ΔGR＝（GR-GRmin）/（GRmax-GRmin）

ΔSP＝（SP-SPmin）/（SPmax-SPmin）

MD=Aexp（-CΔGR）

式中，GR、SP分別為自然伽馬值和自然電位值，MD表示平均粒度，A.C為常數。

圖4—8　自然伽馬（GR）與岩性的對應關係（據沙32井）

在一些特別的情況下可能有不同的例子，但從總體來說，自然伽馬電位與平均粒度有較好的對應關係。

在靜水條件下，隨可容納空間的增大，水深增加，距離岸線也將更遠，沉積物的供應必將減少，在相同的時間段內沉積的地層厚度也將減小，粒度將變細，自然伽馬曲線值將增大。可容納空間越大，水深越深，離岸線越遠，沉積物的旋迴厚度就越小。反之，旋迴厚度就變大，而粒度就將變粗，自然伽馬曲線也將減小。

圖4—9　泥質含量與地層電阻率交會圖

圖4—10　自然伽馬曲線與平均粒度對應關係（據黃智輝，1986）

這樣，當一個旋迴的可容納空間大時，應對應較大的沉積物厚度，而實際情況是旋迴厚度較小。而可容納空間變小時，應對應的沉積厚度較小，而實際情況是旋迴厚度大。大的沉積可容納空間對應細粒的沉積物和較小的旋迴厚度，而小的沉積可容納空間對應粗粒的沉積物和較厚的旋迴厚度。我們設計了一個指數模型來逼近這種變化，因為指數模型能展開為多項式形式，是眾多的逼近函數的首選函數。校正的公式如下：

H′＝Hexp（a·Δ）

式中H′是校正後的旋迴厚度，H是校正前的旋迴厚度，a是校正的常數，△可表示為：

Δ＝MIN（ΔGR，ΔSP）

式中MIN是取極小的函數，由於各種因素對測井曲線的影響使得相應的測井曲線值升高，取極小可以減小偶發因素的影響。

參數a的選擇要謹慎，要依據井區的環境資料或根據經驗，目的是儘可能地補償水深的影響，而又避免引進不合理的變化因素。

圖4—11是校正前後的圖解曲線，由於旋迴多，僅作出包絡線，左邊第一欄是自然伽馬曲線，第二欄是校正前的可容納空間變化曲線，第三欄是校正後的可容納空間變化曲線，校正後對應於較好的泥岩段可容納空間變化增大更明顯，更能反映相對海平面的變化，而未校正的曲線中間部分的下降顯然不如校正後中—上部的一致偏高更接近於實際。

圖4—11　校正前後的可容納空間變化曲線

4.2.2.3　區域場與局部場的分離

上述處理得到的曲線含有各種變化成份，而且長波信號的幅度遠大於短波信號，變化快的低幅信號顯示不清或根本顯示不出來。從圖4—11中也可以清楚地看到這一點。為了提高分析信號的能力，採用了提取區域場與局部場的方法。

設待處理信號為Y（X），有下列實現模型：

新疆塔里木盆地北部應用層序地層學

式中，Ai是係數，N是多項式的次數，已有的序列作為它的解，用最小二乘法求解超定方程，得到係數，就可用上述公式計算區域場和局部場。

Z（X）＝Y′（X）-Y（X）

式中Z（X）表示X點的局部場值，Y′（X）表示X點的實際信號值，Y（X）是X點區域場值。

圖4—12是經過各階多項式提取後的區域場與局部場，右邊第一欄為原始曲線，左邊一、二、三、四欄分別是二、三、四、五次多項式代表的區域場與局部場。從圖中可以清楚地看到，各短波信號明白無遺。而未處理前的信號由於長波信號太強使短波顯示不出來。這種方法實現了在同一個曲線上顯示各信號的能力。

圖4—12　區域場與局部分離

AvaReader如何採用濾波,過濾掉雜訊

1、給定均值濾波窗口長度，對窗口內數據求均值。

2、作為窗口中心點的數據的值，之後窗口向後滑動1，相鄰窗口之間有重疊。

3、邊界值不做處理，即兩端wid_length2長度的數據使用原始數據即可過濾掉雜訊。

平滑濾波的濾波方法

圖像的雜訊濾波器有很多種，常用的有線性濾波器，非線性濾波器。採用線性濾波如鄰域平滑濾波，對受到雜訊污染而退化的圖像復原，在很多情況下是有效的。但大多數線性濾波器具有低通特性，去除雜訊的同時也使圖像的邊緣變模糊了。而另一種非線性濾波器如中值濾波，在一定程度上可以克服線性濾波器所帶來的圖像模糊問題，在濾除雜訊的同時，較好地保留了圖像的邊緣信息。

鄰域平滑濾波原理

鄰域平均法[2]是一種利用Box模版對圖像進行模版操作（卷積運算）的圖像平滑方法，所謂Box模版是指模版中所有係數都取相同值的模版，常用的3×3和5×5模版如下：

鄰域平均法的數學含義是：

（式4-1）

式中：x,y=0,1,…,N-1;S是以（x,y）為中心的鄰域的集合，M是S內的點數。

鄰域平均法的思想是通過一點和鄰域內像素點求平均來去除突變的像素點，從而濾掉一定雜訊，其優點是演算法簡單，計算速度快，其代價會造成圖像在一定程度上的模糊。

中值濾波原理

中值濾波[2]就是用一個奇數點的移動窗口，將窗口的中心點的值用窗口內的各點中值代替。假設窗口內有五點，其值為80、90、200、110和120，那麼此窗口內各點的中值及為110。

設有一個一維序列f1,f2,…,fn,取窗口長度（點數）為m（m為奇數）,對其進行中值濾波，就是從輸入序列中相繼抽出m個數fi-v,…,fi-1,fi,fi+1,…,fi+v(其中fi為窗口中心值,v=(m-1)/2),再將這m個點按其數值大小順序排序，取其序號的中心點的那個數作為濾波輸出。數學公式表示為：

Yi=Med{fi-v,…,fi-1,fi,fi+1,…,fi+v} i∈N v=(m-1)/2 （式4-2）

Yi稱為序列fi-v,…,fi-1,fi,fi+1,…,fi+v的中值

例如，有一序列{0，3，4，0，7}，重新排序後為{0，0，3，4，7}則Med{0，0，3，4，7}=3。此列若用平滑濾波，窗口也取5，那麼平滑濾波輸出為（0+3+4+0+7）/5=2.8。

把一個點的特定長度或形狀的鄰域稱作窗口。在一維情況下，中值濾波器是一個含有奇數個像素的滑動窗口。中值濾波很容易推廣到二維，此時可以利用二維形式的窗口。

對於平面圖像採用的二維中值濾波可以由下式表示：

（式4-3）

式中：A為窗口，{fij}為二維數據序列,即數字圖像各點的灰度值。

對於本系統，由於採集到的是24位真彩色圖像，每個像素點分別有R、G、B三個灰度分量，故要在窗口內分別找到這三個分量的中值，分別用這三個中值去代替窗口中心像素點的R、G、B三個灰度分量的值。

Gabor濾波器

Gabor函數是一個用於邊緣提取的線性濾波器，十分適合紋理表達和分離，其頻率和方向表達同人類視覺系統類似。另外從生物學實驗發現，Gabor濾波器而已很好地近似單細胞的感受野函數。

另外，Gabor函數分為實部和虛部，用實部進行濾波後圖像會平滑；虛部濾波後可用來檢測邊緣（也許這是一個傳說），不過的確可以實踐一下。

Gabor變換時基於傅里葉變換的局限性提出來的，見下圖：

其較為詳細的說明可參見博文

為解決傅里葉變換中存在的局部性問題，其變換的基本思想是：把信號劃分成許多小的時間間隔，用傅里葉變換分析每一個時間間隔中的頻域信息，其處理方法是對傅里葉變換函數設置一個滑動窗口，再對其進行傅里葉變換。其表達如下圖所示：

原傅里葉變換：

比較可見，Gabor變換通過增加一個窗函數來實現了局部化的操作，充分體現了時域與頻域的局部化信息，同時，能夠在整體上提供信號的全部信息。

plc模擬量輸入濾波程序和方案？

1，硬體配置濾波， 如果是200PLC打開系統塊，再Analog里設定濾波時間和頻率 如果是300400PLC打開硬體配置，再相關模塊里設定濾波時間和頻率，這個一般是過濾高頻的雜波 2，然後再程序里，編程實現： 均值濾波：我一般用最後五次採樣的平均值，採樣時間間隔和幾次求平均值可以自己定。 中值濾波：我沒用過，可以嘗試。 峰值濾波：直接取多次採樣的最高或最低值，也是特殊情況有用的。 總結：你首先要觀察你的測量量的特性，否則濾波是低效、盲目的。