阅读说明：本技术 一种傅里叶及窗口傅里叶变换的联合相位重构算法 (Fourier and window Fourier transform combined phase reconstruction algorithm ) 是由 刘丙才 陈瑜 潘永强 朱学亮 田爱玲 王红军 岳鑫 于 2020-04-23 设计创作，主要内容包括：本发明公开了一种傅里叶及窗口傅里叶变换的联合相位重构算法,其特征在于,针对单帧载频干涉条纹图,利用傅里叶空频分析法滤除背景光强及其他无关项的干扰,利用傅里叶逆变换处理将其转换到空域内；随后对其进行二维窗口傅里叶变换处理,将得到的窗口傅里叶频谱分布中的最大值进行提取,滑动窗口函数获得整幅干涉条纹的窗口傅里叶脊分布,包裹相位可对窗口傅里叶脊求相位角计算求取。随后对该包裹相位进行二维解包裹处理,获得所需连续相位。本发明与传统WFT相比,克服了背景光强干扰问题,并通过正一级频谱移中处理化简了载频条纹在面形拟合时存在的消倾斜等繁杂过程,显著提高了WFT相位重构技术的精度,可进一步应用于大口径光学元件测量中。(The invention discloses a combined phase reconstruction algorithm of Fourier and window Fourier transform, which is characterized in that aiming at a single-frame carrier frequency interference fringe pattern, a Fourier space-frequency analysis method is utilized to filter the interference of background light intensity and other irrelevant items, and inverse Fourier transform processing is utilized to convert the interference into a space domain; and then, carrying out two-dimensional window Fourier transform processing on the interference fringes, extracting the maximum value in the obtained window Fourier spectrum distribution, obtaining the window Fourier ridge distribution of the whole interference fringes by a sliding window function, and calculating and solving a phase angle of the window Fourier ridge by a wrapping phase. The wrapped phase is then two-dimensional unwrapped to obtain the desired continuous phase. Compared with the traditional WFT, the method overcomes the problem of background light intensity interference, simplifies the complex processes of inclination elimination and the like existing in the process of surface shape fitting of carrier frequency stripes by processing in the primary frequency spectrum shift, obviously improves the precision of the WFT phase reconstruction technology, and can be further applied to the measurement of large-caliber optical elements.)

1. A Fourier and window Fourier transform joint phase reconstruction algorithm is characterized in that: the algorithm comprises the following steps:

the method comprises the following steps: carrying out two-dimensional Fourier transform on the collected or simulated initial carrier frequency fringes to obtain the frequency spectrum distribution of the initial carrier frequency fringes;

the carrier frequency fringe intensity is expressed as:

wherein the content of the first and second substances,I_d(x, y) is a background light intensity distribution of the interference fringes, b (x, y) is a modulation degree distribution of the interference fringes,for a phase distribution function containing phase information of the wave surface to be measured, f_x、f_yRepresenting the complex conjugate for the spatial carrier frequency in the x, y directions;

step two: and (2) aiming at the single-frame carrier frequency stripe, carrying out two-dimensional Fourier transform processing on the single-frame carrier frequency stripe, wherein the frequency domain distribution corresponding to the carrier frequency stripe is represented as follows:

F(f₁,f₂)＝A(f₁,f₂)+C(f₁-f_x,f₂-f_y)+C*(f₁+f_x,f₂+f_y) (2)

wherein, the equal-sign right side is from left to right: the direct current component, a positive primary frequency spectrum and a negative primary frequency spectrum, and the size of the carrier frequency controls the distance between the positive primary frequency spectrum and the negative primary frequency spectrum;

step three: selecting a proper type of a filtering window according to the characteristics of the primary frequency spectrum distribution of different fringe patterns; then, the central bandwidth of the corresponding filter window is selected, and the value of the central bandwidth is equal to the central position (f) of the primary spectrum_x,f_y) Keeping consistency, repeatedly debugging and extracting complete primary frequency spectrum and filtering background light intensity and irrelevant items, moving primary frequency spectrum component to central position to obtain interference fringe component with background light intensity item filtered, and marking as C (f)₁,f₂)；

Step four: and (3) performing inverse Fourier transform processing on the positive-level frequency spectrum extracted in the step three, transferring the positive-level frequency spectrum from a frequency domain into a space domain, and obtaining a distribution component containing phase information, which is marked as c (x, y) and is expressed as:

step five: performing window Fourier transform phase extraction processing on the spatial domain component of the positive primary spectrum containing the phase information obtained in the step four, wherein for the WFT phase extraction technology, the input end is the c (x, y), and the corresponding window Fourier transform is expressed as

According to multiple times of debugging of simulation and experiment, selecting a default window type used by a WFT algorithm as a Gaussian window, wherein one-dimensional expression of the default window type is

Where σ represents the extension of the gaussian window function in the x-direction, this can be understood as the size of the gaussian window.

Rotating or point-multiplying the one-dimensional window function of the formula (5) to obtain the corresponding two-dimensional distribution, wherein the corresponding two-dimensional Gaussian window function can be expressed as

g(x,y)＝g_x(x)g_y(y) (6)

Wherein, g_x(x)＝g_y(y) g (x), define σ_x、σ_yFor this purpose, the two-dimensional window function extends in the x, y direction. In the precision analysis of the window size in the previous period, the optimal Gaussian window size when the stripe linear phase error suppression degree can be maximized is determined to be sigma_x＝σ_y＝5pixel；

C (x, y) and a two-dimensional Gaussian window function g (x, y) are brought into a formula (4) and are arranged to obtain

Sf₀(u,v；ξ,η)＝c(x,y)G_x,y(u,v；ξ,η) (7)

G_x,y(u,v；ξ,η)＝G_x,y(u,v；ξ)G_x,y(u,v；η) (8)

Wherein

Definition G_x,y(u,v；ξ)、G_x,y(u, v; η) are one-dimensional complex gain factors along the x, y axes, respectively, and G_x,y(u, v; ξ) is a two-dimensional complex gain factor.

In formulas (9) and (10)Respectively, local frequencies along the x, y axes at the pixel points (u, v) which are numerically equal to the first partial derivative of the phase;respectively, the local curvature along the x, y axis at the pixel point (u, v), which is numerically equal to the phaseThe second partial derivative of (d);

when satisfying (ξ) ═ ω_x,ω_y) Time, two-dimensional gain factor G_x,y(u, v; ξ) takes a maximum value of

Determining an argument value at the maximum spectrum value (or ridge value) in the coverage area of the window, assigning the argument value as a local frequency (or phase derivative) in the window, accumulating and superposing ridge values in all the windows through a sliding window to obtain a window Fourier ridge value distribution of the whole carrier frequency stripe, and performing phase angle calculation on the window Fourier ridge value distribution to obtain a wrapping phase:

step six: and unpacking the wrapped phase by using a two-dimensional discrete cosine phase unwrapping algorithm to obtain a continuous phase, and reconstructing three-dimensional surface distribution by using a Zernike polynomial.

Technical Field

The invention relates to the field of optical interferometry, in particular to a Fourier and window Fourier transform combined phase reconstruction algorithm.

Background

Optical interferometry as a non-contact measurement method is defined as an important development direction of precision measurement technology in the new century. Phase extraction is a key step in optical interferometry. The traditional phase extraction technology mainly comprises a phase shift method and a Fourier transform method, but the phase shift method and the Fourier transform method have limitations, the phase shift method is easily disturbed by atmosphere when measuring a large-caliber optical element, and the data processing is complex; the Fourier transform method effectively overcomes the mechanical vibration and air flow influence existing in the phase shift method, but the Fourier transform is a global analysis technology, so that the noise in a local area is not thoroughly filtered, and the phase reconstruction precision is reduced. The window Fourier transform phase extraction algorithm proposed in recent years solves the problem of local analysis to a certain extent, enables interference fringes not to be influenced by noise through local (window) analysis, and becomes one of the current common methods for measuring large-aperture optical elements.

When a certain interference fringe is processed by using a window Fourier phase extraction technology, fringe components contained in a default window are in linear distribution, but actual carrier frequency fringes are in nonlinear distribution, so that an algorithm generates a linear phase error during operation. At present, for the linear phase error processing of interference fringes, the size of a window is generally reduced, so that the interference fringes in the window are approximately considered to be linear, but the size of a small-size window in a frequency domain is increased rapidly, so that the low-frequency components of the fringes in a region are interfered by background light intensity, and the algorithm has long running time under the processing of the small-size window, which is not beneficial to practical application.

Therefore, it is an urgent technical problem to be solved by those skilled in the art to provide a new fourier and windowed fourier transform joint phase reconstruction algorithm.

Disclosure of Invention

In order to solve the technical problems, the invention provides a combined phase reconstruction algorithm of Fourier transform and window Fourier transform

The invention provides a Fourier and window Fourier transform combined phase reconstruction algorithm, which comprises the following steps:

step 1: carrying out two-dimensional Fourier transform on the collected or simulated initial carrier frequency fringes to obtain the frequency spectrum distribution of the initial carrier frequency fringes;

the carrier frequency fringe intensity is expressed as:

wherein the content of the first and second substances,I_d(x, y) is a background light intensity distribution of the interference fringes, b (x, y) is a modulation degree distribution of the interference fringes,for a phase distribution function containing phase information of the wave surface to be measured, f_x、f_yRepresenting the complex conjugate for the spatial carrier frequency in the x, y directions;

step 2: the single-frame carrier frequency stripe is subjected to two-dimensional Fourier transform processing, due to the introduction of carrier frequency components, the positive and negative primary frequency spectrums in the original exponential phase field are subjected to position offset, and a state of separating a positive one from a negative one is presented in a corresponding frequency domain, and in this state, the corresponding frequency domain distribution of the carrier frequency stripe can be expressed as:

F(f₁,f₂)＝A(f₁,f₂)+C(f₁-f_x,f₂-f_y)+C*(f₁+f_x,f₂+f_y) (2)

wherein, the equal-sign right side is from left to right: the direct current component, the positive primary frequency spectrum and the negative primary frequency spectrum, and the size of the carrier frequency controls the distance between the positive primary frequency spectrum and the negative primary frequency spectrum.

And step 3:selecting a proper type of a filtering window according to the distribution characteristics of the primary spectrum of different fringe patterns, wherein the commonly used filtering windows such as a Hanning window, a Hamming window and a Blackman filtering window can be selected according to the distribution characteristics of the corresponding primary spectrum; then, the central bandwidth of the corresponding filter window is selected, and the value of the central bandwidth is equal to the central position (f) of the primary spectrum_x,f_y) Keeping consistency, repeatedly debugging and extracting complete primary frequency spectrum and filtering background light intensity and irrelevant items, moving primary frequency spectrum component to central position to obtain interference fringe component with background light intensity item filtered, and marking as C (f)₁,f₂)。

And 4, step 4: and (3) performing inverse Fourier transform processing on the positive-level frequency spectrum extracted in the step (3), transferring the positive-level frequency spectrum from the frequency domain into a space domain, and obtaining a distribution component containing phase information, which is marked as c (x, y) and is expressed as c (x, y)

And 5: performing window Fourier transform phase extraction processing on the spatial domain component of the positive-order frequency spectrum containing the phase information obtained in the step 4, wherein for the WFT phase extraction technology, the input end is c (x, y) above, and the corresponding window Fourier transform can be expressed as

According to multiple times of debugging of simulation and experiment, selecting a default window type used by a WFT algorithm as a Gaussian window, wherein one-dimensional expression of the default window type is

Where σ represents the extension of the gaussian window function in the x-direction, this can be understood as the size of the gaussian window.

Rotating or point-multiplying the one-dimensional window function of the formula (5) to obtain the corresponding two-dimensional distribution, wherein the corresponding two-dimensional Gaussian window function can be expressed as

g(x,y)＝g_x(x)g_y(y) (6)

Wherein, g_x(x)＝g_y(y) g (x), define σ_x、σ_yFor this purpose, the two-dimensional window function extends in the x, y direction. In the precision analysis of the window size in the previous period, the optimal Gaussian window size when the stripe linear phase error suppression degree can be maximized is determined to be sigma_x＝σ_y＝5pixel；

Sf is obtained by carrying c (x, y) and a two-dimensional Gaussian window function g (x, y) into a formula (4) and arranging₀(u,v；ξ,η)＝c(x,y)G_x,y(u,v；ξ,η) (7)

G_x,y(u,v；ξ,η)＝G_x,y(u,v；ξ)G_x,y(u,v；η) (8)

Wherein

Definition G_x,y(u,v；ξ)、G_x,y(u, v; η) are one-dimensional complex gain factors along the x, y axes, respectively, and G_x,y(u, v; ξ) is a two-dimensional complex gain factor.

In formulas (9) and (10)Respectively, local frequencies along the x, y axes at the pixel points (u, v) which are numerically equal to the first partial derivative of the phase;respectively, local curvatures along x and y axes at the pixel points (u, v), which are numerically equal to the second partial derivative of the phase;

when satisfying (ξ) ═ ω_x,ω_y) Time, two-dimensional gain factor G_x,y(u, v; ξ) takes a maximum value of

And determining an argument value at the maximum spectrum value (or ridge value) in the coverage area of the window, assigning the argument value as a local frequency (or phase derivative) in the window, accumulating and superposing ridge values in all the windows through a sliding window to obtain the window Fourier ridge value distribution of the whole carrier frequency stripe, and calculating the phase angle of the window Fourier ridge value distribution to obtain the wrapping phase.

Step 6: and unpacking the wrapped phase by using a two-dimensional discrete cosine phase unwrapping algorithm to obtain a continuous phase, and reconstructing three-dimensional surface distribution by using a Zernike polynomial.

Compared with the related art, the invention has the following beneficial effects:

1. the wrapping phases extracted in the traditional WFT technology all contain carrier frequency components, and artificial inclination elimination processing is needed during subsequent surface shape fitting.

2. The invention solves the problem of background light intensity interference of the WFT phase extraction technology under the condition of selecting a small-size window, and uniformly inhibits the linear phase error of carrier frequency interference fringes under the combined action of the optimal window size, so that the phase extraction precision under the FT-WFT combined algorithm processing is superior to that of the traditional WFT technology.

Drawings

FIG. 1: based on the FT-WFT combined phase reconstruction algorithm flow chart provided by the invention;

FIG. 2: generating a single-frame carrier frequency interference fringe pattern through computer simulation;

FIG. 3: a Fourier transformed carrier frequency interference fringe spectrum distribution diagram;

FIG. 4: positive first-order spectral distribution diagram extracted by rectangular filtering window

FIG. 5: a shifted positive level spectral profile;

FIG. 6: the positive frequency spectrum is converted into an intensity distribution diagram in a space domain through inverse Fourier transform;

FIG. 7: the window Fourier ridge value distribution diagram obtained after the interference fringe component obtained by filtering the background light intensity is processed by WFT

FIG. 8: calculating a phase angle of the window Fourier ridge value to obtain a wrapping phase distribution diagram;

FIG. 9: the generated wrapped phase is subjected to two-dimensional discrete cosine phase unwrapping processing to obtain a continuous phase distribution map;

FIG. 10: processing the wrapped phase distribution map obtained in the figure 2 by a traditional WFT phase reconstruction algorithm;

FIG. 11: the wrapping phase is subjected to Zernike polynomial fitting reconstruction to obtain a three-dimensional surface profile distribution diagram;

Detailed Description

The invention will be further explained with reference to the drawings and the embodiments.

A fourier and windowed fourier transform joint phase reconstruction algorithm, as shown in fig. 1, includes the following detailed steps:

step 1: carrying out two-dimensional Fourier transform on the single-frame carrier frequency interference fringe pattern to obtain frequency spectrum distribution; simulating a carrier frequency interference fringe pattern with the size of 300 × 300 pixels, and adding white gaussian noise with the mean value μ being 0 and the variance σ being 1 to obtain an initial interference fringe pattern, as shown in fig. 2; performing two-dimensional Fourier transform on the spectrum to obtain two-dimensional spectrum distribution of the spectrum, and obtaining a corresponding three-dimensional spectrum three-dimensional distribution top view as shown in FIG. 3;

step 2: extracting and transferring a primary spectrum; selecting a filtering window as a rectangular window according to the primary spectrum distribution characteristics, determining the size of the rectangular window to be 24 pixels through computer debugging, performing dot multiplication operation on an interference fringe Fourier spectrum and the filtering window in a frequency domain to obtain a fundamental frequency spectrum for filtering background light intensity and irrelevant items, as shown in FIG. 4, and moving the fundamental frequency spectrum to the central position of an area, as shown in FIG. 5;

and step 3: performing two-dimensional Fourier inverse transformation to obtain spatial distribution; performing two-dimensional inverse Fourier transform processing on the shifted positive-order frequency spectrum to obtain an intensity distribution diagram in a space domain, as shown in FIG. 6;

and 4, step 4: performing two-dimensional Fourier transform on the two-dimensional fringe component obtained by IFT to obtain a window Fourier spectrum; performing two-dimensional window Fourier transform on the carrier frequency fringe signal after the frequency spectrum processing, wherein the window type is selected as a Gaussian window, the size of the window is 5 pixels, so as to obtain a window Fourier spectrum of the signal, searching for an independent variable value which can make the spectrum value in a window area reach the maximum through a maximum spectrum value, and recording the independent variable value as the local frequency of the point, wherein the ridge value distribution is shown in FIG. 7; calculating the phase angle of the ridge value to obtain the wrapping phase, as shown in fig. 8;

and 5: extracting a ridge value (maximum spectral value) in each window and calculating a phase angle of the ridge value to obtain a wrapping phase; unpacking the wrapped phase obtained in the step 4 by using a discrete cosine phase unwrapping algorithm to obtain a continuous phase, as shown in fig. 9;

step 6: discrete cosine unwrapping processing is carried out to obtain a continuous phase and reconstruct a surface shape; directly performing WFT technical processing on the initial interference fringe image, selecting a small-size window sigma of 5 pixels for phase extraction, and obtaining wrapped phase distribution as shown in FIG. 10; comparing the wrapping phase extracted by the invention, the wrapping phase obtained by the traditional WFT technology has the phenomenon that the phase changes in certain areas are not obvious and continuous, and part of low-frequency components are considered to be interfered by background light intensity; the three-dimensional surface profile shown in fig. 11 is obtained by reconstructing the surface profile using Zernike polynomials.

The wrapping phases extracted in the traditional WFT technology all contain carrier frequency components, and artificial inclination elimination processing is needed during subsequent surface shape fitting.

The invention solves the problem of background light intensity interference of the WFT phase extraction technology under the condition of selecting a small-size window, and uniformly inhibits the linear phase error of carrier frequency interference fringes under the combined action of the optimal window size, so that the phase extraction precision under the FT-WFT combined algorithm processing is superior to that of the traditional WFT technology.

The above description is only an embodiment of the present invention, and not intended to limit the scope of the present invention, and all modifications of equivalent structures and equivalent processes, which are made by using the contents of the present specification and the accompanying drawings, or directly or indirectly applied to other related technical fields, are included in the scope of the present invention.