Phase response reconstruction for non-minimum phase systems using frequency-domain magnitude values

It is often more complicated to measure the phase response of a large system than the magnitude. In that case, one can attempt to use the Kramers-Kronig (KK) relations for magnitude and phase, which relates magnitude and phase analytically. The advantage is that then only the magnitude of the frequency response needs to be measured. We show that the KK relations for magnitude and phase may yield invalid results when the transfer function has zeros located in the right half of the complex s-plane, i.e. the system is non-minimum phase. In this paper we propose a method which enables to determine these zeros, by using specific excitation signals and measuring the resulting time responses of the system. The method is verified using blind tests among the authors. When the locations of the zeros in the right half of the complex s-plane are known, modified KK relations can be successfully applied to non-minimum phase systems. We demonstrate this by computing the phase response of the electric field, excited by a point dipole source inside a closed cavity with Perfect Electrically Conducting (PEC) walls. Also, the effects of simulated measurement noise are considered for this example.


Introduction
The frequency response of a system is an important property to define the relation between the input and output and is used for characterizing and designing systems. The frequency response generally consists of a magnitude and a phase. This is under the assumption that the system behaves (approximately) Linear Time Invariant (LTI).
Applications where measuring the frequency response is important are among others in the Defence industry. An example is the evaluation of navy ships for hardening against the Nuclear Electro Magnetic Pulse (NEMP) threat. The frequency response of the electric field transmission, from a RF source outside, is measured at key positions inside the ship. Using the measured frequency response, a full-threat (N)EMP can then be numerically simulated in the time domain by the Inverse Fast Fourier Transform (IFFT). This has considerable advantages in terms of costs, complexity and safety measures compared to generating a full-threat EMP and measuring the time response of the electric field inside the ship. A problem which emerges in this application, is that the phase of the frequency response cannot be measured easily. However, the phase information is also required to accurately reconstruct the time response of the electric field inside the ship when a (N)EMP threat is simulated. Other applications, where the phase is important but difficult to measure, can be found in the areas such as optical device characterization and biomedical imaging [1].
Fortunately, the magnitude of the frequency response can often be measured accurately. To avoid the necessity to use the phase information, one could in some cases use a minimum phase method as described in e.g. [2]. In [3] [4] [5] the KK relations for magnitude and phase are used to accurately calculate the phase from the magnitude of the frequency response. This is under the assumption that the corresponding system is minimum phase. It is equivalent to requiring that the transfer function of the system does not have zeros in the right half of the complex s-plane. However, the minimum phase method is not always applicable and some examples of non-minimum phase systems are described in [6]. One of these examples is the optical reflection coefficient of optical materials. Another example describes the frequency response of mechanical systems such as an atomic force microscope which use a flexible cantilever to probe a surface.
In [3] [5] modifications of the KK relations, using the Blaschke product, are described. These modifications allow the use of the KK relations for magnitude and phase, even when the system is non-minimum phase. However, an important assumption here is that the locations of the zeros of the transfer function in the right half of the s-plane are known. It is stressed that for many applications the locations of these zeros are unknown, and difficult to obtain.
Here we propose a method which enables us to accurately estimate the locations of the zeros in the right half of the s-plane. This is done using specific excitation signals and calculating the time responses of the system. For this method, no a priori knowledge of the transfer function is required. After using our new zero-search method to precisely locate the zeros, we then successfully reconstruct the phase from the measured magnitude of the frequency response using modified KK relations. This paper is organized as follows: Section 2 gives a detailed mathematical description of the EM cavity model, where we apply the theory described in [7] to a closed PEC cavity with a point dipole source. Section 3 discusses the KK relations for magnitude and phase and is concluded with some numerical results where the KK relations for magnitude and phase are applied to the cavity model from Section 2. In Section 4 modified KK relations using the Blaschke product are introduced. In Section 5 we describe our new zero-search method including numerical results. It is concluded with a successful application of modified KK relations to the cavity model from Section 2. In Section 6 we investigate the 2 consequences of noise on the accuracy of the proposed zerosearch method and the phase reconstruction method using the modified KK relations. An example with added noise using the cavity model from Section 2 is discussed. Finally, we present conclusions and recommendations.

Geometry
We start with defining the geometry of the problem under consideration, which we will use to test the phase reconstruction using the KK relations for magnitude and phase. The geometry consists of a closed cavity with PEC walls and enclosed current source . An example of such a cavity is a rectangular cavity with dimensions , and , is shown in Fig. 1.

Eigenmode expansion of the electric field
The boundary conditions of the electric field on the perfectly conducting cavity walls, Ω , are given by where is the inward normal vector on Ω . The charge free domain Ω inside the cavity is filled with a homogeneous isotropic time-invariant insulator. Thus, the electric field ( , ) at location inside the cavity and at time , satisfies the equation In this equation, ∇ 2 denotes the vector Laplace operator, and denote the permittivity and permeability of the insulator, respectively and ( , ) denotes the impressed electric current density inside Ω.
Since we consider the field distribution inside the cavity, it is advantageous to describe the electric field, ( , ), in terms of its truncated modal expansion [7]. For a rectangular cavity, it is convenient to make a distinction between transverse magnetic and transverse electric modes. This means that either the electric or the magnetic field is forced to be transverse to the axis in our Cartesian coordinate system. The modal expansion is split up in a transverse magnetic and a transverse electric part where denotes the number of TM modes denotes the number of TE modes used in the modal expansion, ( ) is the 3-D electric field distribution of the ν-th mode of the cavity and ( ) is the corresponding time dependent expansion coefficient. The total number of modes used is = + . For an exact representation, an infinite number of modes should be considered: → ∞. However, for numerical purposes the number of modes taken should be finite but sufficiently large so that the electric field is computed appropriately.
The 3-D electric field distributions of the modes ( ) can be found by solving the source-free 3-D Helmholtz equation in the frequency domain and enforcing the boundary conditions in (1) where 2 = 2 is the -th wavenumber corresponding to the ν-th mode which will resonate with frequency .
Here, the and mode field patterns of a rectangular cavity are calculated analytically [8]. Note, however, that the method described in this paper can also be applied to a closed cavity with any geometry, if the modes of the cavity can be found by solving (4) analytically or numerically.

Cavity transfer function
The cavity transfer function defines the frequency response of the electric field at an arbitrary location due to a local current source inside the cavity. To obtain the cavity transfer function we consider a source located at = and an observation point located at = inside the cavity and find an explicit relation between the impressed current density ( , ) and the observed field ( , ).
To keep the model simple, we assume that ( , ) in (2) is a single dipole source, oriented in the -direction as shown in Fig. 1. Furthermore, we assume that the dipole has length [m], and has a uniform current distribution. Thus, we define where ̂ is a unit vector in the x-direction, ( ) [A] the impressed current source function, which is the input current of the dipole and ( ) is the Dirac delta function. The function ( − ) is defined such that the current distribution is localised within an infinitesimal volume hence ↓ 0. In [7], it is discussed to rewrite (2) as a system of second order Ordinary Differential Equations (ODEs) by substituting (5) in (2) and using the orthogonality of the cavity mode field patterns, where indicates the energy stored in the -th mode, which is computed as where is the permittivity of the medium inside the cavity and * denotes complex conjugation. While assuming that ( ) is constant across the infinitesimal volume enclosing ( , , ), we substitute (5) into (6) which leads to Equation (8) where [ ( )] denotes the x component of the ν-th mode field pattern at the location of the source. After some mathematical steps, a commonly known state equation is obtained from (9) with input ( ), input vector , state matrix and state vector ( ) defined as Solving this equation for ( ) leads to the time dependent expansion coefficients in (3). Now ( , ), the electric field inside the cavity observed at position , can be written in terms of the state vector ( ) which we rewrite as ( , ) = ( ) ( ).
The Laplace transform of the electric field is defined as follows from which it is easily verified that ̃( , ) = ( ) ̃( ). Then (10) can be converted to the frequency domain by substituting the differential operator / with the complex frequency variable = + , β > 0 . The electric field is then written in terms of the excitation current source as with the transfer function ( , ) given by where is the identity matrix. From (14) and (15) It is observed that the transfer function in (16) has 2 poles all present in mirrored pairs. Furthermore, we have 2 − 1 zeros with one zero at = 0 and the remaining zeros present in mirrored pairs. These zeros can be determined from (16) by combining the sum of fractions such that we obtain a single fraction and then finding complex values of s for which the numerator becomes zero. Because the cavity is lossless, the poles are always located at the imaginary frequency axis, = ± . However, the zeros can have nonzero real part, i.e. ≠ 0 in = ± . Note that the locations of the zeros depend on the shape of the cavity and the location of the observation point and source, i.e. they depend on [ ( )] , [ ( )] * , and . As a result, the transfer function of the x, y and z component of ( , ) share the same poles, however their zeros may be different.

Kramers-Kronig relations for magnitude and phase
We assume the transfer function of a system has a causal time response. The causal response function in Fourier space is then explicitly given by ( ) being real implies that (− ) = * ( ). We now construct the analytical continuation of by replacing with the complex frequency variable = + , where ≥ 0. It follows that ( ) is an analytical function in the left half of the s-plane. The Cauchy integration theorem for the causal transfer function ( ) yields [9] where the closed contour consists of the imaginary axis and the infinite semicircle enclosing the right half of the s-plane.
We assume that is oriented counterclockwise. Let us also assume that ( ) → 0 for | | → ∞. Considering then again imaginary values of = lim ↓0 + which leads to the KK relations [9] Re where denotes the principal value. Thus, the real and imaginary parts are Hilbert transform [10] pairs of each other Causality yields via the KK relations the imaginary part of the frequency response if the real part is known and vice versa. Therefore, one may expect that either the magnitude or the phase is enough to reconstruct analytic frequency responses. These are defined as where ( ) is the phase. For completeness, their relation to the real and imaginary part is given by Recalling that (− ) = * ( ), an implication of ( ) = * ( ), yields the symmetry relations It is more the rule than the exception that in experiments only the magnitude can be measured. The question is whether the phase, and eventually the analytic frequency response, can be reconstructed via KK(-like) relations. The theory related to the last question is discussed in, for example, [2] where also the concept of minimum phase is introduced. Here we partly follow [3] which appears to be most useful for our purposes. The reader may consult the references for more details. In order to separate the phase, the complex logarithm of the response is studied. We introduce the infinitesimal parameter and consider where > 0 . At this point it seems to be possible to formulate dispersion relations for the logarithm of the magnitude and the phase. Indeed, with some caveats to be addressed below, the phase follows as or, explicitly showing arguments and the prescription, Note that ensures that the integration path does not contain any poles or zeros of ( ) in case they lay on the imaginary axis. The first caveat which has to be made is that the complex response does not support zeros in the right half of the s-plane. In that case the logarithm is not defined and to satisfy the KK relations the appearing logarithm should be analytic for Re > 0. The requirement that the logarithm of the magnitude and the phase form a Hilbert pair is called the minimum phase condition, see e.g. [2]. In [11] this minimum phase condition has been imposed. The second caveat is that the complex response does not support zeros at infinite complex frequency i.e. | | → ∞, ( ) → 0.

Numerical results
To test the KK relations for magnitude and phase presented in subsection 3.1 with the cavity transfer function as defined in subsection 2.3 , we present three examples here. For the different examples, we vary the location where the electric field is observed, . Also, the dimensions of the cavity are changed. For all examples the magnitude of the transfer function is used to reconstruct the phase using (26). We compare the actual phase of the transfer function to the reconstructed phase and see if the minimum phase condition is violated, which is dependent on the location of the zeros in the transfer function. Note that there is also a zero at infinite complex frequency because [ ( , )] → 0 for | | → ∞, as is observed from (16). For all examples studied here however, it has been observed that the phase error caused by the zero at infinity is negligible. In the examples we only consider the x component of the electric field, which means that we only must compute [ ( , )] . The resulting transfer functions were verified using numerical results obtained by using the Method of Moments (MoM) solver in FEKO [12].
To observe the effect of phase errors on the time response of the electric field when using the reconstructed phase, we excite the cavity fields using a pulsed current generated by a point dipole source. The time response of the current pulse is defined as with 0 = 6.5 × 10 4 V/m , 1 = 4 × 10 7 s −1 , 1 = 6 × 10 8 s −1 and = 1 × 10 −8 s . The current pulse, ( ) , is subsequently filtered through a first order high pass filter with crossover frequency of 0.6 MHz, such that the resulting current source inside the cavity, ( ), has zero mean. The pulse has a -3 dB frequency bandwidth of approximately 8 MHz and a -10 dB frequency bandwidth of approximately 20 MHz.
We then compare the time response of the electric field using the directly calculated phase to the time response of the electric field using the reconstructed phase. For the first example we set the cavity size parameters in Fig. 1 to = 0.8 m, = 0.9 m and = 1 m. The source is positioned at For all the following calculations, in (26) is set to 0.25/ MHz . This value was experimentally determined to give good results. Note that the length of the dipole source as defined in (5) only scales the transfer function by a constant factor while the phase response is not effected. For calculating [ ( , )] , all modes for which /2 > 500 MHz are discarded, which means that a total of = 32 (16 TE and 16 TM) modes are used for the calculations.
We first compare the phase directly extracted from the model with the reconstructed phase. In Fig. 2-a we see that the reconstructed phase starts to deviate for frequencies higher than 430 MHz, which is due to the fact that two zeros are present in the right half of the s-plane at = 2 × (0.055 ± 4.4) × 10 8 . These zeros violate the minimum phase condition needed for (26) to hold.
The resulting time response is shown in Fig. 2-b. We see that the reconstructed time response is very similar to the directly calculated time response. This is because the bandwidth of the pulse is relatively small, thus the cavity modes for which the reconstructed phase is starting to deviate have a much smaller magnitude than the lower frequency modes which are excited more significantly by the pulse. As a second example, we take the observation point at The reconstructed phase is shown in Fig. 3-a. Similar effects are observed compared to the previous case shown in Fig. 2. However, for this case a total of six zeros are present in the right half of the s-plane, which are now located at: = 2 × (2.3 ± 3.5) × 10 8 , = 2 × (0.061 ± 4.3) × 10 8 and = 2 × (0.31 ± 5) × 10 8 . The zeros of the transfer function are located at different points compared to the previous case, since we changed the observation point, , which also moves some of the zeros further away from the axis in the right half of the s-plane. This, combined with the fact that there are significantly more zeros in the right half of the s-plane, explains why the reconstructed phase is deviating much more severely compared to the previous case.
The resulting time response shows a much larger difference between the reconstructed and directly calculated results, as is presented in Fig. 3-b. This is as expected, because the difference between the reconstructed phase and the directly calculated phase is already in the order of 90° at the lowest order mode resonance frequency, i.e. 220 MHz. However, as the lowest order mode has the largest contribution to the resulting time responses, a simple phase correction of 90° should already significantly improve the reconstructed time response. We see that apart from this 90° phase shift, the overall shape and most importantly the envelope and peak values of the time signal are preserved. a b Fig. 3. (a) Magnitude, phase and reconstructed phase of ( ) as a function of frequency for From this we conclude that some error in the reconstructed phase might be allowed. However, this error largely depends on the locations of the zeros of the transfer function. In practice this means that we do not know what phase correction we must apply, unless we are able to directly determine the locations of the complex zeros.
As third example, we increase the size of the cavity and change its shape by setting = 8 m, = 9 m and = 3 m. The location of the source and the observation point remain unchanged from the second example. Note that because the cavity is much larger now, the modes have much larger wavelengths and as a result, the corresponding resonance frequencies are much lower. Thus, more modes are excited by the EMP and we need to consider a larger number of modes in the mode expansion. To bound computation time, we limit the number of modes used in the mode expansion to 50 TE-and 50 TM-modes, hence resulting in a total of = 100 modes. The results for the configuration are shown in Fig. 4. Note that the transfer function is now plotted between 40 MHz and 130 MHz.
In Fig. 4-a we see that the discrepancy between the exactly calculated and reconstructed phase is now more severe compared to the cases shown in Fig. 2-a and Fig. 3

-a.
This is caused by the fact that now a total number of 16 zeros are located in the right half of the s-plane compared to only 2 zeros for the example shown to Fig. 2 and 6 zeros for the example shown in Fig. 3. a b Fig. 4. (a)  Next, we calculate the resulting time response of the electric field ( , ). The result is shown in Fig. 4-b where we compare the time responses calculated using the exact phase and the reconstructed phase. We clearly observe the discrepancy between the exact and reconstructed phase and exact and reconstructed time response, which is caused by the zeros lying in the right half of the s-plane and the zero and infinite complex frequency. For this case, a simple constant phase shift is not sufficient, and we need an advanced method to correct the phase error.

Blaschke product and modified Kramers-Kronig relations
At this point, we resume the theory on KK relations for magnitude and phase [3]. We address the possibility of zeros of the response function in the right half of the s-plane. The logarithm of the magnitude is not defined at such zeros. The symmetry relations, after analytic continuation for complex = + , are written as It implies that zeros of appear in pairs in the s-plane. Zeros with a vanishing real part form an exception. There may also be zeros for infinite complex frequency. In order to deal with the zeros of the response at finite frequency, we follow [3] which has been preceded by the seminal papers [13] [14]. It is specifically assumed that the complex numbers , = 1, … , with Re > 0 are zeros of . The Blaschke product is then defined as Next the modified response function, ̂( ), is constructed as The expansion of ̂( ) is possible if the integral satisfies i.e. converges. The Blaschke product should converge as well, which can be reformulated as [13] ∑ Re This condition is obviously fulfilled for finite . The magnitude of the responses ( ) and ̂( ) are equal | ( + )| = |̂( + )|.
It can be verified by first recalling that the zeros of ̂( ) in the right half s-plane appear in complex conjugated pairs. We can therefore write the product for real as with index labelling the pairs. Since we get | ( )| = 1. Hence, the validity of relation (33) is determined. The modified response function ̂( ) by applying the Blaschke product has by construction no zeros in the right half of the s-plane. Therefore, we obtain the KK relations for its magnitude and phase ̂ The phases of the responses are related by where arg ( ) can be seen as a phase correction term which corrects the phase error caused by zeros of ( ) located in the right half of the s-plane. The arg ( ) term is obtained by first calculating The complete modified KK relations are written as and ln|̂( + )| In practice, the locations of the zeros may actually not be known. In Section 5 we propose a method which enables us to find the locations of the zeros using only amplitude measurements of the time response while using specific excitation signals. This allows us to calculate the phase correction term arg ( ) in (37), needed to get the correctly reconstructed phase for the transfer function.

Additional modifications for zeros at infinity
In subsection 3.1 we discussed that transfer functions which have a zero at infinite complex frequency i.e. | | → ∞, ( ) → 0, do not satisfy the KK relations for magnitude and phase. Here we propose to introduce once more a Blaschke factor for zeros at infinity. We first assume that possible other zeros are already being taken care of by the procedure given in subsection 4.1. Next, we assume that for some positive integer . Analogous to the use of the Blaschke product, we define a modified response At this point we specify the Blaschke factor as which has a multiple zero at = −1 but none in the right half of the s-plane. Furthermore, it is easily seen that and consequently Therefore, we obtain the KK relation for the modified response ̅ For imaginary arguments of ( ) we get log| ( )| = 1 2 log( 2 + 1) and arg ( ) = arctan2( , 1) (mod 2π).
Hence, we obtain the modified KK relation When also incorporating the Blaschke product and corresponding phase correction for zeros laying in the right half of the s-plane, the modified KK relation finally becomes where arg ( ) is defined in (40).

Method for obtaining zero locations using time domain measurements
In subsection 2.3 it is noted that the cavity transfer function in some cases does not satisfy the minimum phase condition. In that case a severe phase error can emerge when the unmodified KK relations are used to reconstruct the phase response from the magnitude of the frequency response. The locations of the complex zeros within the right half of the splane need to be determined in order to correct errors in the reconstructed phase by using modified KK relations, as discussed in subsection 4.1. In this section we propose a method to estimate the locations of the complex zeros in the right half of the s-plane using specific excitation signals and measuring the resulting time responses of the system.
We define a transfer function of a system in the domain The system has input ( ) and output ( ), and represents an arbitrary linear time-invariant single-input single-output system, such as a cavity with input electric field strength of an excitation source and as output the measured electric field strength. Let us further assume that our zero-search input signal in the time domain ( ) is given by where ( ) is the Heaviside step function. With > 0, ( ) is an exponentially increasing cosine wave, which is made causal by multiplication with ( ) . The zero-search input signal ( ) is switched off at = by multiplication with (− + ), which ensures that the Laplace transform of ( ) exists for > 0. The function ( ) is defined such that the signal is zero for = 0.
The one-sided Laplace transform of ( ) is given by which has poles at = ± and = ± . Note that the Laplace transforms of ( ) and ( ) have equal poles. When the first pole pair is placed at = ± , with and corresponding to the -th zero pair of the transfer function, the k-th zero pair of the transfer function and the corresponding pole pair of the zero-search input signal cancel each other. In the time domain this can be detected by analysing the output of the system, ( ; , ), and observing if the exp( ) cos( ) term in ( ) is suppressed by the zero present in the transfer function.
We will illustrate this using a simple example. Let us define the transfer function which has one real zero at = 1, and one pole pair at = ± . Note that the zero is positioned in the right half of the s-plane on the real frequency axis, and as a result this is a nonminimum phase system. As a zero-search input signal we use ( ) from (55) with = 0 in order to find the location of the zero on the real frequency axis. After performing the inverse Laplace transform on ( ) ( ) we obtain the system response ( ; , = 0) to the input ( ) for 0 < < ( ; , = 0) = ( 2 − 1) exp( ) ( 2 + 1 ) -2 (cos( ) + sin( )) + (cos( ) − sin( )) ( 2 + 1 ) . (58) Note that for = 1 the first term in ( ; , = 0) containing the zero-search input signal exp( ) , becomes zero. This is because ( ) in (56) has a pole at = which, for = 1, cancels the zero at = 1 in ( ). Thus, we sweep until we observe the first term, containing exp( ) , to become zero in ( ; , = 0) and thereby the location of the zero is determined. The fact that exp( ) increases exponentially in time makes this relatively easy to detect in practice, e.g. by measuring the magnitude of the system response. The exp( ) cos( ) signal is detected in the output by integrating the magnitude of the system response over time as follows When the exponential function present in the output ( ; , ) is suppressed by the (complex) zero in the transfer function, we observe a decrease in ( , ) for a certain combination of and . This minimum provides us the estimated location of the pertaining complex zero.
Due to the exponential term in (55), the input ( ) will increase rapidly when increases. To ensure that the output signal ( ; , ) does not becomes too large numerically we scale the input according to where is a positive constant.

Blind tests for verification of the zero-search method
We now test the zero-search method on transfer functions having one or more zero pairs in the right half of the s-plane at varying locations. To enforce our confidence in the method, we perform several blind tests. This means that no a priori knowledge about the zero locations is used for finding them. The blind testing method used here is done through the following procedure: The second author has a transfer function of which only he knows the locations of the zeros in the right half of the s-plane. The second author then calculates the impulse response of the transfer functions and hands it to the first author. Note that the first author only has the numerically calculated impulse response of the system and no information about the corresponding transfer function, thus he does not know the locations of the zeros. Note that for computing ( ; , ) in (59) the convolution between the impulse response and the zero-search input signal defined in (55) can be used. The first author then uses the zero-search method described in Section 5 to attempt finding the locations of the zeros in the right half of the complex s-plane.
The transfer function used in these blind tests is defined as where labels the number of pairs and = + with real and positive , , . This form of the transfer function is chosen such that all zeros are in the right half of the s-plane are mirrored across the axis with respect to the locations of the poles. It is very similar to the product expansion of the causal S-matrix describing scattering of the Maxwell field [13].
We discuss the result of one of the blind tests in detail. The example contains two zero pairs in the right half of the splane where the parameter values used are 1 = 0.65, 1 = 5, 2 = 1.3, 2 = 10 and = 1. This results in the zeros of (61) being located at = 0.65 ± 5 and = 1.3 ± 10 . For the analysis we set = 10[s]. The inverse exponential scaling constant is set to = 10, which ensures that the input signal | ( )| never exceeds 1. In Fig. 5, as defined in (59) is plotted for varying and . The estimated location of the first zero is ( = 0.652, = 5.01), which is close to the actual location. Next, we observe the location of the second zero estimated at ( = 1.29, = 9.98) which is also close to the actual location. As the zeros form complex conjugated pairs, we only need to localize 2 out of 4 zeros. The accuracy of the estimated zero locations here is mainly limited by the sample density of and used to find the zeros.
More of these cases, with either one or two zero pairs were blindly tested for varying values of , , . All these tested cases provided similar accuracy in terms of the estimated zero locations. Thus, we can conclude that this method looks very promising for finding the zero locations without any a priori knowledge of the transfer function.

Fig. 5.
plotted as a function of and .

Numerical results for the cavity transfer function
As a next step, we test our method to search for the complex zeros of the cavity transfer function defined in (16). The zero-search method is first applied to the example shown in Fig. 3. The magnitude of the transfer function for the x component of the electric field is displayed in Fig. 6-a. is plotted as a function of and in Fig. 6

-b
After visual inspection, three complex zeros located in the right half of the s-plane are identified. Other zeros are located on the axis, as can be seen e.g. around ( = 0, = 2 × 10 9 ), a b Fig. 6. (a) The magnitude response of ( ) in Ω/ in the complex s- plane. (b) plotted as a function of and .
but these do not pose a problem as we use the parameter to ensure that these zeros are positioned left of the integration path in (26). Note that all zeros come in complex pairs, mirrored around the axis. We aim to find the locations of the zeros located in the right half of the s-plane, which are marked by an arrow in Fig.  6-a. To find the locations of these zeros we excite the cavity using the input function ( ) defined in (55). The output signal, ( ; , ), is then analysed by calculating the integral for ( , ) as defined in (59).
To this end, we will sweep both and . The time duration of the excitation is set to = 200 ns. The inverse exponential scaling constant is set to = 200 × 10 −9 , which ensures that | ( )| never exceeds 1. The results are shown in Fig. 6-b. We clearly observe that the locations of the zeros are marked by a decrease in ( , ). Comparing Fig. 6-a to Fig. 6-b we see that the decrease in ( , ) appears at exactly the same place where the transfer function becomes zero in the s-plane. Hence, for this case the zeros can be accurately located by varying and and observe the local minima of ( , ) where the gradient of ( , ) becomes zero. Note that because the zeros form complex conjugated pairs, the other three zeros for negative are now also automatically found. All three zero pairs, present in the right half of the splane, are determined. The accuracy of our localisation method is mainly limited by the sample density of and used to find the zeros. With the zero locations known, the corrected phase is calculated using (53) and the result in Fig.  7 is obtained.
Comparing Fig. 7 with Fig. 3-a we observe that the Blaschke product does indeed correct the phase error caused by zeros located in the right half of the s-plane and the phase error caused by the zero at infinite complex frequency. Thus, we demonstrated that by using modified KK relations, the phase can be successfully determined from the magnitude of the frequency response.
It is verified that the phase error caused by the zero at infinite complex frequency is negligible within the bandwidth considered here, and thus (53) and (41) provide similar numerical results. This is likely due to the large number of higher order modes, which causes the magnitude of the complex transfer function to go to zero slowly for increasing frequencies. Other cases however might benefit from using (53) instead of (41) Fig. 7. Magnitude of ( ) and the reconstructed phase of ( ) corrected using the Blaschke product (modified KK relations).
when the magnitude of the complex transfer function goes to zero more rapidly for increasing frequencies.

Noise sensitivity
In this section we investigate the influences of noise on the accuracy of the proposed zero-search method and the phase reconstruction method using the modified KK relations. This will be illustrated using the cavity transfer function defined in (16) for the example shown in Fig. 4.
To apply the proposed phase reconstruction method using the modified KK relations in combination with the zero-search method, two separate measurements must be performed. The first measurement is performed in the time domain, where the time responses of the system are measured to find the zeros in the right half of the s-plane. The second measurement is performed in the frequency domain where the magnitude of the frequency response of the system is measured, which is used to reconstruct the phase using the modified KK relations. In practice, both measurements will be subject to measurement noise.
First, we will simulate the effects of measurement noise on the zero-search method. Gaussian noise with zero mean and a standard deviation of 10 −4 V/m is added to the measured time response, ( ; , ), in (5). The magnitude of the transfer function for the x component of the electric field is displayed in Fig. 8-a.
is plotted as a function of and in Fig. 8-b. a b Fig. 8. (a) The magnitude response of ( ) in Ω/ in the complex s- plane. (b) plotted as a function of and .
We observe that all zeros except the zero at ( = 0.1 × 10 6 , = 7 × 10 8 ) can be located with high accuracy. The zero that cannot be located using our zero-search method is positioned very close to the poles on the imaginary frequency axis. This causes the influence of this zero to be nearly cancelled by a nearby pole, making locating the zero using the measured time response very difficult. Note that we consider a lossless cavity here. In practice the poles will usually be positioned to the left of the imaginary frequency axis due to electromagnetic losses in the cavity. Our hypothesis is that this should make locating zeros close to the imaginary frequency axis easier in practice than encountered in this exercise.
In the example shown here, the noise only influences the accuracy of locating the zeros positioned far away from the imaginary frequency axis, as the measured signal strength is very low around these zeros. The average signal to noise ratio over time (RMS) of the measured time response is approximately 2 when = 7.7 × 10 7 and = 3.3 × 10 8 , which is just sufficient to be able to locate the corresponding zero. When the noise floor is too high, zeros cannot be located as the they disappear below the noise floor. This problem can be circumvented by applying a stronger excitation signal to the system, e.g. by decreasing in (60).
Next, we simulate the effects of measurement noise on the phase reconstruction method using the modified KK relations. Complex Gaussian noise with zero mean and a standard deviation of 3.16 Ω/m for both real and imaginary part is added to the frequency response, ( ′ + ), in (41).
With all but one of the zero locations known, the corrected phase is calculated using (53) and the result in Fig. 9 is obtained. Note that all noise outside the frequency band of interest, 50 MHz---125 MHz, is filtered out using a bandpass filter.
Comparing the result in Fig. 9 to the result in Fig. 4-a we see that, despite the fact that noise has been added here, the phase error is significantly reduced by using modified KK relations. We see 5 distinct locations where the absolute phase error becomes larger than 20°. The first 4 spikes in the absolute error are caused by a low signal to noise ratio. However, at these locations the transfer function is very small and thus these errors will not have a significant impact on the reconstructed time response. Fig. 9. Magnitude of ( ) and the reconstructed phase of ( ) corrected using the Blaschke product (modified KK relations). Also, the absolute error of the reconstructed phase is displayed.
The last error spike around 110 MHz is caused by one of the missing zero's in the correction factor, which we could not find using the proposed zero-search method. The effect on the reconstructed time response depends on the bandwidth of the input signal as only frequencies above 110 MHz are affected.
In this section we aimed to apply a realistic amount of measurement noise to determine the noise sensitivity of our zero-search and phase reconstruction method. In future research we need to further investigate the influences of noise on our methods and draw up requirements for measuring equipment used for performing experiments.
The localisation of unknown zeros of a transfer function using time domain responses, and the subsequent reconstruction of the phase response may be interpreted as an inverse problem. The locations of the zeros as well as the number of zeros in the right half of the s-plane are, in our approach, parameters which need to be estimated appropriately. Note that we search for local minima in , indicating the zero locations, and not for a global minimum regarding an inverse problem. In practice however, false minima might occur, e.g. due to measurement noise or errors. As a result, the system might seem to have zeros that do not exist. More research is envisioned on this topic. In the examples discussed here we did not observe this problem. To solve an inverse problem, normally one would define a cost function or an associated error norm which needs to be minimized. Here, we did not yet define an error norm to localize the zeros but use a visual approach instead.

Conclusions and recommendations
The electromagnetic field inside a PEC cavity is efficiently modelled using an eigenmode expansion when the EM fields can be sufficiently well approximated using a reasonable number of modes. With this eigenmode expansion we determine a transfer function from an excitation point dipole source to an observation point of the electric field. The locations of the zeros of this transfer function, within the complex s-plane, depend on the shape of the cavity and the locations of the observation point and excitation source.
If only the magnitude of the transfer function is known (by measurements), the KK relations may be used to reconstruct the phase from the magnitude. This method provides correct results when the minimum phase condition is satisfied, which is equivalent to the requirement that the transfer function has no zeros in the right half of the s-plane. We have observed that the transfer function of a cavity, as described here, may exhibit several zeros in the right half of the s-plane. As a result, the minimum phase condition is violated. This introduces a substantial error in the reconstructed phase when using the KK relations, which in turn introduces an error in the reconstructed time response of the electric field. However, we have observed that for some cases the overall shape and most importantly the envelope and peak values of the reconstructed time signal are preserved compared to the exact analytical calculations, despite the minimum phase condition being violated.
In other cases, the phase error present in the reconstructed phase is not acceptable. In order to handle these cases we used the Blaschke product, from which modified KK relations are derived. The Blaschke product allows for calculating a phase correction term which corrects the errors introduced due to the violation of the minimum phase condition. This modification requires the locations of the zeros in the right half of the s-plane to be accurately known. If the locations of these zeros are known, the Blaschke product and corresponding phase correction can be calculated. Subsequently, the reconstruction of the phase using modified KK relations is performed appropriately. In our numerical studies, however, within the limited bandwidth considered here, the influence of zeros at infinity turned out to be negligible in the case of the cavity transfer function.
Generally, the locations of the zeros of the transfer function in the right half of the s-plane are unknown, if only the magnitude of the frequency response is measured. However, it is possible to locate these zeros by analysing a specific time response of the system. When we choose the zero-search input signal such that it has a pole that coincides with a zero of the transfer function, a (significant) decrease in energy in the measured time response can be observed. This is due to the fact that a pole of the zero-search input signal gets cancelled by a zero of the transfer function. To this end an exponentially increasing cosine zero-search input signal is used as defined in (55). By varying the exponential growth constant and the oscillation frequency of the zero-search input signal, the zero locations in the right half s-plane of the transfer function are obtained. The usability of this approach was tested by using a blind testing method where only the time response of the system is given to the first author. These blind tests provided positive results for all verified cases, as all zeros were located successfully without a priori knowledge about the transfer function.
Additionally, the method is applied to the cavity transfer function in two examples. For the first example, all of the zeros in the right half of the s-plane are located, and the phase is accurately calculated from the magnitude of the frequency response by using modified KK relations. For the second example, one of the zeros could not be located because it is positioned too close to poles on the imaginary frequency axis. Our hypothesis is that in real world measurements, locating zeros close to the imaginary frequency axis is easier since losses are present and therefore poles are positioned to the left of the imaginary frequency axis.
Measurement noise will have an impact on the accuracy of the zero-search and phase reconstruction method. However, under realistic noise conditions the methods still give good results, as was demonstrated for the cavity transfer function example here.
It is important to indicate that this paper demonstrates a proof of concept for applying modified KK relations to nonminimum phase systems. More research is currently being invested by the authors in the practical realisation of the method and measurements have been planned for this purpose. This will also incorporate more in-depth research on the influences of measurement noise.