## Abstract

Recent development of parallel phase-shifting interferometry (PPSI) enables accurate measurement of time-varying phase maps. By combining a high-speed camera with PPSI, it became possible to observe not only time-varying but also fast phenomena including fluid flow and sound in air. In such observation, one has to remove static phase (time-invariant or slowly-varying phase unrelated to the phenomena of interest) from the observed phase maps. Ordinarily, a signal processing method for eliminating the static phase is utilized after phase unwrapping to avoid the 2*π* discontinuity which can be a source of error. In this paper, it is shown that such phase unwrapping is not necessary for the high-speed observation, and a time-directional filtering method is proposed for removing the static phase directly from the wrapped phase without performing phase unwrapping. In addition, experimental results of simultaneously visualizing flow and sound with 42 000 fps are shown to illustrate how the time-directional filtering changes the appearance. A MATLAB code is included within the paper (also in https://goo.gl/N4wzdp) for aiding the understanding of the proposed method.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Since investigation of time-varying phenomena (including heat, flow and vibration) is important in various fields of science and engineering, a lot of optical imaging methodologies has been developed to observe them [1–17]. Parallel phase-shifting interferometry (PPSI) [18–20] is one of the recent methods which enables quantitative observation of time-varying phenomena [21–31]. The authors have been applied PPSI with a polarized high-speed camera [32, 33] to visualize sound in air [34–36] and been developing signal processing techniques related to that [37–44]. Recently, it was shown that multiple physical phenomena can be simultaneously visualized by PPSI [45], which indicates further applicability of PPSI for observing complicated phenomena.

In the observation of time-varying phenomena with PPSI, static phase often deteriorates the visibility of the phenomena of interest. In this paper, static phase refers to time-invariant (or slowly varying) phase which is not related to the phenomena of interest. As usual in interferometric measurement, the observed phase is wrapped into [−*π, π*) owing to the ambiguity of absolute phase. Then, static phase results in a discontinuous pattern which obstacles the phenomena to be visualized. This pattern can be quite problematic when the phase modulation caused by the investigated phenomena is small comparing to 2*π* [34–36] as shown in the subsequent sections. Thus, a method of removing static phase is necessary for investigating such phenomena.

Ordinarily, static phase removal is accomplished via phase unwrapping [46–48]. After phase unwrapping, a smooth phase map without the discontinuity is obtained that allows the usual bias removal techniques to be executed. This process is schematically illustrated in Fig. 1, where the static phase is removed by time-directional high-pass filtering (the mathematical notations in the figure correspond to those in the subsequent sections). However, phase unwrapping may be erroneous due to the presence of noise and/or fine patterns. Its computational requirement may also be troublesome when the number of pixels in the phase map is large and/or the frame rate is high (for example, sound typically requires frame rate of more than 10 000 fps that easily produces a huge amount of data). Since every phase unwrapping method has a trade-off between computational cost and accuracy, an accurate method can be exceedingly time consuming. On the other hand, a computationally cheaper method may generate unwrapping error resulted in the 2*π* discontinuity which causes spurious step responses within the process of static phase removal that contaminates the final result. Therefore, a method of avoiding unwrapping error without computational effort is desired especially for high-speed imaging.

In this paper, a time-directional filtering method directly applicable to wrapped phase maps is proposed to overcome the above-mentioned difficulty associated with phase unwrapping. To avoid the risk of unwrapping error without computational effort, the proposed method directly filters a sequence of wrapped phase maps so that phase unwrapping is not required anymore. Such direct filtering is realized by a trick of utilizing the wrapping operator and an inverse filtering technique as described in Section 3. The proposed method is applied to simultaneous imaging of flow and sound [45], and the necessity of the time-directional filtering for such experiments is illustrated. In addition, its MATLAB code is included within the paper for easier understanding (this code with a part of the measured data can be downloaded from https://goo.gl/N4wzdp).

## 2. Observing transient phenomena by parallel phase-shifting interferometry

In this section, the problem of static phase removal is stated in the context of PPSI, and the time-difference operation is explained as the conventional static phase removal method [45]. An example of measured transient phase, utilized for testing the proposed method in Section 4, is also introduced here to demonstrate the importance of static phase removal in such measurement.

#### 2.1. Parallel phase-shifting interferometry (PPSI)

As well known, phase-shifting interferometry (PSI) is the standard method for accurately obtaining a phase map *ϕ* from multiple interference patterns. Phase-shifted fringe images {*I*^{[}^{m}^{]}} are obtained by optically shifting the phase inside the interferometer:

*B*is the background illumination,

*A*is the amplitude of fringe, and

*δ*is the phase-shifting step. From a set of phase-shifted fringes {

*I*

^{[}

^{m}^{]}}, the algorithm objective phase

*φ*is retrieved by using a phase retrieving [49–51].

When PSI is considered for observing time-varying phenomena, the objective phase map *φ* is time-varying, and thus the fringe images also vary from moment to moment:

*n*is the time index (hereafter, subscript represents time). In this case, the traditional PSI which temporally shifts the phase by

*δ*cannot be used because both

_{n}*ϕ*and

_{n}*δ*cause time variation to

_{n}*I*, and separately treating them is difficult.

_{n}In order to observe time-varying phenomena by PSI, the phase shifts must occur simultaneously. PPSI enables such simultaneous shifting of phase and observing multiple phase-shifted fringe images in parallel. Therefore, after observing the time-varying phenomena as a time-sequence (i.e., video) of interferograms $\left\{{I}_{n}^{[m]}\right\}$, the time-varying phase maps *φ _{n}* can be recovered by applying a phase retrieval algorithm to each frame one-by-one. In [34–36,45], a polarized Fizeau interferometer with a polarized high-speed camera [33] was utilized to realize PPSI in high-speed so that fast-varying phenomena including sound and fluid flow can be observed.

#### 2.2. Static phase removal for observing transient phenomena by time-difference

In most of cases, observed phase maps contains static phase *σ* as follows:

*n*th time frame becomes

*φ*+

_{n}*σ*which is contaminated by the additional phase

*σ*. When |

*σ*| ≫ |

*φ*|, which is an usual situation for measurement of weak density variation as sound in air [34–36], the phenomenon of interest

_{n}*φ*may not be visible due to the presence of static phase

_{n}*σ*. Therefore, removing the static phase is necessary for improving the visibility of

*φ*.

_{n}For removing the static phase, a straightforward method is to subtract the successive phase in time-direction as described in [45]. After retrieving two successive phase maps *φ _{n}*

_{−1}+

*σ*and

*φ*+

_{n}*σ*at time indexed by

*n*− 1 and

*n*, subtraction of these two phase maps removes the static phase

*σ*and results in the phase of interest only:

*φ*as demonstrated in the next subsection. Note that the computational cost of this process is extremely cheap that is a desired property in high-speed imaging which easily produces a huge amount of data.

#### 2.3. An example of time-varying phase measured by PPSI

The measured phase of sound in air around a small whistle blown by a gas is utilized as an example of static phase removal [45]. The experimental setup and the whistle are shown in Fig. 2. The frame rate of the high-speed camera was set to 42 000 fps so that the high-speed phenomena, fluid flow and sound, can be observed at reasonable speed. Although the sound emitted from the whistle was included in the wrapped phase map, only the fluid flow is visible as shown in the right figure of Fig. 2. This is because the static phase, which appeared as the slanting pattern, totally conceals the sound wave contained in the phase map (variation of the phase related to the sound is in the order of 0.1 rad).

For quantitatively demonstrating the effect of the static phase, variation of the phase in time direction at several pixels are shown in Fig. 3. The positions of selected three pixels are illustrated as red dots on the left of the figure, while the variation of the phase in time at these pixels are shown on the right. The small periodic waves contained in these pixels are the sound-related components which were vertically shifted by the static phase. Note that the amount of the shift is different for each pixel because the static phase *σ*(*x, y*) varies in accordance with the position of the pixel. This dependency on the position means that the spatial pattern of emitted sound wave cannot be seen by adjusting the color of the wrapped phase map because a color range suitable for one pixel (say, ±1 rad which makes the sound component visible at Pixel B) results in saturation at other pixels (such as Pixel A and Pixel C). Therefore, static phase must be removed for revealing the sound-related components contained in the phase maps.

As introduced in the previous subsection, one of the most simple methods for static phase removal is the time-directional difference in Eq. (4). The effect of taking difference of successive phase maps is illustrated in Fig. 4 (the equations in this figure will be used in the next section). As in Eq. (4), the static phase was canceled out by taking the difference in time direction, and the spatial distribution of the sound wave emitted from the whistle became clearly visible. This result demonstrates the necessity of static phase removal for observing transient phenomena especially when the phase values corresponding to the phenomena of interest is small comparing to 2*π*. However, time difference has two disadvantages (insufficiency of reduction of slowly varying components and magnification of rapidly varying components) which will be introduced in the next section. Although understanding of these disadvantages is important for interpreting the processed results, theoretical consideration of them has not been presented previously.

## 3. Proposed time-directional filtering method

In this section, we propose a time-directional filtering method whose computational cost is extremely cheap thanks to its direct applicability to the wrapped phase maps, i.e., phase unwrapping is not required. The assumption behind the proposed method are the following: (1) observation speed (frame rate) is sufficiently high comparing to the observing phenomena; and (2) time-invariant phase is required to be eliminated perfectly (time-invariant phase must not remain in the processed result). Before introducing the proposed method, the time-difference operation is firstly related with a high-pass filter, and then the wrapping operation is related with the filter so that the logical connection to the proposed method becomes apparent.

#### 3.1. Filtering of general time sequences

Here, before proceeding to some specific discussion, filtering as a general topic is briefly reviewed (readers who are already familiar with filtering can skip this subsection). Let ${\psi}_{n}^{\text{Fast}}$ be a fast varying component, and let ${\psi}_{n}^{\text{Slow}}$ be a slowly varying component. When the rates of time variation, or frequencies, of these two components differ sufficiently, then each component can be recovered from their mixture ${\psi}_{n}={\psi}_{n}^{\text{Fast}}+{\psi}_{n}^{\text{Slow}}$ by using a suitable filter.

In the situation of high-speed imaging using PPSI as in Fig. 3, the fast varying component ${\psi}_{n}^{\text{Fast}}$ corresponds to the phenomenon of interest, and ${\psi}_{n}^{\text{Slow}}$ corresponds to the static phase. Then, the problem addressed in this paper is to recover the phenomenon of interest ${\psi}_{n}^{\text{Fast}}$ from the observed mixed signal *ψ _{n}*. This situation is schematically illustrated in Fig. 5 (a), where the phenomenon of interest is the fast varying but small component, and larger static phase (consisting of the constant bias and slowly varying component) in the observed mixed signal reduces the visibility of the phenomenon of interest as in Fig. 5 (b) (see Fig. 3 for comparison with actual data).

Although separating each component is not easy in the time domain, it can be quite easy in the frequency domain when the frequencies of the components are sufficiently different. A frequency domain representation of the mixed signal in Fig. 5 (b) is shown in Fig. 5 (c). Since the frequencies of each component in Fig. 5 (c) are different, the static components can be eliminated by multiplying zero to them. Such process of multiplying a transfer function Λ in the frequency domain (or *z*-domain) is called filter:

*z*is a complex number, and Ψ(

*z*) is the

*z*-transform of

*ψ*(see [52–54] for details), (one may just consider

_{n}*z*-transform as a variant of the Fourier transform because Eq. (6) includes the Fourier transform as a special case, $z={e}^{\dot{\mathbb{I}}\omega}$, and the discussion in this paper does not rely on the property of

*z*-transform much). By multiplying a high-pass filter as in Fig. 5 (c), the static components are eliminated, and only the phenomenon of interest remains. Then, the separated signal can be recovered as in Fig. 5 (d) by inverting the filtered signal into time domain.

A filtering process can be directly written in time domain by convolution. Let *λ _{n}* be the time-domain representation (impulse response) of the transfer function Λ(

*z*). Then, Eq. (5) can be written in time domain as

#### 3.2. Time-directional difference operation as a high-pass filter

As shown in Eq. (4) and Fig. 4, the time-directional difference operation $\mathcal{D}$ can remove the static phase and improve visibility of the transient phenomena. Since $\mathcal{D}$ is the operation of generating the sequence {*φ _{n}* −

*φ*

_{n}_{−1}} from {

*φ*}, it can be written by the convolution in Eq. (7) as

_{n}*h*is the impulse response corresponding to the time-directional difference whose non-zero elements are

*h*

_{0}= +1 and

*h*

_{1}= −1, the last equality holds from the definition of convolution in Eq. (7), and (

*h*∗

*φ*)

*is shortened to*

_{n}*h*∗

*φ*. As convolution is the time-domain representation of filtering, its frequency-domain representation provides the characteristics of the time-difference operation as a filter.

_{n}By taking *z*-transform, the transfer function *H* of the impulse response *h* is obtained as

*z*-plane) is where $\dot{\mathbb{I}}=\sqrt{-1}$, and

*ω*is the normalized frequency (

*ω*=

*π*and

*ω*= 2

*π*correspond to the Nyquist and sampling frequencies, respectively) [52–54]. This transfer function, which reveals the amount of magnification imposed by the time difference $\mathcal{D}$ for each frequency, is illustrated in Fig. 5 (c). Since |sin(0)| = 0, the time-difference operation eliminates the time-invariant phase perfectly. It can also be seen that the operation reduces the amount of slowly varying phase corresponding to small |

*ω*| and amplify the varying phase whose frequency is

*π*/3 < ω < 5

*π*/3 because 2 |sin(

*π*/6)| = 2 |sin(5

*π*/6)| = 1. That is, every component whose frequency is not

*π*/3 is attenuated or magnified by taking the difference.

This characteristics results in two unwanted effects: (1) insufficient attenuation of slowly varying components; and (2) modification of the magnitude of the phenomenon of interest (i.e., quantitative value of measured phase is distorted), which can be seen from Fig. 5 (d). For circumventing these effects, one should apply a better high-pass filter $\tilde{h}$, which attenuate the slowly varying components effectively and retain the magnitude of the phenomenon of interest as in Fig. 5 (c), instead of just difference *h*. However, a general filter $\tilde{h}$ requires a phase unwrapping process which might be exceedingly time-consuming for high-speed imaging.

#### 3.3. Calculating time-difference without phase unwrapping

In the above discussions in Sections 2.2 and 3.2, the effect of wrapping and unwrapping of the phase have been omitted for the sake of simplicity. It must be considered carefully because a general filter cannot be applied to a wrapped phase directly.

As every phase has 2*π* ambiguity, it is usually treated by the wrapping operator $\mathcal{W}$ defined by

*π, π*) by adding 2

*πk*, where

*k*is a specific integer depending on the value of

*φ*. That is, any phase map

_{n}*φ*takes a value within [−

_{n}*π, π*) by the wrapping operation $\mathcal{W}[{\phi}_{n}]$ owing to modulo of 2

*π*. Phase unwrapping is a complicated process which recovers the original phase

*φ*only from its wrapped counterpart $\mathcal{W}[{\phi}_{n}]$ up to a global constant. After unwrapping, Eqs. (4) and (8) can be directly executed. However, every phase unwrapping method is either computationally demanding or inaccurate under the presence of observation noise, and therefore avoiding the unwrapping process is preferable whenever it is possible.

_{n}Fortunately, the time difference $\mathcal{D}$ can completely avoid phase unwrapping when the frame rate of the observation is sufficiently high comparing to the phenomena of interest. More specifically, the observed phase is assumed to satisfy $-\pi \le \mathcal{D}{\phi}_{n}<\pi $ (for all *x* and *y*) which is the well known condition called Itoh’s condition [55]. Then, the following identity holds:

*n*except for the 2

*π*discontinuity; and (2) $\left|\mathcal{D}\mathcal{W}[{\phi}_{n}]\ge \right|\pi $ at the discontinuity which can be eliminated by the another wrapping operation [55]. Therefore, phase unwrapping is unnecessary for the static phase removal using time difference whenever Itoh’s condition holds in time direction. Note that this condition can be met in general by increasing the frame rate until it becomes sufficient.

However, the above trick for avoiding phase unwrapping is only valid for time difference *h*, and a general filter $\tilde{h}$ cannot be applied to wrapped phase maps directly. This can be seen as a trade-off between the computational complexity and the effectiveness of filtering because one has to compromise on the characteristics of a filter for avoiding phase unwrapping, and vice versa. Therefore, a method for applying a good filter $\tilde{h}$ without phase unwrapping is desired.

#### 3.4. Proposed method: Nonlinear filtering of wrapped phase without unwrapping

The goal of the proposed method is to apply an user-defined filter $\tilde{h}$, which eliminates the static phase without changing the information related to the phenomena of interest, to the unwrapped phase *φ _{n}* without phase unwrapping. In other words, the goal is to calculate $\tilde{h}*{\phi}_{n}$ from $\tilde{h}$ and $\mathcal{W}[{\phi}_{n}]$ without applying an unwrapping process to $\mathcal{W}[{\phi}_{n}]$. To do so, the user-defined filter $\tilde{h}$ is assumed to eliminate the zero-frequency component completely.

Then, the following combination of the filters and the wrapping operator $\mathcal{F}$ is proposed:

where*h*

^{−1}is the impulse response of the “inverse filter” of

*h*whose transfer function is This operation $\mathcal{F}$ is the proposed filter which is nonlinear because of the composition of the wrapping operator $\mathcal{W}$. By inputting a wrapped phase $\mathcal{W}[{\phi}_{n}]$ into the proposed nonlinear filter $\mathcal{F}$, it can be confirmed that the desired result $\tilde{h}*{\phi}_{n}$ is obtained as follows:

*h*∗ is canceled by

*h*

^{−1}∗, and the desired filter $\tilde{h}$ is applied without the complication associating with the wrapping effect thanks to Eq. (12).

Although the idea of utilizing *h*^{−1}∗ after Eq. (12) for avoiding the effect of *h*∗ is simple, its calculation has a serious issue caused by the numerical instability. Since the magnitude of the transfer function of the “inverse filter” in Eq. (14) diverges to the infinity at the zero frequency (*ω* = 0, or *z* = 1), it is not stable (thus, the inverse filter of *h* is usually treated as it does not exist) [52–54]. This instability of the inverse filter is caused by the existence of the pole (diverging point) at *z* = 1. Therefore, a special treatment for circumventing it is necessary to execute Eq. (15) stably.

#### 3.5. Numerically stable form of the proposed filter

The difficulty of the proposed method is caused by the pole of *H*^{−1}(*z*) at *z* = 1(*ω* = 0), i.e., any numerical error, such as the quantization error, is rapidly expanded by *H*^{−1} toward infinity. Although the easiest treatment to this issue is to move the pole toward the origin that can solve the instability, such strategy is not considered in this paper because it changes *H*^{−1}, which makes Eq. (16) different from unity, and makes design of $\tilde{h}$ difficult (in such case, $\tilde{h}$ must be designed by carefully taking the amount of change in Eq. (16) into account).

When using a computer, almost every floating-point calculation causes numerical error. Then, *h*^{−1}∗ should not be applied solely even though *h*∗, which cancels the expansion effect of *h*^{−1} theoretically, is applied in advance. On the other hand, *h*∗ must be applied solely for using Eq. (12). That is, only $\tilde{h}\ast $ can be controlled for the stabilization purpose of the proposed filter. Therefore, a design constraint is imposed upon $\tilde{h}\ast $ so that the instability of *h*^{−1}∗ is theoretically eliminated.

As in the previous subsection, the user-defined filter $\tilde{h}$ is assumed to eliminate zero-frequency component completely (note that this assumption is not restrictive because the purpose of the filter is to remove the static component). This is the design constraint which ensures that $\tilde{h}$ has a zero at *z* = 1. This also means that $\tilde{h}$ admits the following decomposition:

*z*= 1. Therefore, the design constraint ($\tilde{h}$ must have zero at

*z*= 1) eliminates the instability of $\tilde{h}*{h}^{-1}*$ as whenever $\tilde{h}$ is stable. This corresponds to the pole-zero cancellation where the unstable pole at

*z*= 1 is eliminated by the zero at the same position [52–54].

Based on the above equations, as long as $\tilde{h}$ is stable and has a zero at *z* = 1, the proposed filter $\mathcal{F}$ in Eq. (13) can be calculated by

*h*, and then Eq. (12) is applied to

*h*for circumventing the wrapping effect. In any case, the proposed method should always be implemented in the form of Eq. (19) instead of Eq. (13).

#### 3.6. Summary and MATLAB code of the proposed filtering without unwrapping

In summary, the proposed method can be written as the following procedure:

- Input the coefficients of an user-defined filter $\tilde{h}$ which has a zero at
*z*= 1, - Decompose $\tilde{h}$ into $\stackrel{\u2323}{h}$ and
*h*by factorizing the filter coefficients, - Apply
*h*∗ and the wrapping operator $\mathcal{W}$ to the wrapped phase to be filtered, - Apply $\stackrel{\u2323}{h}*$ to the phase obtained by the previous step.

*z*= 1 by subtracting the mean value of the coefficients corresponding to the numerator of the transfer function. When $\stackrel{\u2323}{h}$ is an infinite impulse response (IIR) filter, representing it by a serially-cascaded second-order (or biquad) sections with an appropriate pole-zero paring may improve the numerical stability [52–54].

An example of the implementation of the proposed method in MATLAB language is shown in Fig. 6. The input variable <m>wrappedPhase</m> is assumed to be a three-dimensional array whose size is *N*_{v} × *N*_{h} × *N*_{t}, where *N*_{v} and *N*_{h} are respectively the vertical and horizontal number of pixels, and *N*_{t} is the number of time sequence. That is, time-directional filtering is implemented as a filter applied to the third dimension. The filter coefficients b and a (corresponding to the coefficients of numerator and denominator of the transfer function, respectively) can be obtained from filter designing functions, e.g., <m>butter</m> and <m>cheby2</m>. When the filter only has the numerator coefficients b, which means the filter is finite impulse response (FIR), the denominator coefficient a should be a scalar and set to 1 (*a* = 1). If the filter coefficients are in sos form or <m>digitalFilter</m> class because of using <m>filterDesigner</m> or <m>designfilt</m>, the corresponding coefficients can be obtained by <m>sos2tf</m> or <m>tf</m>, respectively. Another way of treating them with less computational error is to directly convert into the zero-pole-gain form instead of using <m>tf2zp</m> in the 11th line of Fig. 6.

Since any unwrapping process is not required, the proposed method is a quite computationally efficient procedure because it consists of only two kinds of filtering with the wrapping operation. The increase in the amount of the computational cost per sample by using the proposed method (comparing to the ordinary filtering method using $\tilde{h}$ illustrated in Fig. 1) is one subtraction in the time-difference filter *h* and three multiplications (divisions) with three additions (subtractions) in the wrapping operation $\mathcal{W}$ because <m>mod(a,n) = a − n._{*}floor(a./n)</m>. On the other hand, the decrease in the amount of the computational cost per sample is the whole unwrapping process whose computation usually consists of a massive amount of additions and multiplications [46–48]. Therefore, the proposed method can greatly reduce the total amount of computational cost by avoiding the unwrapping process. In addition, since its procedure can be applied to each pixel independently, the proposed method can be speeded up by parallel implementation when a multi-core processor is available.

## 4. Experimental results

In this section, filtered results of the observed wrapped phase introduced in Section 2.3 are shown for illustrating how the filters with different magnitude responses change the appearance of the observed transient phenomena. It is also demonstrated that the time-directional filter can assist in interpreting the observed phenomena especially when two or more phenomena such as flow and sound are simultaneously observed.

#### 4.1. Frequency characteristics of the observed data and filters

For eliminating the static phase, high-pass filters with different cut-off frequencies and a band-pass filter whose center frequency coincide with that of the phenomenon of interest (sound) were utilized in this experiment. The magnitude responses of the filters are shown in the top of Fig. 7. The colored lines correspond to the third-order maximally flat magnitude (or Butterworth) filters, while the gray line represents the time difference *h*, and the black line is the response without filtering. For the high-pass filters, the cut-off frequencies *f*_{c} are also indicated. The band-pass filter was centered at 8700 Hz, which is about the frequency of observed sound, and the bandwidth was set to 2000 Hz.

These filters were chosen subjectively based on the frequency characteristics of the data to be processed. The spectrum of observed phase obtained by PPSI is shown in Fig. 7. Spectra of filtered results are also indicated with the colors which are identical with those of the filter responses. These spectra were calculated by taking the Fourier transform in time direction for each pixel and summing their absolute values up, where the magnitude was normalized with the maximum value of the unwrapped phase without filtering. In these spectra, both flow components and sound components are mixed together with the static phase and measurement noise as the observed phase contains both the gas flow and sound (see Fig. 4). Therefore, for visualizing the sound, a filter which retains the sound component and eliminates the other components should be used.

The sound component in these spectra appears as a single peak at the resonant frequency of the whistle (around 8700 Hz) because the sound emitted from the whistle in Fig. 2 was sinusoidal. An appropriate filter targeting the sound should have a flat response around the frequency of sound and become nearly zero at the other frequencies (see Fig. 5 for a schematic example). The maximally flat magnitude filters were chosen here because their responses are flat. Then, the cut-off frequencies are the parameter of the filters that affect the visibility of the phenomena. The cut-off frequencies of the high-pass filters were chosen so that the amount of attenuation of low-frequency components becomes like a step-by-step manner as show in the spectra of Fig. 7 (the cut-off frequency was manually increased so that the difference of the area of attenuated components became similar). The band width of the band-pass filter was chosen in a similar manner.

#### 4.2. Filtered results of the wrapped phase maps

For illustrating the effects of the filters, the appearances of the filtered results obtained by the proposed method (spectra in Fig. 7) are shown in two different ways: adjusted color and constant color. Firstly, the appearances with automatically adjusted color range, calculated from the root mean square (RMS) values of the filtered results, are shown in Fig. 8 (i.e., each range of the color was determined by the averaged level of the pixel values of each video). Five consecutive instances of the filtered results are shown as a representative of whole video because the frequency of sound was about 8700 Hz whose period is 4.86 frames for 42 000 fps. Each row of the figure represents (from top to bottom) the observed wrapped phase, the high-pass filtered results with different cut-off frequencies and the band-pass filtered result, where the time evolves from left to right at 23.81 *µ*s interval. Here, the well-known zero-phase filtering technique [52–54], which can be found as <m>filtfilt</m> function in MATLAB, was utilized for avoiding the group delay (the magnitude responses in Fig. 7 includes this effect). An executable code and a part of the measured data can be downloaded from https://goo.gl/N4wzdp.

Although the gas flow can be clearly recognized, the sound cannot be seen in the wrapped phase because of the smallness of its magnitude as in Figs. 3 and 7. Therefore, the time-directional filtering was necessary for revealing the sound from the wrapped phase. As in the figures, the high-pass filters removed the static phase (the slanting pattern of the wrapped phase) and remained the time-varying phenomena of interest. Note that although phase unwrapping was not performed in advance, the proposed method successfully eliminated the 2*π* discontinuity (artificial noise related to the slanting pattern cannot be seen in the figures).

Since the order of magnitudes of the flow and sound was different as in Fig. 7, a color range suitable for both flow and sound did not exist (sound is not visible if one set the color suitable for flow, or flow may be totally saturated for the color suitable for sound). In such situation, time-directional filtering can be utilized to balance the magnitudes of the time-varying phenomena so that both phenomena are visualized simultaneously. By removing the low-frequency components using high-pass or band-pass filters, the level of the gas flow was reduced and the sound became more apparent, where the amount of reduction of the flow components was determined by the filter responses in Fig. 7. The band-pass filter eliminated the flow components more than the high-pass filters because it removed not only low- but also high-frequency components. Although completely separating the two phenomena by only a filter was difficult, these filters enabled simultaneous observation of flow and sound from multiple points of view by changing the filter responses. Changing the responses is easy for the proposed method because any stable filter $\tilde{h}$ having a zero at *z* = 1 can be used. Moreover, applying several filters to a single data is effortless thanks to the cheap computational cost of the proposed method realized by the direct filtering of wrapped phase maps without unwrapping.

For confirming the appearances of the sound in each setting, the same figures with fixed color range are shown in Fig. 9. The appearances of the sound are almost the same because all filters used in this section had flat responses at the frequency of the sound as in Fig. 7. This result indicates the advantage of utilizing such flat magnitude filter and the importance of possibility of utilizing a general filter that is the key point of the proposed method. Note that the appearance for the case *f*_{c} = 200 Hz is different from the others. This is because some slowly varying components remained in the result owing to the low cut-off frequency. Since those components seems to be uniformly changing all pixels, they should be the vibration of the mirror reflecting the object light (see Fig. 2). Eliminating such slowly varying components not related to the phenomena of interest is important for evaluating the observed data.

#### 4.3. Comparison between the ordinary method and the proposed filter

For illustrating the advantage of direct filtering of the wrapped phase, the proposed method is compared with the ordinary filtering method shown in Fig. 1 which requires phase unwrapping in advance. The filtered results using the band-pass filter are shown in Fig. 10, where the unwrapping method of Goldstein et al. [46] was utilized in the ordinary method. As in the figures, the ordinary method resulted in unwrapping error which spreads along time due to the time-directional filter applied after the phase unwrapping. Filtering based on phase unwrapping always has risk of such spread of unwrapping error because phase unwrapping is a difficult problem especially when observation noise is contained in the data. On the other hand, the proposed method does not have such risk thanks to the direct filtering which avoid phase unwrapping.

## 5. Conclusions

In this paper, a time-directional filtering method which can be directly applied to wrapped phase was proposed. By considering Itoh’s condition in terms of filtering, the wrapping effect was circumvented without utilizing an explicit unwrapping process which may be time consuming. The inverse filter is combined with Itoh’s condition to eliminate the unwanted effect of the circumvention, and then an user-defined filter is applied. Instability related to the inverse filter is canceled by the zero contained in the user-defined filter, which also avoids the unnecessary computation. The effectiveness of the time-directional filtering is demonstrated by the experimental data, which consist of two time-varying phenomena (gas flow and sound), measured by PPSI. In addition, an example code in MATLAB language was provided to demonstrate how the proposed method can be implemented through the zero-pole-gain representation of the filter.

Since the proposed method allows control of the filter design (in contrast to the time difference), it can be applied to any situation that the static phase must be removed from a time sequence of wrapped phase maps. In addition, it is possible to execute the proposed filter on low-performance devices because its computational cost is almost same as the usual filtering because unwrapping is not required but just an additional modulo operation (wrapping operation) is necessary. Further consideration of the proposed method in terms of the fixed-point arithmetic, which is important for hardware implementation, is remained as the future works.

## Funding

Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Research Activity Start-up (17H07191), Grant-in-Aid for JSPS Research Fellow (16J06772).

## References and links

**1. **A. Hettwer, J. Kranz, and J. Schwider, “Three channel phase-shifting interferometer using polarization-optics and a diffraction grating,” Opt. Eng. **39**, 960–966 (2000). [CrossRef]

**2. **S.-H. Baik, S.-K. Park, C.-J. Kim, and S.-Y. Kim, “Shockwave visualization using holographic interferometer,” Opt. Rev. **7**, 535–542 (2000). [CrossRef]

**3. **G. S. Settles, *Schlieren and Shadowgraph Techniques: Visualizing Phenomena in Transparent Media* (Springer, 2001). [CrossRef]

**4. **H. S. Ko, K. Okamoto, and H. Madarame, “Reconstruction of transient three-dimensional density distributions using digital speckle tomography,” Meas. Sci. Technol. **12**, 1219 (2001). [CrossRef]

**5. **G. Pedrini, I. Alexeenko, W. Osten, and H. J. Tiziani, “Temporal phase unwrapping of digital hologram sequences,” Appl. Opt. **42**, 5846–5854 (2003). [CrossRef] [PubMed]

**6. **G. Pedrini, W. Osten, and M. E. Gusev, “High-speed digital holographic interferometry for vibration measurement,” Appl. Opt. **45**, 3456–3462 (2006). [CrossRef] [PubMed]

**7. **J.-M. Desse, “Recent contribution in color interferometry and applications to high-speed flows,” Opt. Lasers Eng. **44**, 304–320 (2006). [CrossRef]

**8. **S. Suzuki, K. Sakaue, and K. Iwanaga, “Measurement of energy release rate and energy flux of rapidly bifurcating crack in homalite 100 and Araldite B by high-speed holographic microscopy,” J. Mech. Phys. Solids **55**, 1487–1512 (2007). [CrossRef]

**9. **R. Hruschka, S. O’Byrne, and H. Kleine, “Diode-laser-based near-resonantly enhanced flow visualization in shock tunnels,” Appl. Opt. **47**, 4352–4360 (2008). [CrossRef] [PubMed]

**10. **P. Bon, G. Maucort, B. Wattellier, and S. Monneret, “Quadriwave lateral shearing interferometry for quantitative phase microscopy of living cells,” Opt. Express **17**, 13080–13094 (2009). [CrossRef] [PubMed]

**11. **W. Sun, J. Zhao, J. Di, Q. Wang, and L. Wang, “Real-time visualization of Karman vortex street in water flow field by using digital holography,” Opt. Express **17**, 20342–20348 (2009). [CrossRef] [PubMed]

**12. **D. D. Aguayo, F. M. Santoyo, M. H. D. la Torre-I, M. D. Salas-Araiza, C. Caloca-Mendez, and D. A. G. Hernandez, “Insect wing deformation measurements using high speed digital holographic interferometry,” Opt. Express **18**, 5661–5667 (2010). [CrossRef] [PubMed]

**13. **P. K. Panigrahi and K. Muralidhar, *Schlieren and Shadowgraph Methods in Heat and Mass Transfer* (Springer, 2012). [CrossRef]

**14. **S. M. noz Solís, F. M. Santoyo, and M. del Socorro Hernández-Montes, “3D displacement measurements of the tympanic membrane with digital holographic interferometry,” Opt. Express **20**, 5613–5621 (2012). [CrossRef]

**15. **T. Aratani, T. Inage, Y. Miwa, M. Ota, and K. Maeno, “Application of high-speed camera to 4D-CT density measurement of unsteady shock-vortex flow discharged from two inclined and cylindrical holes,” J. JSEM **13**, s64–s68 (2013).

**16. **J. S. Pérez-Huerta, T. Saucedo-Anaya, I. Moreno, D. Ariza-Flores, and B. Saucedo-Orozco, “Digital holographic interferometry applied to the investigation of ignition process,” Opt. Express **25**, 13190–13198 (2017). [CrossRef] [PubMed]

**17. **B.-C. Cui, J.-L. Wang, K.-N. Yao, and T. Chen, “Measurement of high-dynamic temperature field using high-speed quadriwave lateral shearing interferometer,” Optoelectron. Lett. **14**, 124–128 (2018). [CrossRef]

**18. **Y. Awatsuji, M. Sasada, and T. Kubota, “Parallel quasi-phase-shifting digital holography,” Appl. Phys. Lett. **85**, 1069–1071 (2004). [CrossRef]

**19. **J. E. Millerd, N. J. Brock, J. B. Hayes, M. B. North-Morris, M. Novak, and J. C. Wyant, “Pixelated phase-mask dynamic interferometer,” Proc. SPIE **5531**, 304–314 (2004). [CrossRef]

**20. **M. Novak, J. Millerd, N. Brock, M. North-Morris, J. Hayes, and J. Wyant, “Analysis of a micropolarizer array-based simultaneous phase-shifting interferometer,” Appl. Opt. **44**, 6861–6868 (2005). [CrossRef] [PubMed]

**21. **T. Tahara, K. Ito, T. Kakue, M. Fujii, Y. Shimozato, Y. Awatsuji, K. Nishio, S. Ura, T. Kubota, and O. Matoba, “Parallel phase-shifting digital holographic microscopy,” Biomed. Opt. Express **1**, 610–616 (2010). [CrossRef]

**22. **T. Kakue, R. Yonesaka, T. Tahara, Y. Awatsuji, K. Nishio, S. Ura, T. Kubota, and O. Matoba, “High-speed phase imaging by parallel phase-shifting digital holography,” Opt. Lett. **36**, 4131–4133 (2011). [CrossRef] [PubMed]

**23. **O. Matoba, H. Inokuchi, K. Nitta, and Y. Awatsuji, “Optical voice recorder by off-axis digital holography,” Opt. Lett. **39**, 6549–6552 (2014). [CrossRef] [PubMed]

**24. **P. Xia, Y. Awatsuji, K. Nishio, and O. Matoba, “One million fps digital holography,” Electron. Lett. **50**, 1693–1695 (2014). [CrossRef]

**25. **T. Kakue, Y. Endo, T. Nishitsuji, T. Shimobaba, N. Masuda, and T. Ito, “Digital holographic high-speed 3D imaging for the vibrometry of fast-occurring phenomena,” Sci. Rep. **7**, 10413 (2017). [CrossRef] [PubMed]

**26. **T. Fukuda, Y. Wang, P. Xia, Y. Awatsuji, T. Kakue, K. Nishio, and O. Matoba, “Three-dimensional imaging of distribution of refractive index by parallel phase-shifting digital holography using Abel inversion,” Opt. Express **25**, 18066–18071 (2017). [CrossRef] [PubMed]

**27. **D. I. Serrano-García, A. Martínez-García, N.-I. Toto-Arellano, and Y. Otani, “Dynamic temperature field measurements using a polarization phase-shifting technique,” Opt. Eng. **53**, 112202 (2014). [CrossRef]

**28. **E. Shoji, R. Nakaoku, A. Komiya, J. Okajima, and S. Maruyama, “Quantitative visualization of boundary layers by developing quasi-common-path phase-shifting interferometer,” Exp. Therm. Fluid Sci. **60**, 231–240 (2015). [CrossRef]

**29. **E. Shoji, A. Komiya, J. Okajima, H. Kawamura, and S. Maruyama, “High-speed phase-shifting interferometry using triangular prism for time-resolved temperature measurement,” Appl. Opt. **54**, 6297–6304 (2015). [CrossRef] [PubMed]

**30. **A. Safrani and I. Abdulhalim, “Full-field parallel interferometry coherence probe microscope for high-speed optical metrology,” Appl. Opt. **54**, 5083–5087 (2015). [CrossRef] [PubMed]

**31. **M. Ney, A. Safrani, and I. Abdulhalim, “Three wavelengths parallel phase-shift interferometry for real-time focus tracking and vibration measurement,” Opt. Lett. **42**, 719–722 (2017). [CrossRef] [PubMed]

**32. **T. Onuma and Y. Otani, “A development of two-dimensional birefringence distribution measurement system with a sampling rate of 1.3 MHz,” Opt. Commun. **315**, 69–73 (2014). [CrossRef]

**33. **T. Yatagai, B. J. Jackin, A. Ono, K. Kiyohara, M. Noguchi, M. Yoshii, M. Kiyohara, H. Niwa, K. Ikuo, and T. Onuma, “Instantaneous phase-shifting fizeau interferometry with high-speed pixelated phase-mask camera,” Proc. SPIE **9660**, 966018 (2015).

**34. **K. Ishikawa, K. Yatabe, N. Chitanont, Y. Ikeda, Y. Oikawa, T. Onuma, H. Niwa, and M. Yoshii, “High-speed imaging of sound using parallel phase-shifting interferometry,” Opt. Express **24**, 12922–12932 (2016). [CrossRef] [PubMed]

**35. **K. Ishikawa, K. Yatabe, Y. Ikeda, Y. Oikawa, T. Onuma, H. Niwa, and M. Yoshii, “Interferometric imaging of acoustical phenomena using high-speed polarization camera and 4-step parallel phase-shifting technique,” Proc. SPIE **10328**, 103280I (2017).

**36. **Y. Oikawa, K. Yatabe, K. Ishikawa, and Y. Ikeda, “Optical sound field measurement and imaging using laser and high-speed camera,” in “Proc. 45th Int. Congr. Noise Control Eng. (INTER-NOISE 2016),” (2016), pp. 258–266.

**37. **K. Yatabe and Y. Oikawa, “PDE-based interpolation method for optically visualized sound field,” in “IEEE Int. Conf. Acoust., Speech Signal Process. (ICASSP),” (2014), pp. 4738–4742.

**38. **K. Yatabe and Y. Oikawa, “Optically visualized sound field reconstruction based on sparse selection of point sound sources,” in “IEEE Int. Conf. Acoust., Speech Signal Process. (ICASSP),” (2015), pp. 504–508.

**39. **K. Yatabe and Y. Oikawa, “Optically visualized sound field reconstruction using Kirchhoff–Helmholtz equation,” Acoust. Sci. Tech. **36**, 351–354 (2015). [CrossRef]

**40. **N. Chitanont, K. Yaginuma, K. Yatabe, and Y. Oikawa, “Visualization of sound field by means of Schlieren method with spatio-temporal filtering,” in “IEEE Int. Conf. Acoust., Speech Signal Process. (ICASSP),” (2015), pp. 509–513.

**41. **N. Chitanont, K. Yatabe, K. Ishikawa, and Y. Oikawa, “Spatio-temporal filter bank for visualizing audible sound field by schlieren method,” Appl. Acoust. **115**, 109–120 (2017). [CrossRef]

**42. **K. Yatabe and Y. Oikawa, “Convex optimization based windowed Fourier filtering with multiple windows for wrapped phase denoising,” Appl. Opt. , **55**, 4632–4641 (2016). [CrossRef] [PubMed]

**43. **K. Yatabe, K. Ishikawa, and Y. Oikawa, “Compensation of fringe distortion for phase-shifting three-dimensional shape measurement by inverse map estimation,” Appl. Opt. **55**, 6017–6024 (2016). [CrossRef] [PubMed]

**44. **K. Yatabe, K. Ishikawa, and Y. Oikawa, “Acousto-optic back-projection: Physical-model-based sound field reconstruction from optical projections,” J. Sound Vib. **394**, 171–184 (2017). [CrossRef]

**45. **K. Ishikawa, R. Tanigawa, K. Yatabe, Y. Oikawa, T. Onuma, and H. Niwa, “Simultaneous imaging of flow and sound using high-speed parallel phase-shifting interferometry,” Opt. Lett. **43**, 991–994 (2018). [CrossRef] [PubMed]

**46. **D. C. Ghiglia and M. D. Pritt, *Two-Dimensional Phase Unwrapping: Theory, Algorithms and Software* (Wiley, 1998).

**47. **Q. Kemao, S. H. Soon, and A. Asundi, “A simple phase unwrapping approach based on filtering by windowed fourier transform,” Opt. Laser Technol. **37**, 458–462 (2005). [CrossRef]

**48. **Q. Kemao, W. Gao, and H. Wang, “Windowed Fourier-filtered and quality-guided phase-unwrapping algorithm,” Appl. Opt. **47**, 5420–5428 (2008). [CrossRef] [PubMed]

**49. **K. Yatabe, K. Ishikawa, and Y. Oikawa, “Improving principal component analysis based phase extraction method for phase-shifting interferometry by integrating spatial information,” Opt. Express **24**, 22881–22891 (2016). [CrossRef] [PubMed]

**50. **K. Yatabe, K. Ishikawa, and Y. Oikawa, “Simple, flexible, and accurate phase retrieval method for generalized phase-shifting interferometry,” J. Opt. Soc. Am. A **34**, 87–96 (2017). [CrossRef]

**51. **K. Yatabe, K. Ishikawa, and Y. Oikawa, “Hyper ellipse fitting in subspace method for phase-shifting interferometry: practical implementation with automatic pixel selection,” Opt. Express **25**, 29401–29416 (2017). [CrossRef]

**52. **L. B. Jackson, *Digital Filters and Signal Processing: With MATLAB Exercises* (Springer, 1996). [CrossRef]

**53. **P. O’Shea, A. Z. Sadik, and Z. M. Hussain, *Digital Signal Processing: An Introduction with MATLAB and Applications* (Springer, 2011). [CrossRef]

**54. **S. M. Alessio, *Digital Signal Processing and Spectral Analysis for Scientists: Concepts and Applications* (Springer, 2016).

**55. **K. Itoh, “Analysis of the phase unwrapping algorithm,” Appl. Opt. **21**, 2470 (1982). [CrossRef] [PubMed]