本發(fā)明屬于數字全息技術領域，尤其涉及一種基于頻譜能量的數字全息圖諧波探測與消除方法。

背景技術：

利用數字全息技術進行顯微成像或精密測量時，數字全息圖的質量對于還原像的質量或測量精度影響很大，激光散斑、參考光傾斜、諧波的存在會大幅降低成像質量或測量精度。數字全息需要物光與參考光干涉記錄干涉圖，相對于普通照相，記錄光路中需要更多的光學元件，如分束器、合束器、空間濾波器等。光波經過這些器件時，其前后表面會發(fā)生反射、折射等現象。在棱鏡或其他鏡面之間多次反射的光波會發(fā)生干涉，給全息圖附加周期性的正余弦分布的直條紋，被稱為諧波。這種諧波是由于分束器兩平面不嚴格平行，或某個或幾個光學器件表面不平行造成的多次反射的具有一定夾角的平面波干涉引起的。這一干擾信息會嚴重干擾全息圖的準確定位并造成全息再現時最終物光波的畸變或測量失敗。要消除這些諧波的影響，需要對全息圖進行預處理，消掉全息圖中的這些諧波條紋。現有全息圖預處理技術使用傅里葉變換在頻譜面上進行濾波消掉諧波的頻譜后再逆傅里葉變換獲得無諧波的全息圖。這一諧波消除方法需要設計相應的頻譜濾波器進行濾波，去掉諧波頻譜分量。這一濾波方法無法準確消除諧波成分，會剩余諧波分量，或者誤傷有用的成分。而且，濾波器的濾波范圍和濾波幅度也會影響濾波效果：濾波幅度較小時，不能完全除掉諧波信息，濾波幅度較大時，會增加錯誤信息。濾波范圍小時，無法消除諧波影響；范圍較大時，也會濾掉有用的物光頻譜信息。濾波法的第二個技術缺陷是計算工作量大，濾波后還要進行一次逆傅里葉變換，獲得沒有諧波的全息圖，兩次傅里葉變換比一次傅里葉變換增大一倍的計算量。這也導致了濾波法的第三個缺陷：無法應用到連續(xù)的視頻圖像的諧波消除中。由于傅里葉變換計算量大，每幅圖像都進行兩次傅里葉變換需要大量的計算時間，即使是一次傅里葉變換，通常情況下也無法達到視頻播放的速度。

綜上所述，現有技術存在的問題是：

現有全息圖預處理技術存在設計濾波器過程繁瑣，成本高，尋找頻譜困難；難以確定濾波的范圍；難以確定濾波幅度；容易丟失有用的物光信息；需要對濾波后的頻譜圖進行逆傅里葉變換獲得處理后的全息圖，增加一倍的計算量；無法應用到在線連續(xù)圖像處理中。

技術實現要素：

針對現有技術存在的問題，本發(fā)明提供了一種基于頻譜能量的數字全息圖諧波探測與消除方法。

本發(fā)明是這樣實現的，一種基于頻譜能量的數字全息圖諧波探測與消除方法，所述基于頻譜能量的數字全息圖諧波探測與消除方法包括圖像的傅里葉變換、諧波頻譜的尋找、讀取諧波頻譜的坐標和復數值、計算諧波在水平豎直兩個方向上的頻率和干擾強度系數及初相位、利用參數構建諧波分布、圖像校正六個步驟完成有諧波干擾的圖像的處理。

進一步，所述基于頻譜能量的數字全息圖諧波探測與消除方法具體包括以下步驟:

第一步，圖像的傅里葉變換：使用計算機內的快速二維傅里葉變換算法完成對圖像的傅里葉變換；諧波在原圖像上是有一定偏置的正余弦周期性分布，其傅里葉頻譜的模(強度)是三個對稱分布的δ函數，其中一個在頻譜面上的原點，另外兩個關于原點對稱分布；

第二步，諧波頻譜的尋找：由于δ函數是無限小的點，在頻譜面上的諧波的頻譜是三個很亮的點，零頻的點在原點，另兩個點對稱分布在原點的兩側；先分析出只包含該亮點的大體區(qū)域，再用依次求最大值的算法尋找；

第三步，讀取諧波頻譜的坐標和復數值：

用尋找最大值的方法得到諧波的非零頻譜的水平和豎直方向的頻譜坐標，分別設為u1和v1；再讀出該像素頻譜函數值的實部和虛部，分別設為er和ei，設其模為e1；另一非零頻頻譜坐標對稱分布，頻譜函數值相等，不用讀?。煌瑯?，原點的頻譜坐標為(0，0)，無需讀取，但需要讀取原點的頻譜的函數值，這一函數值對應全息圖的零級頻譜，設為e0，其函數值為正實數；

第四步，計算諧波在水平豎直兩個方向上的頻率和干擾強度系數及初相位：

假設諧波頻譜在頻譜面上水平方向上的坐標為u1，在豎直方向上的頻譜坐標為v1；求出水平方向的諧波頻率和在豎直方向上的諧波頻率fx和fy；

假設諧波和全息圖零級項對應的頻譜函數值的模分別為e1和e0，則由這兩個參數求得諧波的干擾強度系數m；

假設諧波頻譜對應函數值的實部和虛部分別設為er和ei；求出諧波對應的初相位

第五步，利用參數構建諧波分布：

利用得到的水平方向的諧波頻率和在豎直方向上的諧波頻率和干擾強度系數及初相位四個參數組建諧波的強度分布；

第六步，圖像校正：

求出諧波的強度分布后，通過圖像相減完全去掉強度圖中的諧波成分，得到消除諧波影響后的干涉圖。

進一步，所述圖像的傅里葉變換具體包括：

設圖像的強度分布為i(x,y)，使用計算機內的快速二維傅里葉變換算法完成對圖像的傅里葉變換，即：

e(u,v)＝fft2[i(x,y)](1)

在這些頻譜中，含有物光信息的頻譜是連續(xù)的，且分布范圍較大；諧波在圖像上是有一定偏置的正余弦函數分布；諧波的傅里葉頻譜的模是三個對稱分布的δ函數，其中一個在頻譜面上的原點，另外兩個關于原點對稱分布。

進一步，用依次求最大值的算法尋找頻譜坐標的方法包括：

對全息圖，零頻分量的能量遠大于高頻分量；對數字全息圖進行傅里葉變換后，零級分量最強，利用數值算法對整個頻譜面復數組取模，再求最大值，只能得到零級頻譜；

先將零級頻譜賦值為0，再找最大值，找到兩個諧波譜中的一個，記下其坐標和復數值；最后把零頻的值再恢復到原值。

現有技術是通過人眼觀察頻譜位置，依據觀察的位置設計濾波器，比較繁瑣。

進一步，所述一定偏置來自于光的直透項或常數項。

進一步，所述計算諧波在水平豎直兩個方向上的頻率和干擾強度系數，以及初相位，具體包括：

假設諧波頻譜在頻譜面上水平方向上的坐標為u1，在豎直方向上的頻譜坐標為v1，求出水平方向的諧波頻率為：

fx＝u1(2)

在豎直方向上的諧波頻率為：

fy＝v1(3)

假設諧波和全息圖的零級項對應的頻譜函數值的模分別為e1和e0，則由這兩個參數求得諧波的干擾強度系數為：

m＝2e1/(e0-2e1)(4)。

假設諧波對應頻譜的頻譜函數值的實部與虛部分別為er和ei，則諧波的初相位是：

式中arg()表示對括號內的復數取幅角。

進一步，利用所述參數構建諧波分布，具體包括：

利用公式(2)、(3)、(4)、(5)得到的四個參數組建諧波的強度分布，方程如下：

式中<>符號表示對所有像素取平均值，坐標x和y以像素為坐標單位。

進一步，所述圖像校正，具體包括：

由公式(6)求出諧波的強度分布后，通過圖像相減完全去掉強度圖i中的諧波成分：

i1就是消除諧波影響后的干涉圖；

存在多個諧波或連續(xù)譜時，依次找到其頻譜坐標或連續(xù)的頻譜坐標范圍，記錄坐標和相應的頻譜函數值，求出每一項諧波的干擾強度系數和初相位，構建每一項強度分布，最后依次從帶有諧波的圖像中減掉；

對于非諧波的周期性干擾，展開成周期性諧波疊加的形式進行求解消除非諧波的影響；對于非周期性非簡諧波干擾，通過傅里葉級數展開成不同諧波的疊加來求解。

本發(fā)明的優(yōu)點及積極效果為：

不需要設計濾波器，只利用零級傅里葉變換譜和諧波譜的能量關系和諧波譜的復數值即可。

本發(fā)明無需獲得諧波對應頻譜的范圍。本發(fā)明能夠精確消除諧波成分，不會對原物光的信息有任何影響。本發(fā)明不需要對濾波后的頻譜圖進行逆傅里葉變換，可以使計算速度提高近50％。本發(fā)明利用能量關系獲得諧波參數后，構建諧波強度分布，直接相減消掉諧波，計算簡單。如果是連續(xù)記錄的圖像，只對一幅圖像進行一次傅里葉變換，求出諧波參數后，可以對多幅帶有諧波的圖像進行快速消除，需要的運算只是一次圖像相減。

本發(fā)明計算出諧波參數后，反演出諧波強度分布，直接進行圖像相減，完成諧波的消除。只需對一幅圖像進行一次傅里葉變換即可，不需要對所有圖像進行傅里葉變換和逆傅里葉變換，大大節(jié)省了計算時間。特別是在視頻顯示方面，濾波法需要對大量圖片進行傅里葉變換和逆傅里葉變換，計算負擔大，而本發(fā)明中提出的方法，只進行一次傅里葉變換后，直接進行圖片相減即可，極大地減輕了視頻處理的計算負擔。本發(fā)明的方法不僅適用于全息圖諧波的探測與消除，也適用于其它數字圖像和視頻的諧波探測與處理，處理方法相同。

附圖說明

圖1是本發(fā)明實施例提供的基于頻譜能量的數字全息圖諧波探測與消除方法流程圖。

圖2是本發(fā)明實施例提供的原始物光強度和相位分布圖；

圖中：(a)強度分布；(b)相位分布。

圖3是本發(fā)明實施例提供的棱鏡前后表面平行與不平行情況對比圖；

圖中：(a)平行時全息圖；(b)不平行時全息圖。

圖4是本發(fā)明實施例提供的傅里葉頻譜圖。

圖5是本發(fā)明實施例提供的校正后的全息圖。

圖6是本發(fā)明實施例提供的未校正還原像強度與相位分布圖；

圖中：(a)強度圖；(b)相位圖。

圖7是本發(fā)明實施例提供的校正后還原像強度與相位分布圖；

圖中：(a)強度圖；(b)相位圖。

具體實施方式

為了使本發(fā)明的目的、技術方案及優(yōu)點更加清楚明白，以下結合實施例，對本發(fā)明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

本發(fā)明實施例提供的基于頻譜能量的數字全息圖諧波探測與消除方法包括圖像的傅里葉變換、諧波頻譜的尋找、讀取諧波頻譜的坐標和復數值、計算諧波在水平豎直兩個方向上的頻率和干擾強度系數及初相位參數、利用參數構建諧波分布、圖像校正六個步驟完成有諧波干擾的圖像的處理。

下面結合附圖對本發(fā)明的應用原理作詳細描述。

如圖1所示，本發(fā)明實施例提供的基于頻譜能量的數字全息圖諧波探測與消除方法包括以下步驟：

s101:圖像的傅里葉變換：使用計算機內的快速二維傅里葉變換算法完成對圖像的傅里葉變換；諧波在圖像上是有一定偏置(來自于光的直透項或常數項)的正余弦周期性分布，其傅里葉頻譜的模是三個對稱分布的δ函數，其中一個在頻譜面上的原點，另外兩個關于原點對稱分布。

s102:諧波頻譜的尋找：由于δ函數是無限小的點，在頻譜面上的諧波的頻譜的模是三個很亮的點，零頻的點在原點，另兩個點對稱分布在原點的兩側；先觀察出包含該亮點的大體區(qū)域，再用求最大值的算法尋找。

s103:讀取諧波頻譜的坐標和復數值：用依次尋找最大值的方法得到諧波的非零頻譜的水平和豎直方向的頻譜坐標，分別設為u1和v1；再讀出該像素頻譜函數值，其實部和虛部分別設為er和ei，則該頻譜函數的模為復數值的模e1＝|er+jei|，也就是該頻譜的能量或強度；另一非零頻頻譜坐標對稱分布，頻譜函數值相等，不用讀取；同樣，原點的頻譜坐標為(0，0)，無需讀取，但需要讀取原點的頻譜的函數值，這一函數值對應全息圖的零級頻譜，為實數，設為e0。

s104:計算諧波在水平豎直兩個方向上的頻率和干擾強度系數及初相位：假設諧波頻譜在頻譜面上水平方向上的坐標為u1，在豎直方向上的頻譜坐標為v1；求出水平方向的諧波頻率和在豎直方向上的諧波頻率；假設諧波和全息圖零級頻譜對應的頻譜函數值的模分別為e1和e0，則由這兩個參數求得諧波的干擾強度系數m；假設諧波頻譜對應函數值的實部和虛部分別設為er和ei；求出諧波對應的初相位

s105:利用參數構建諧波分布：利用得到的參數求出水平方向的諧波頻率和在豎直方向上的諧波頻率和干擾強度系數及初相位四個參數組建諧波的強度分布。

s106:圖像校正：求出諧波的強度分布后，通過圖像相減完全去掉強度圖中的諧波成分，得到消除諧波影響后的干涉圖。

用依次求最大值的算法尋找頻譜坐標的方法包括：

對全息圖，零頻分量的能量遠大于高頻分量；對數字全息圖進行傅里葉變換后，零級分量最強，利用數值算法對整個頻譜面復數組取模再求最大值，只能得到零級頻譜；

先將零級頻譜賦值為0，再找最大值，找到兩個諧波譜中的一個，記下其坐標和復數值；最后把零頻的值再恢復到原值。

現有技術是通過人眼觀察頻譜位置，依據觀察的位置設計濾波器，比較繁瑣。

下面結合具體實施例對本發(fā)明的應用原理作進一步描述。

本發(fā)明實施例提供的基于頻譜能量的數字全息圖諧波探測與消除方法，包括：

第一步，圖像的傅里葉變換：

設圖像的強度分布為i(x,y)，使用計算機內的快速二維傅里葉變換算法完成對圖像的傅里葉變換，即：

e(u,v)＝fft2[i(x,y)](1)

在這些頻譜中，含有物光信息的頻譜是連續(xù)的，且分布范圍較大。諧波在圖像上是有一定偏置量(來自于光的直透項或常數項)的正余弦周期性分布，其傅里葉頻譜的模是三個對稱分布的δ函數，其中一個在頻譜面上的原點，另外兩個關于原點對稱分布。

第二步，諧波頻譜的尋找。

由于δ函數是無限小的點，在頻譜面上的諧波的頻譜的模是三個很亮的點，零頻的點在原點，另兩個點對稱分布在原點的兩側。先觀察出只包含該亮點的大體區(qū)域，再用求最大值的算法尋找。

第三步，讀取諧波頻譜的坐標和復數值。

用尋找最大值的方法得到諧波的非零頻譜的水平和豎直方向的頻譜坐標，分別設為u1和v1。再讀出該頻譜函數值，其實部和虛部分別設為er和ei，則該頻譜函數的模為復數值的模e1＝|er+jei|，也就是該頻譜的能量或強度。另一非零頻譜的坐標對稱分布，頻譜函數值相等，不用讀取。同樣，原點的頻譜坐標為(0，0)，無需讀取，但需要讀取原點的頻譜的函數值，這一函數值對應全息圖的零級頻譜，為實數，設為e0。

第四步，計算諧波在水平豎直兩個方向上的頻率和干擾強度系數及初相位參數。

假設諧波頻譜在頻譜面上水平方向上的坐標為u1，在豎直方向上的頻譜坐標為v1，可求水平方向的諧波頻率為：

fx＝u1(2)

在豎直方向上的諧波頻率：

fy＝v1(3)

假設在諧波和全息圖的零級項對應的頻譜函數的模分別為e1和e0，則可以由這兩個參數求得諧波的干擾強度系數為：

m＝2e1/(e0-2e1)(4)

假設諧波頻譜復數值對應的實部和虛部分別為er和ei，則可以求出與諧波對應的初相位為：

式中“arg()”表示對復數取幅角。

第五步，利用參數構建諧波分布。

利用公式(2)、(3)、(4)、(5)得到的四個參數組建諧波的強度分布方程如下：

式中<>符號表示對所有像素取平均值，坐標x和y以像素為坐標單位。

第六步，圖像校正。

由公式(6)求出諧波的強度分布后，就可以通過圖像相減完全去掉強度圖i中的諧波成分：

i1就是消除諧波影響后的干涉圖。

下面結合模擬實驗對本發(fā)明的應用效果作詳細的描述。

在廣義相移數字全息實驗中，參考光一般是濾波準直后的平面光波，如果棱鏡前后兩表面不嚴格平行，透射過分束器棱鏡時會形成等厚干涉斜條紋。隨著參考光強度的增加，干涉條紋會越強，意味著在全息圖的記錄中干擾信息會更強，在全息圖還原時干擾信息對還原結果的影響也會越大，尤其是單步還原和兩步還原算法，會導致還原結果嚴重扭曲，甚至無法還原出結果。

本發(fā)明不需要更多的操作也不需要采集更多的數據，用兩步還原算法所需的兩幅全息圖中的一幅就可以得到需要的信息，通過計算機運算就可以有效降低棱鏡前后兩表面不平行造成的誤差。接下來的內容中，本發(fā)明將對這種新算法進行計算機模擬并比較有校正與無校正的兩種情況下還原像的質量。

在計算機模擬實驗中設置物平面上有一球面波。模擬的物平面和像平面的大小都是512×512像素，像素尺寸為15×15微米。如圖2所示，(a)為物光強度圖，振幅分布滿足高斯分布，(b)為物光相位。原始物光經過菲涅爾衍射后與參考光在記錄面上形成干涉全息圖，見圖3(a)。由于棱鏡前后表面不平行，導致在全息圖上形成與全息圖無關的斜條紋，如圖3(b)所示：

圖3(a)中棱鏡前后表面平行，可以看到在全息圖上沒有諧波。圖3(b)中棱鏡前后表面存在很小的傾斜角度，造成了全息圖上干擾性的斜條紋。在模擬實驗中，本發(fā)明設置的相移值為1.5rad、干擾強度系數m值設置為0.2、諧波初相位設置為3rad，對全息圖圖3(b)進行傅里葉變換，其頻譜函數值的強度分布如圖4所示。

圖4中，頻譜圖上有兩個極大值點。由于零級譜太強，在顯示圖4時，零級譜被我們去掉了。利用這兩個點的坐標值和諧波譜的模及幅角便可以求出諧波的水平方向和豎直方向的空間頻率，以及校正公式中需要的m值，初相位

m的計算值為0.2032，初相位的計算值為3.0018rad，這與本發(fā)明的設定值誤差很小，將fx、fy、m、的值代入校正算法中，可以將有干擾信息的全息圖中干擾信息有效的消除，圖5給出了一幅校正后的全息圖，完全消除了諧波的影響。下面本發(fā)明將對無校正時的還原圖像(圖6)與校正后的還原結果(圖7)進行對比。

在圖6與圖7之間的對比中，可以看出未校正的還原物光波無論是強度分布還是相位分布都受到了很大影響，尤其是相位分布圖發(fā)生了嚴重的扭曲，而在校正后的物光波還原圖中，強度分布與相位分布都得到了很好的校正。

本發(fā)明不需要設計濾波器，只利用全息圖零級傅里葉變換譜和諧波譜的能量關系即可。

本發(fā)明無需獲得諧波對應頻譜的范圍。本發(fā)明不會對原物光的信息有任何影響。本發(fā)明不需要對濾波后的頻譜圖進行逆傅里葉變換。本發(fā)明利用能量關系獲得諧波參數后，構建諧波強度分布，直接相減消掉諧波，計算簡單。

本發(fā)明計算出諧波參數后，反演出諧波強度分布，直接進行圖像相減，完成諧波的消除。只需對一幅圖像進行一次傅里葉變換即可，不需要對所有圖像進行傅里葉變換和逆傅里葉變換，大大節(jié)省了計算時間，至少節(jié)省50％以上的計算資源。特別是在視頻顯示方面，濾波法需要對大量圖片進行傅里葉變換和逆傅里葉變換，計算負擔大，而本發(fā)明中提出的方法，只進行一次傅里葉變換后，直接進行圖片相減即可，極大地減輕了視頻處理的計算負擔。本發(fā)明的方法不僅適用于全息圖像諧波的探測與消除，也適用于其它數字圖像和視頻諧波的探測與處理，處理方法相同。

存在多個諧波或連續(xù)譜時，依次找到其頻譜坐標或連續(xù)的頻譜坐標范圍，記錄坐標和相應的頻譜函數值，求出每一項諧波的干擾強度系數和初相位，構建每一項強度分布，最后依次從帶有諧波的圖像中減掉即可。

對于非諧波的周期性干擾，可以展開成周期性諧波疊加的形式進行求解消除其影響。對于非周期性非簡諧波干擾，可以通過傅里葉級數展開成不同諧波的疊加來求解。

以上所述僅為本發(fā)明的較佳實施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內所作的任何修改、等同替換和改進等，均應包含在本發(fā)明的保護范圍之內。