Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.cplire.ru/html/lab234/pubs/2003CompPhysCp171.pdf Äàòà èçìåíåíèÿ: Mon Feb 28 21:09:04 2005 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 10:46:39 2012 Êîäèðîâêà: IBM-866 Ïîèñêîâûå ñëîâà: ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï

Computer Physics Communications 151 (2003) 171í186 www.elsevier.com/locate/cpc

Deconvolution problems and superresolution in Hilbert-transform spectroscopy based on a.c. Josephson effect
E.L. Kosarev
a,

, A.Ya. Shul'man b , M.A. Tarasov b , T. Lindstr?m

c

a P.L. Kapitza Institute for Physical Problems, Russian Academy of Sciences, ul. Kosygina, 2, Moscow 117334, Russia b Institute of Radioengineering and Electronics of the RAS, 101999 Moscow, Russia c MINA, Department of Physics, Chalmers University of Technology, S-41296 G?teborg, Sweden

Received 31 May 2002

Abstract The analysis of the numerical aspects of Hilbert transform spectroscopy based on the a.c. Josephson effect is presented. The resolving power of Hilbert transform spectroscopy is determined by such factors as the linewidth of the Josephson oscillations (intrinsic or natural resolution limit) and the limitation of the measurement interval (extrinsic or technical resolution limit) like in any spectroscopic technique based on some integral transformations. The deconvolution problem in Hilbert transform spectroscopy is posed and its solution is considered using the approach of the 1st kind integral equation for the spectrum of the incident radiation constructed from the input data of the Hilbert transform spectroscopy--the `hilbertogram'. The program package RECOVERY based on the maximum likelihood method is used for this purpose. This method allows to attain the maximum possible resolution enhancement in output result for a given signal-to-noise ratio in the input experimental data. The samples of numerical simulations and the spectrum of frequency-modulated BWO radiation measured by means of the Josephson junction made from high-Tc superconductor are presented. It is shown also that the integral equation approach allows to recover the sought spectrum beyond the intrinsic resolution limit and to achieve the superresolution. 2002 Elsevier Science B.V. All rights reserved.
PACS: 85.25.Cp; 85.25.Pb; 95.85.Bh; 29.30.-h; 29.40 Wk Keywords: Hilbert transform spectroscopy; Josephson effect; Deconvolution; RECOVERY; Superresolution

1. Introduction It was shown in [1] that the change in (I -V ) characteristic of a Josephson junction irradiated by a wide-band incoherent electromagnetic field is related to the spectrum of the incident radiation by the Hilbert transformation. Based on this relation a spectroscopic technique for submillimeter wavelength range was suggested and later named as Hilbert (-transform) spectroscopy [2]. Since 1980, the mainstream of publications in this filed was connected
Preliminary version of this paper was reported to the 2nd International Conference "Modern Trends in Computational Physics", Dubna, Russia, July 24í29, 2000 * Corresponding author. E-mail addresses: kosarev@kapitza.ras.ru (E.L. Kosarev), ash@cplire.ru (A.Ya. Shul'man).

0010-4655/02/$ í see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S0010-4655(02)0 0 7 0 1 - 4


172

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

with the choice of junction type suitable for the Hilbert-transform spectroscopy implementation [3í5]. Recently, the Hilbert transform technique was implemented to the response of high-Tc superconductor junction at the 1-THz frequency [6]. However, not much attention has been paid to the numerical and pure spectroscopic aspects of the Hilbert spectroscopy technique. The aim of the present work is to fill this gap. The deconvolution problem in Hilbert spectroscopy is considered using the approach to the Hilbert spectrum (hilbertogram) as the first kind integral equation relative to the spectrum of the incident radiation. We use for this purpose the package RECOVERY, which is based on the maximum likelihood method. This method and the program were chosen because the method attains the maximum possible resolution enhancement for a given signalto-noise ratio in the input experimental data [7,8]. Just a minor modification has been necessary to take account of the singular character of the Hilbert transformation kernel. The solutions of the integral equation have been obtained in the course of numerical simulations and compared with the results of the usual Hilbert transform procedure using the FFT algorithm. It is also shown that the integral equation approach allows the recovery of the sought spectrum beyond the intrinsic resolution in the case of the experimentally measured response of a Josephson junction made from a high-Tc superconductor and to realize the superresolution in accordance with the definition and criterion stated in [7]. It is worthwhile to define more exactly the terms deconvolution and superresolution as applied to the present case. Because the action of the point-spread (or instrumental) function (PSF) of a spectroscopic device is usually described by the convolution integral (see e.g. [9]) the deconvolution procedure should remove distortions in the measured spectra introduced by the PSF. This can be achieved by various techniques but we keep the term deconvolution for the restoration of altered but measured details in the spectrum under study. The definition of these details may be given, e.g., by the known Rayleigh's resolution criterion for the original spectrum or by the effective bandwidth of its Fourier transformation. The term superresolution can be defined as a procedure for recovering the spectral details that have been projected to zero set or suppressed below the noise level owing to the convolution with the PSF [7,12]. In the case of the Hilbert-transform spectroscopy we define as deconvolution any way that allows the minimization of spectrum distortions originating from the limitation of the measurement interval. The superresolution is defined as the procedure for removing to some extent the effect of the non-zero bandwidth of the Josephson oscillations. The results of the present work show that the integral equation approach to spectrum recovery in Hilbert spectroscopy provides simultaneous solutions to both deconvolution and superresolution problems.

2. Basics of Hilbert-transform spectroscopy 2.1. Josephson effect in brief The readers which are interested in a more detailed description of this topic can refer to the textbooks [13í15]. The superconducting state is characterized by the complex order parameter (condensate wave function) = exp(i),

where is related to the energy gap for normal quasiparticle excitations and is the phase of the order parameter. The Josephson effect arises as two superconductors are separated by a thin region (barrier) where is suppressed, but a tunneling process for Cooper pairs is still possible. In this case the current I of Cooper pairs (or superconducting current) across the barrier is expressed by the formula I = Ic sin , = L - R , (1)

where is the difference of the phases of left and right superconductors and Ic is named as the critical current of the d.c. Josephson effect. While I < Ic, that is the d.c. current provided by the external circuit is sufficiently


E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

173

low there is no voltage drop across the junction and the charge transfer through the junction is due to the pure superconducting current (supercurrent). Besides the superconducting current that accompanies the spatial inhomogeneity of the phase, there are situations when the phase difference is time-dependent. In this case there is the voltage drop V across the barrier and the famous Josephson relations apply: d = 2eV , dt V = 2eV /h. ï h ï (2) (3)

The immediate consequence of the two Eqs. (1), (2) is that if there is a constant-in-time V then the superconducting current oscillates with the Josephson radian frequency V given by Eq. (3). These `self-sustaining' oscillations are the basis of all spectroscopic applications of the Josephson effect. Interacting non-linearly with the external time-dependent perturbations they produce changes in d.c. currentvoltage (I -V ) characteristic of the Josephson junction or signals with some intermediate frequencies. The simplest model of the Josephson junction is the Resistively Shunted Junction (RSJ) model. The total current across the junction is considered as the sum of the superconducting current and the current of normal quasiparticle excitations (electrons): V/R + Ic sin( ) = I, where R is the junction resistance in the normal (non-superconducting) state. This equation and the Josephson relations Eq. (2) constitute the set of equations describing all the properties of a small-area junction. It is generally accepted to use the dimensionless units. Let us denote i = I/Ic ; Vc = Ic R ; v = V/ Vc ; = c t. c = 2eVc /h; ï = /c ;

Then we can rewrite the set of the equations in the dimensionless form: d d + sin = i, = v. (4) d d If we could provide the time-independent bias voltage to the junction (the voltage-driven junction) then the d.c. I -V characteristics (the time-averaged current) would be a simple linear relation i = v because the supercurrent ï is the pure harmonic function of time in this case. The solution of these equations in the more realistic case of time-independent current i (the current-driven junction) gives a non-linear I -V curve that can be written in the form: i = sign(v) v 2 + 1, ïï (5) where v is the time-averaged bias voltage and the bar indicates a time-averaging procedure. ï It must be stressed that the voltage in the junction is in this case a periodic (but not harmonic) function of time. Its fundamental frequency is determined by the Josephson relation (3) taken for time-averaged bias voltage, that is ï v = v . ï Let now the monochromatic external current i( ) with radian frequency be supplied to the junction. Then ~ the d.c. I -V characteristic is changed. The respective expression describing the current response under condition v = const was obtained by Kanter and Vernon [16] in the quadratic approximation: ï ~ 1 i 2 () . (6) 2 - 2 4i v ïï Evidently, the above-written expression is invalid near the `resonance' bias where the condition v = is fulfilled. ï More detailed analysis shows that the `current step' (Shapiro step) is formed at v = . The magnitude of this step ï i(v, ) = ïï


174

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

is proportional to the first grade of the amplitude i(). Inside this interval of the d.c. current the frequency of the ~ Josephson oscillations becomes locked to the external frequency (see e.g. [13]). Let us consider now the KanteríVernon formula (6) and ask whether it can generate the Hilbert transformation. At first glance, it looks like the kernel of the Hilbert transformation. If we had the right to describe the junction response to the a.c. current with a continuous spectrum Si () in the form


i(v) - Si () i(v, ) d, ïï ïï
0

where - denotes the Cauchy principal value of the integral, then we might hope to obtain the spectrum S() of the incident radiation by means of the inverse Hilbert transformation. However, near the resonance voltage v = the ï response is not quadratic but linear in the a.c. current amplitude. Hence, this situation has to be considered more thoroughly. 2.2. Theoretical background of Hilbert-transform spectroscopy The KanteríVernon expression (6) was obtained by neglecting the fluctuation current if ( ) in the junction. In fact, such fluctuations are always present and must be inserted into the expression for the current. Thus, we have to rewrite the first equation of our previous set (4): d + sin d The solution of spectrum of the the response to =ï + if . i (7)

(7) shows that the presence of the fluctuations gives rise to the finite width of all harmonics in the Josephson oscillations and eliminates the singularity in Eq. (6) [13]. The analytical expression for the monochromatic perturbation has the form [1]: v+ ï i 2 () ~ v- ï + 8i v (v + )2 + 2 (v - )2 + ïï ï ï
2

i(v, ) = ïï

.

(8)

It was derived in the RSJ model using the results of [13] under the condition that the fluctuation current if is -correlated in time. Here is the fluctuation-induced linewidth of the Josephson oscillations. This expression is valid for v, . Now we are able to establish that the a.c. current with a spectral intensity distribution Si () ï should result in the response determined by the formula [2]:


i(v) = ïï
0

1 d i(v, ) = ïï - 8i v ïï



dSi ()
-

-v ï ( - v)2 + ï

2

.

(9)

Some comments are necessary with regard to the relation considered above and the spectrum of the electromagnetic coupling of the junction to the radiation is usually realized spectrum Si () has to be related with the spectrum S() of th Si () = K() S(),
2

~ between the spectrum Si () of the a.c. current i radiation incident on the Josephson junction. The by means of various kinds of antennas. Thus, the e incident radiation by the expression

where K() is the transfer function of the antenna. Since we are interested in the radiation spectrum S(), the function K() must be taken into account for the case of sufficiently wide-band external radiation. It is seen that the expression in brackets in Eq. (9) tends to the Hilbert transformation of the Si () if the limit 0 is taken. Entering the function g(v) = ï 8 i v i(v) ïï ï ï (10)


E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

175

and applying the Hilbert transformation to it in the form: S() = Hv ï we obtain [2] 1 S() =




1 g(v) ï g(v) - dv ï ï v- ï
-



(11)

d Si ( )
-

. ( - )2 + 2

(12)

The estimation of the spectrum Si () is noted in Eqs. (11), (12) with S(). It is easy to recognize in Eq. (12) the convolution of the spectral density of the external a.c. current with the Lorentz distribution. The latter represents the sideband J( - v, v) of the known spectrum of the intrinsic ïï Josephson oscillations in the current-driven junction SJ ( , v) ï (v) ï (v) ï + J( - v, v) + J( + v, v) ïï ïï ( - v)2 + 2 (v) ( + v)2 + 2 (v) ï ï ï ï (13)

obtained theoretically in the RSJ model with -correlated thermal noise. If Si () varies slowly over frequency intervals of the order of then the Lorentz distribution in the integrand may be replaced by ( - ) and the function obtained S() becomes evidently equal to the sought spectrum Si (). On the other hand, if Si () is a more narrow band function of the frequency than the spectrum of the Josephson oscillations and can be taken as ( - 0 ) then the last expression gives us the intensity of the Josephson oscillations J(0 - v, v) at the ïï ï frequency 0 as a function of the bias v . It should be noted that the bandwidth of the Josephson oscillations is really bias-dependent as it is explicitly shown in Eq. (13) by the second argument of J( ‘ v, v). In the scope of the same RSJ model the fluctuation ïï current if originates from the equilibrium thermal fluctuations in the normal resistance R of the junction and the result is (in dimensional form) [13]: (V) (v)c = ï 1 2e h ï
2 2 Rd (V) I2 kT 1 + 2c . ï R 2I (V)

ï Here Rd = d V/d I is the differential resistance of the junction. It is seen from this expression that depends in general on the bias voltage V . This circumstance makes the transfer from Eq. (9) to Eq. (12) by means of the Hilbert transformation (11) questionable. Below we point out the conditions when this V -dependence of should not be an obstacle in the Hilbert-spectroscopy implementation. It is necessary to note that the real Josephson junctions can more or less differ from the ideal RSJ model in their characteristics. On the one hand, the RSJ model does give a good description of some Josephson junctions. The discussion of the RSJ model and its applicability for low-Tc superconductor junctions can be found in [14]. As for high-Tc -superconductor Josephson junctions, there are the successful experimental tests of the RSJ model with reference to high-frequency applications and noise properties in [4,5,17,18] and in [19] for SQUIDs. On the other hand, it is certainly of importance to know to what extent the results obtained by the Hilbert spectroscopy technique are independent of the junction's "non-ideality". This problem is now the subject of current investigations, especially as regards high-Tc based junctions or more sophisticated Josephson structures like coherent chains of junctions. In particular, the suggestion to use Eq. (12) for measuring the spectrum of Josephson oscillations was formulated in [2] as a preliminary proposal but it was already applied to a quantitative analysis in [4,5] without any substantiation.


176

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

We deduce here, following [20], an extension of Eq. (9) that allows us to remove the limitation of the Josephson oscillation spectrum by the Lorentzian shape. It is easy to see that Eq. (9) with account for Eq. (13) can be presented in the form: 1 i(v) = ïï 8i v ïï


-

1 J( - v, v) ïï - d dSi () -
-



.

(14)

Eq. (14) can be transformed to the form: i(v) = ïï 1 8i v ïï dSi () -
-

J( - v , v) ï 1 - dv v -v ï
-



.

(15)

Under condition (v)/v 1 over all the essential range of the bias v , the second argument of the spectral density ïï ï J in Eq. (15) does not hinder us from calculating the Hilbert transformation according to Eq. (11) at arbitrary but smooth dependence on v . As a result, we arrive at the generalization of the convolution in Eq. (12) in the ï following manner 1 S(v ) = ï


dSi ()J ( - v , v). ïï
-

(16)

It follows from (16) that Hilbert-transform spectroscopy can be used for measuring either the spectrum of the incident radiation or the spectrum of the Josephson oscillations itself independently of the real spectral distribution of the Josephson oscillations if one of the involved spectral distributions is much narrower the other. Eq. (16) constitutes the basis and establishes the restrictions of experimental investigations where Hilbert spectroscopy is employed to measure the shape of the Josephson oscillation spectra and the dependence of on V . If we are interested in measuring the spectra of external radiation sources with necessary account for the finite bandwidth of the Josephson oscillations, the expressions (15), (16) allow us to deal with the deconvolution problem in functional spaces of either measured or Hilbert transformed functions. 2.3. What is to be measured and how The function g(v) of Eq. (10) must be formed as the result of measurement procedures. To this end the I -V ï ï ï characteristic I(V) and the response function I(V) are to be measured as a function of the bias voltage V . The range of the bias V -sweep has to include the [0íVmax] interval with sufficiently high Vmax in order to provide effectively the infinite limits during an evaluation procedure of the integration over V in Eq. (11). It is convenient to introduce the term hilbertogram for the function g(v) on close analogy with the interferogram in the Fourierï transform spectroscopy. The resolution of Hilbert spectroscopy depends on the linewidth of the Josephson oscillations, on the length of the swept V -interval, and on the spacing of the data points. The linewidth depends on the intrinsic voltage fluctuations in the junction and on the level of external interferences. The latter must be suppressed as low as possible in the wide frequency range beginning from the line frequency and up to TV-station frequencies. The experimental equipment must allow to chop the incident radiation and to measure the phase of the response Iï using a lock-in amplifier. The function g(v) defined by Eq. (10) can be obtained from the measured data by ï straightforward procedure [21]. The simplest case takes place for the d.c. voltage-driven junction when two curves, the d.c. current Iï and the response Iï, must be measured as functions of bias V . For the current-driven junction it ï ï is necessary to calculate the current response I(V) from the measured voltage response V(I) by means of the known expression ï I =- V/Rd . (17)


E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

177

ï Fig. 1. The example of the experimental data as a function of the junction bias in ²V. (1) The I(V) characteristic, current in ²A, voltage in ²V. (2) The differential resistance of the junction Rd in . (3) The value of V(V) in a.u. (4) The hilbertogram function g(V) in a.u.

ï The differential resistance Rd (V) can be directly measured or numerically calculated from the measured I(V) characteristic. ï In Fig. 1, for an example, all constituents of the source data set I(V), Rd (V) and V(V) are presented. The hilbertogram g(V) is generated from these data by using Eqs. (10) and (17). This sample was measured with a high-Tc Josephson junction irradiated by the backward-wave oscillator (BWO) radiation at the frequency of about 410 GHz. 2.4. Transfer to the integral equation with the Hilbert kernel To calculate the Cauchy principal value in the Hilbert transformation 1 f(x ) H{f }(y ) = - dy x-y
-

(18)

the two numerical methods are usually applied: the direct computation by means of a quadrature formula or using the fast Fourier-transform (FFT) technique [23í25]. Both require quite a large region of the bias voltage to be swept in order to avoid the lack of resolution and/or distortion of the measured spectrum. The method for computing the Hilbert transformation (18) proposed by Weideman [26] is reduced to the change of variable x = tan , 2 - < < ,

and it also requires knowledge of the integrand f(x ) over an essential part of its definition. In this study we use another approach. Applying the inverse Hilbert transformation, we can invert Eq. (11) and rewrite it as a convolution integral equation for S() with the Hilbert kernel
-1 Hï v

S() 1 d = g(v). ï S() - - -v ï
-



(19)


178

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

The measured radiation spectrum S() and its estimation S() belong usually to the class of functions with compact support and finite energy, i.e.


S() = 0

and S() = 0 for || >1 ,1 > 0,
-

S() d< .

The Hilbert transformations g(v) of such functions have, on the contrary, long-range tails and slowly decrease as ï 1/v as v tends to infinity. The application of the integral equation approach allows us to avoid measuring g(v) in ï ï ï the region of v, far beyond the existence domain of S(). It will be shown in Section 5 that the measurement of ï g(v) can be restricted by a region v where the magnitude of g(v) still exceeds the noise level considerably. ï ï ï To solve the convolution integral equation (19) with the Hilbert kernel we use the program package RECOVERY [8] which is based on the maximum likelihood method. It was shown in [7] that the maximum likelihood based algorithm allows us to reach superresolution in the case of non-singular kernel of the integral equation. The present work demonstrates that it is also correct for the singular kernel of the Hilbert transformation.

3. Definition of superresolution for the Hilbert-transform spectroscopy There is a natural connection of Eq. (19) with the KramersíKronig relations in electrodynamics and optics [22] and there are computer programs and algorithms specially devoted to compute them [23í26]. All methods implemented in these papers can be outlined as linear methods because all of them use a linear expansion of the data on different sets of basis functions. According to paper [7], any kind of linear restoration methods can only modify the amplitudes of the Fourier harmonics but cannot generate the new ones that are absent in the input data or were lost in the input noise. Therefore, linear methods can not achieve a superresolution. From here onwards we shall define the superresolution ratio SR as the ratio of the Fourier spectrum width of the output signal to the input one SR = ( )out /( )in . The definition of for any signal is given by the formula


(20)

1 = Smax ()

S() d,
0

(21)

where S() is the spectral intensity of the signal with respect to the positive frequency semiaxis and Smax () is the maximum value of the function S() on the interval 0 < < . It must be emphasized that the definition (20) of SR does not formally coincide with the definition given in [7] because the kernel of the Hilbert transformation K(x - y) = 1 x-y (22)

has infinity bandwidth both in the direct x -space and in the Fourier-transform -space. But the new definition (20) can be used with equal ease for the Hilbert transform kernel and for ordinary kernels like Lorentzian or Gaussian ones. It was shown in [7] that the maximum possible superresolution ratio is mainly determined by the signal-to-noise ratio in the experimental data and it can be easy predicted by the simple formula SRmax = Es 1 log2 1 + 3 En dB . 10

In the above formula Es --energy of the measured signal, En --energy of noise in the input data, and dB = 10 lg(Es /En ) is the signal-to-noise ratio of the input data.


E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

179

In view of well-known formulae of the Hilbert transformation H sin(a x ) = cos(ay ), H cos(a x ) =- sin(ay ), a>0

the Fourier harmonics of any function f(x ) do not change their magnitudes after the Hilbert transformation. It follows from here that ( )
out

= ( )in,

and according to our definition of superresolution (20) we obtain always the superresolution ratio SR = 1 for the Hilbert transformation. It is a very important general conclusion and at first glance we can not in principle obtain SR > 1 in Hilbert-transform spectroscopy. However it would be the wrong conclusion. In reality, we always have the finite bandwidth of the Josephson oscillations J( - v, v) in Eqs. (14)í(16). Thus, beside the Hilbert transformation in (11) we have the convolution ïï integrals of the input spectrum Si () either with the Hilbert transform spectrum of Josephson oscillations in Eq. (14) or with J( - v, v) itself in Eq. (16). As a result, the spectral bandwidth of the real kernel in the integral ïï equations for Si () is always finite and the conception of superresolution in Hilbert-transform spectroscopy is thereby well-defined. This gives, in principle, a chance to increase the bandwidth of the output result, using the recovering procedure. In contrast to the methods designed for direct computation of the Hilbert transformation of the input experimental data we suggest here to solve the convolution integral equation (19). We use for this purpose program Dconv2 from the package RECOVERY [8] based on the maximum likelihood method (MLM). This is a nonlinear method and as it was shown in [7] this method can in fact achieve superresolution. The detailed description of the recovery procedure has been presented in [7,8] and some recent examples of the RECOVERY program package applications are published in [10,11]. There is a frequency-domain analysis of the superresolution resulting from the expectation-maximization algorithm [12], which is in essence the same as used in the RECOVERY package, which is in its turn an improvement of Tarasko's algorithm [27]. The interesting paper of Kuz'menko et al. [28] should be noted here, devoted to infrared spectroscopy where the KramersíKronig relations are also reduced to a system of integral equations of the first kind solved by the nonlinear method of minimization.

4. Point spread functions for the Hilbert transform kernels and start point for the iterations in the RECOVERY programs Because the Hilbert transform kernel (22) has odd symmetry, it is convenient to use input data having also odd symmetry I (-V) =- I (V ). The kernel function (22) has the nonintegrable singularity at x = y and we can not directly use this function as a point spread function required by the deconvolution program Dconv2. Instead of (22) we accept as the PSF PSF1 (xi ,x0 ) = 1 , x i - x 0 + 0 .5 xi ,x0 [-NK ,NK - 1]. (23)

Here 2NK is the total number of points in the digitized PSF. The PSF defined in Eq. (23) has no singularity at the point x0 , the position of symmetry point of the PSF. The substitution of the PSF in form (23) is implemented in the new Dconv2_n recovering program. We do also use in our recovering program another version of the PSF y PSF2 (x ; D) = , y = (x - Npsf /2)/D , (24) 1 + y2


180

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

where D is a free parameter for achieving superresolution. The PSF2 is proportional to the inverse Hilbert transform of the Lorentzian peak L(x ) = 1 x2 + D
2

(25)

and the PSF2 has also the same asymptotic at y with respect to (22). It was implemented in both standard programs Dconv and Dconv2 of the RECOVERY package that the maximum of the recovered spectrum was placed at the point of the maximum of the input spectrum. This was changed in the Hilbert transform recovering program Dconv2_n where the maximum of the recovered spectrum is placed at the point of the sign change of the input hilbertogram data. One more point is essential in using the Dconv2 program in Hilbert-transform spectroscopy and it should be mentioned here. It was implemented in Dconv2 of the standard RECOVERY package [8] where iterations always start from the standard initial approximation S0 () = 1. We have found that it is more effective if the initial approximation S0 () = H g(V) , V = /(2f0 ), f0 = 0.4835 GHz/²V

is equal to the linear Hilbert transform of the input data g(V). This is implemented as an option in the Dconv2_n recovering program.

5. Simulation results Before demonstrating the application of the Dconv2_n program to real input data in Hilbert spectroscopy we have to test this method on simulated data. As a rule we can only use an finite interval of the input hilbertogram support 1 2 , 1 < ,2 < (usually 1 = 0) rather than - < < required by the definition of the Hilbert transformation. We call this as data limitation. The effect of the data limitation on the output result for two different techniques of the Hilbert transformation is presented in the Fig. 2. One case is the FFT based linear method combined with the optimal filtering [29] and another case is the nonlinear method based on the Dconv2_n program. The input hilbertogram having the signal-to-noise ratio S/N = 20 dB and the initial spectrum to be determined by Hilbert spectroscopy are shown in Fig. 2(A) of this picture. The results of the Hilbert transformation obtained by the two above-mentioned techniques are shown in Fig. 2(B) of this picture. We see the drastic difference between results of both methods. The implementation of the Hilbert transformation of the hilbertogram in Eq. (11) by FFT method has a large bias near = 2 , whereas the solution for the integral Eq. (19) by Dconv2_n program does not give bias and the level of the output noise equals approximately to the level of the input noise. The wrong "tail" of the "by FFT" output curve is a trace of the logarithmic divergence of the Hilbert transformation at the place of the first kind discontinuity resulting from the periodical continuation of the hilbertogram due to the FFT implementation. A striking example of superresolution in Hilbert spectroscopy is presented in Fig. 3. Here the sought-for spectrum consists of three narrow peaks having equal amplitudes and non-uniformly distributed over the frequency axis. If we do not use the information about the real profile of the Josephson generation line (curve 1) we can only obtain as the recovery result the curve 2 from part (A) of this figure. The curves 2 and 3 have the same frequency band width and according to our definition (20) the superresolution ratio should be equal to unity SR = 1. At the same time, if we do use the approximate estimation of the Lorentzian profile width D = 45 instead of the true value D = 50 we can obtain the superresolution restoration result (curve 2 of the part (B)). The superresolution ratio in this example equals SR = 1.872. We have used the PSF2 from Eq. (24) in the RECOVERY program Dconv2_n.


E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

181

Fig. 2. Data limitation effect. (A) 1--input Lorentzian peak (25) with D = 60, 2--the inverse Hilbert transform of the Lorentzian peak with the normal distributed noise having signal-to-noise ratio S/N = 20 dB added to this curve. This curve is truncated on right. (B) 1--restoration of the input Lorentzian peak by the nonlinear method using the Dconv2_n program and the PSF 1 from Eq. (23) (thin noise curve) and the original peak inside (thick curve), 2--restoration by the linear FFT based method with optimal filtering.

Fig. 3. Superresolution in Hilbert spectroscopy. (A) 1--the Josephson oscillation line Lorentzian peak 1 with the 3-line spectrum (see bottom of (B), 3--input data for the this curve is also truncated on the right side. (B) 1--the same curve as 2 from (A), (Gaussian-like peaks), the three narrow peaks inside is the sought-for spectrum. The iteration number N = 300.

with the Lorentzian profile D = 50, 2--convolution of the RECOVERY program, S/N = 30 dB. Similarly to Fig. 2 2--the RECOVERY result of 3-line spectrum restoration superresolution ratio achieved is SR = 1.872, 2 = 0.98,


182

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

The last simulation example presented here is the frequency modulation (FM) test. This case gives us an opportunity to simulate a signal with a continuous spectrum of quite a complicated character. The time-dependent function describing the oscillation with harmonically modulated frequency can be written as sin t . (26) The Fourier series of this almost periodic function may be found in many textbooks (see, e.g., [30, p. 39]). However, in our case it is more convenient to deal with the covariance function of f(t ) and its spectral density. Using the definition of the covariance function ( ) and its spectral decomposition that is appropriate for periodic and for almost periodic functions
T/2

f(t ) = cos 0 t +

1 ( ) = lim T T

dt f (t )f (t + ),
-T/2 T/2

2 S () = lim T T
0

d ( ) cos( ), S (k )e
-ik

(27)

( ) =
k (-,)

,

we obtain for the function (26): 1 2 ( ) = J0 sin cos(0 t), 2 2

12 S (k ) = Jk [ 4

mk

+ nk ].

(28)

Here Jk (x )--Bessel function with integer index k, 0 , > 0, - < < , k = 0, 1, 2, ... , the set of {k } is determined as the solutions of the equations | ‘ 0 | = k , integer parameters m and n are defined by the expressions |k + 0 | |k - 0 | , n= , m= and mn is the Kronecker characteristic function
mn

=

1, m = n, 0, m = n. / 1.


We shall be interested in the case of very low modulation frequency, that is 0 / (29) Thus, the summation over k in Eq. (27) can be approximated by the integration over in the following way ( ) = 2
k >0

S (k ) cos k 2
0

d J2 2 2 ²

2 + J



cos( ).

(30)

The obtained expression (30) defines the continuous approximation of the frequency modulated (FM) oscillation spectral density, viz. S () = J2 2 ²
2 + J



,

²=

| - 0 | | + 0 | , = .

(31)

Under conditions (29) this formula can be yet more simplified using Langer's asymptotic representation for Bessel functions [31]. Taking the limit 0 in (31) and neglecting the exponentially small second term in the brackets, we obtain


E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

183

S () =



1 2

2

- ( - 0 )
(1- (

2 -1/2
2

,
]

| - 0 | < , | - 0 | >

, . (32)

exp[-

2 |-0 | 3

4|-0 |[1-

) (-0 )2 2 ]1/2 -0 )2

3/2

It is easily seen that the spectrum (32) of the FM oscillation is exponentially small outside the region | - 0 | < and we assume further the following expression as the spectral density of the FM oscillation S0 , | - 0 | < , -0 2 ) 1-( (33) S () = . 0, | - 0 | The expression (33) can be derived also by the simple time averaging of oscillating spectral component of the harmonic signal
/

S() = 4

dt ( - 0 -
-/

cos t).

In Fig. 4(A) the Josephson oscillation line with the Lorentzian profile at D = 50 (curve 1), the convolution of the Josephson oscillation line with the FM spectrum (33) (curve 2) and the simulation of experimental data (hilbertogram) with signal-to-noise ratio S/N = 40 dB (curve 3) are presented. In Fig. 4(B) the RECOVERY result of FM spectrum restoration using the point spread function (24) at D = 50 (curve 1) together with the original FM spectrum (33) (curve 2) are presented. The superresolution ratio achieved is SR = 1.955 after 50 iterations.

Fig. 4. Frequency modulation test. (A) 1--the Josephson oscillation line with the Lorentzian profile D = 50, 2--convolution of the Lorentzian peak 1 with the FM spectrum Eq. (33) presented by curve 2 in (B), 3--the inverse Hilbert transform of the curve 2 with the normal noise S/N = 40 dB added and truncated on the right. (B) 1--the RECOVERY result of FM spectrum restoration using the PSF 2 at D = 50. The superresolution ratio achieved is SR = 1.955, 2 = 2.66, iteration number N = 50.


184

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

6. Real measurement results As a practical test of the integral equation approach we present the spectrum of frequency-modulated BWO measured by the Hilbert-spectroscopy technique. Measurements were carried out using high critical temperature Josephson junctions. The junctions were fabricated by laser ablation of YBaCuO on MgO Yttrium Stabilized Zirconia (YSZ) or sapphire bicrystal substrates and integrated with log-periodic broadband planar antennas. Typical parameters of submicrometer-wide junctions are: normal resistance up to 20 Ohm, critical current about 0.1 mA, characteristic voltage up to 2 mV at liquid helium temperature 4.2 K. Samples with Josephson junctions and antennas were placed on the flat side of an extended hyperhemisphere MgO lens in a cryostat with an optical window. The detector response was measured using the usual technique of a chopper for amplitude modulation of the input signal and a lock-in amplifier for extraction of the detector response. The backward-wave oscillator was used as a source of submm-wave radiation in the range 200í550 GHz. The linewidth of the BWO oscillations is below 1 MHz and for producing a broadband signal we use the harmonic deviation of the cathode voltage of the BWO. The deviation amplitude of 100 V corresponds to the frequency deviation of about 6 GHz at a frequency of around 410 GHz. The total frequency deviation from minimum to maximum frequency equals 12 GHz. Fig. 5 shows results of the sequential stages of data treatment in Hilbert-transform spectroscopy. The hilbertograms of the current-driven Josephson junction to the incident radiation generated by BWO without frequency modulation (thin curve 1) and with the frequency modulation (thick curve 2) are presented in Fig. 5(A). The Hilbert transforms obtained by means of the solution of the integral Eq. (19) with program Dconv2_n are shown in Fig. 5(B) also as thin and thick curves. At first glance, nothing interesting is seen. The smooth thin and thick spectral curves in this graph correspond to the left-hand of Eq. (16) for the monochromatic and FM radiation of BWO, respectively. Due to the large bandwidth of the Josephson oscillations

Fig. 5. Frequency modulation real measurement test. (A) 1--input hilbertogram of BWO spectrum at 410 GHz with NO frequency modulation (thin curve), 2--the same as the curve 1 but with the frequency modulation (thick curve). Signal-to-noise ratio for input data S/N = 39.2 dB. (B) Thin line--the Hilbert transform of curve 1 from (A), thick curve--the Hilbert transform of the curve 2 from (B), twin peak in the center is the FM spectrum obtained by the RECOVERY deconvolution program Dconv_n from the input data shown in (A). Superresolution ratio achieved SR = 1.23 after 200 iterations.


E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

185

measured with non-modulated BWO (see the thin curve) there are no traces of the expected complicated spectrum of FM BWO in the thick curve, in accordance with the convolution integral in Eq. (16). In order to recover the superresolution FM spectrum of BWO modulation, that is to remove to some extent the effect of large bandwidth of the Josephson oscillations, we can use the measured spectrum of the Josephson oscillations (the thin curve in Fig. 5(B)) as the kernel (PSF) in the convolution integral Eq. (16). It turns out, however, that a more successful result can be achieved if we take the hilbertogram measured with monochromatic BWO (the thin curve in Fig. 5(A)) as the kernel of the integral equation (15) where the hilbertogram for FM BWO (the thick curve in Fig. 5(A)) should be used as the left-hand side of Eq. (15). It is worthwhile describing the last procedure in more detail. We can recover the invisible structure hidden in curve 2 in Fig. 5(A) by taking curve 1 from this figure as the PSF function and curve 2 also from this figure as the input data file for the standard deconvolution program Dconv_n from the RECOVERY program package. The recovery result is presented as a twin peak in the center of Fig. 5(B). The distance between two peaks of this twin peaks measured from the plot is equal to 11 GHz, that is in good agreement with the value of 12 GHz of the calibration measurement given above. The comparison of this result with the simulation result in Fig. 4 give us the proof of reliability of the method for data recovery in Hilbert-transform spectroscopy proposed in this paper.

7. Conclusion We have presented an approach to Hilbert transform spectroscopy data recovery as the solution an integral equation of the first kind by a nonlinear method implemented with the RECOVERY program package. This approach solves two problems which are essential in Hilbert transform spectroscopy: the limitation of the measurement interval (extrinsic or technical resolution limit) and the ability to increase the resolving power of the spectroscopy including superresolution reconstruction. The efficiency of this approach is demonstrated by numerical simulations and by the specially designed physical measurement of FM backward-wave oscillator spectrum recovered by this approach. The superresolution ratio achieved is SR = 2 in simulation tests (see Figs. 3 and 4) and SR = 1.23 in real measurements (see Fig. 5). We plan to publish the modifications to the RECOVERY deconvolution programs [8] used in this paper.

Acknowledgements The authors are very grateful to Professor Z. Ivanov and Professor T. Claeson of Chalmers University of Technology for useful discussions and to Dr. E. Stepantsov of Institute of Crystallography Russian Academy of Sciences for the fabrication of bicrystal substrates and samples. We thank Dr. E. Pantos of CLRC, Daresbury Laboratory, UK, for the interest in this work and for his assistance with the preparation of the manuscript. The work was partially supported by INTAS Grants 97-0731 and 01-686. One author (ELK) thanks the Program "Universities of Russia" for support under project No. 0150301028.

References
[1] Y.Y. Divin, O.Y. Polyanskii, A.Y. Shul'man, Incoherent radiation spectroscopy by means of the Josephson effect, Sov. Tech. Phys. Lett. 6 (1980) 454í455. [2] M.A. Tarasov, A.Y. Shul'man, G.V. Prokopenko, V.P. Koshelets, O.Y. Polyansky, I.L. Lapitskaya, A.N. Vystavkin, E.L. Kosarev, Quasioptical Hilbert transform spectrometer, IEEE Trans. Appl. Superconductivity 5 (2) (1995) 2686í2689, Part 3. [3] J.H. Hinken, in: V. Kose (Ed.), Superconducting Quantum Electronics, Springer-Verlag, 1989, p. 151.


186

E.L.Kosarev et al. / Computer Physics Communications 151 (2003) 171í186

[4] Y.Y. Divin, H. Schulz, U. Poppe, N. Klein, K. Urban, V.V. Pavlovskii, Millimeter-wave Hilbert-transform spectroscopy with high-Tc Josephson junctions, Appl. Phys. Lett. 68 (11) (1996) 1561í1563. [5] A. Kaestner, M. Volk, F. Ludwig, M. Schilling, J. Menzel, YBa2 Cu3 O7 Josephson junctions on LaAlO3 bicrystals for terahertz-frequency applications, Appl. Phys. Lett. 77 (19) (2000) 3057í3059. [6] M. Tarasov, E. Stepantsov, Z. Ivanov, A. Shul'man, O. Polyansky, A. Vystavkin, M. Darula, O. Harnack, E. Kosarev, D. Golubev, Methods of submillimeter-wave Josephson spectroscopy, in: Inst. Phys. Conf. Ser., IOP Publ., No. 167, 2000, pp. 611í614. [7] E.L. Kosarev, Shannon's superresolution limit for signal recovery, Inverse Problems 6 (1) (1990) 55í76. [8] V.I. Gelfgat, E.L. Kosarev, E.R. Podolyak, Programs for signal recovery from noisy data using the maximum likelihood principle, Comput. Phys. Comm. 74 (3) (1993) 335í357. [9] S.G. Rautian, Real spectral apparatus, Sov. Phys.--Uspekhi 66 (1) (1958) 245í273. [10] E.L. Kosarev, E.P. Krasnoperov, A new calculation technique of muonium formation rate, Comput. Phys. Comm. 126 (1í2) (2000) 93í100. [11] E.A. Kolganova, E.L. Kosarev, G.A. Ososkov, Superresolution algorithms for data analysis of discrete detectors in nuclear physics, Nucl. Instrum. Methods Phys. Res. A 443 (2000) 464í477. [12] J.-A. Conchello, Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images, J. Opt. Soc. Amer. A 15 (10) (1998) 2609í2619. [13] K.K. Likharev, B.T. Ulrich, Systems with Josephson Contacts, Moscow State Univ. Press, 1978 (in Russian). [14] K.K. Likharev, An Introduction to Dynamics of Josephson Junctions, Nauka, Moscow, 1985 (in Russian). [15] A. Barone, G. Patern?, Physics and Applications of the Josephson Effect, Wiley, 1982. [16] H. Kanter, F.L. Vernon, High-frequency response of Josephson point contacts, J. Appl. Phys. 43 (7) (1972) 3174í3183. [17] R.L. Kautz, R.H. Ono, C.D. Reintsema, Effect of thermal noise on Shapiro steps in high-Tc Josephson weak links, Appl. Phys. Lett. 61 (3) (1992) 342í344. [18] Y.Y. Divin, J. Mygind, N.F. Pedersen, P. Chaudhari, Josephson oscillations and noise temperatures in YBa2 Cu3 O7-x grain-boundary junctions, Appl. Phys. Lett. 61 (25) (1992) 3053í3055. [19] D. Koelle, R. Kleiner, F. Ludwig, E. Dantsker, J. Clarke, High-transition-temperature superconducting quantum interference devices, Rev. Mod. Phys. 71 (3) (1999) 631í686. [20] A. Y. Shul'man, to be published elsewhere. [21] Y.Y. Divin, O.Y. Polyanskii, A.Y. Shul'man, Incoherent radiation spectroscopy based on a.c. Josephson effect, IEEE Trans. Magnetic MAG-19 (3) (1983) 613í615. [22] L.D. Landau, E.M. Lifshits, Electrodynamics of Condensed Matter, Nauka, Moscow, 1970 (in Russian). [23] R. Klucker, U. Nielsen, KramersíKronig analysis of reflection data, Comput. Phys. Comm. 6 (4) (1973) 187í193. [24] S.J. Collocott, Numerical solution of KramersíKronig transforms by a Fourier method, Comput. Phys. Comm. 13 (3) (1977) 203í205. [25] S.J. Collocott, G.J. Troup, Adaptation: numerical solution of the KramersíKronig transforms by trapezoidal summation as compared to a Fourier method, Comput. Phys. Comm. 17 (4) (1979) 393í395. [26] J.A.C. Weideman, Computing the Hilbert transform on the real line, Math. Comp. 64 (210) (1995) 745í761. [27] M.Z. Tarasko, On the one method for solution of the linear system with stochastic matrices, Preprint Physics and Energetics Institute, Obninsk, PEI-156, 1969 (in Russian). [28] A.B. Kuz'menko, E.A. Tishchenko, A.S. Krechetov, Development of the KramersíKronig method for spectra of polarized IR radiation reflected from the surface of low-symmetry crystals, Optics Spectrosc. 84 (3) (1998) 402í409. [29] E.L. Kosarev, E. Pantos, Optimal smoothing of noisy data by fast Fourier transform, J. Phys. E.: Sci. Instrum. 16 (1983) 537í543. [30] A.A. Kharkevich, Spectra and Analysis, GIFM, Moscow, 1962 (in Russian). [31] H. Bateman, A. Erdùlyi, in: Higher Transcendental Functions, Vol. 2, McGraw-Hill, 1953, Section 7.4.4.