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

Journal of Communications Technology and Electronics, Vol. 48, No. 10, 2003, pp. 1124­1136. Translated from Radiotekhnika i Elektronika, Vol. 48, No. 10, 2003, pp. 1225­1237. Original Russian Text Copyright © 2003 by Shul'man, Kosarev, Tarasov. English Translation Copyright © 2003 by MAIK "Nauka / Interperiodica" (Russia).

RADIO PHENOMENA IN SOLIDS AND PLASMA

Hilbert-Transform Spectroscopy Based on the a.c. Josephson Effect. Theory and Computational Technique
A. Ya. Shul'man, E. L. Kosarev, and M. A. Tarasov
Received May 11, 2003

1

Abstract--The limitations of the Hilbert-transform spectroscopy theory involving a resistive­shunted model of the Josephson junction are considered. It is noted that these limitations are especially important when Hilbert-transform spectroscopy is used to analyze the spectrum of Josephson oscillations. Formulas extending the original expressions of Hilbert-transform spectroscopy to a non-Lorentzian form of the Josephson oscillation spectrum are presented. The computational technique applied in Hilbert-transform spectroscopy to recover the radiation spectrum under study from the input data (hilbertogram) via the discrete Fourier transform is analyzed. A method of measured data treatment is stated, in which the Hilbert transformation is considered as an integral equation of the first kind for the radiation spectrum. The integral equation is solved using the maximum likelihood method (one of the methods applied to solve ill-posed problems), which allows one to attain the maximum possible resolution enhancement at a given signal-to-noise ratio. The possibility to resolve spectral peculiarities in the measured spectrum with a typical scale less than the Josephson oscillation linewidth is shown via simulation and measurement of the frequency-modulated (FM) radiation spectrum produced by a backwardwave oscillator (BWO).

INTRODUCTION As it was shown in [1], variations of the I­V characteristic of the Josephson junction affected by broadband incoherent electromagnetic radiation are related to the spectrum of this radiation by the Hilbert transformation. Based on this relation, the submillimeter spectroscopy technique was proposed. It was later called Hilbert-transform spectroscopy [2]. Since 1980, most studies in this area have dealt with the selection of the type of a Josephson junction and quasi-optical elements coupling it with the radiation that are the most suitable for this spectroscopy (see, for example, [3­6]). The Hilbert-transform spectroscopy technique has recently been extended from the gigahertz to terahertz band due to the advent of Josephson junctions based on hightemperature superconductors (HTSCs), which are sensitive to radiation in this band [7]. However, the theory has not been noticeably developed, and numerical and purely spectroscopic aspects of the new spectral measurement method have not been comprehensively analyzed. The purpose of this work is to fill this gap to some extent. Hilbert-transform spectroscopy was incited by the similarity between an odd-resonance shape of the selective response of the Josephson junction to the monochromatic irradiation as a function of the bias voltage and the frequency dependence of the refraction index observed in the anomalous dispersion region,
1

1

In this paper, English scientific terms are suggested by the authors. (Ed.)

which is well-known in optics. Qualitative analysis shows that the Hilbert transformation applied to the curve of such a shape yields a peak whose position on the bias axis depends on the frequency of the radiation incident on the junction. However, the quantitative description of the result of this transformation and substantiation of Hilbert-transform spectroscopy became possible only on the basis of formula (6.49) derived in [8]. In this case, it turned out that insensitivity of the radiation spectrum obtained using the Hilbert transformation to a discrepancy of the real junction from a resistive model is only ensured after the voltage response found in [8] is transformed into the current response used in [1] to derive the basic relationship of Hilbert-transform spectroscopy. These problems are described in detail in Section 1. Here, we only note that, for the determination of the spectrum of the external radiation affecting the junction, the details omitted in [1] are not important. However, the growing tendency to use the Hilbert-transform spectroscopy for determining the self-radiation spectrum of Josephson structures requires to introduce clarity into its basic relationships of Hilbert-transform spectroscopy (Section 1). The discrete Fourier transformation (DFT) is one of the routine numerical methods of applying the Hilbert transformation to measured curves. However, without taking into account specific features of function classes including a measured response and sought spectrum, the direct DFT results in qualitative distortions and quantitative errors of both recovered spectra and calcu-

1124


HILBERT-TRANSFORM SPECTROSCOPY BASED ON THE A.C. JOSEPHSON EFFECT

1125

lated (in the process of simulation) responses of a junction to the radiation with a given spectrum. The sources of these difficulties and possible workaround are described in the first part of Section 2. The Hilbert transformation, like the Fourier or Laplace transformations, is a one-to-one integral transformation. However, when the finite linewidth of the Josephson oscillation is taken into account, the deconvolution problem in Hilbert-transform spectroscopy can be considered most comprehensively by means of solution methods developed for integral equations of the first kind treating spectrum recovering as an ill-posed problem. The second part of Section 2 describes the derivation of the integral equation for the incident radiation spectrum. To solve this equation, we use the maximumlikelihood method realized in the RECOVERY software package [9], which makes it possible to obtain the maximum possible resolution enhancement for a given signal-to-noise ratio and even to attain superresolution, as shown in [10, 11] (Section 3). In Section 4, using numerical simulation and a measured response of a Josephson HTSC junction, we show that the solution of the integral equation allows one to recover the sought spectrum at a resolution exceeding the natural limit, i.e., to attain superresolution in accordance with the definition and criterion formulated in [10]. It is worthwhile to define more exactly the terms of deconvolution and superresolution with reference to the case considered. Since the effect of the point spread (instrumental) function of a spectroscopic instrument is usually described by the convolution integral (see, for example, [13]), the deconvolution procedure means the elimination of distortions of the measured spectrum caused by the point spread function. These distortions can be eliminated using various methods, but we preserve the deconvolution term for procedures which recover measured details that are distorted during measurements. It is possible to characterize these details using, for example, the known Raileigh resolution criterion for the original spectrum or the effective spectral width of its Fourier transform. Superresolution can be defined as a procedure of recovering details which appear to be invisible or suppressed below a noise level due to the spectrum transformation by the point spread function [10­12]. For Hilbert-transform spectroscopy, we call as a deconvolution any method which allows one to minimize spectrum distortions due to a limited measurement interval. By our definition, superresolution is a procedure eliminating to some extent the influence of the finite linewidth of the Josephson oscillation. The results of this work show that, in Hilbert-transform spectroscopy, the approach to the spectrum-recovering problem involving the integral equation makes it possible to simultaneously solve two problems: deconvolution and superresolution.

1. PRINCIPLES OF THE HILBERT-TRANSFORM SPECTROSCOPY Josephson Effect The detailed description of the Josephson effect can be found in [8, 14, 15]. We give only the basic relationships, which are required in the subsequent text. The superconducting state is characterized by a complex order parameter (the wave function of a superconducting condensate): = exp ( i ) , where is related to an energy gap for normal quasiparticle excitations and is the phase of the order parameter. The Josephson effect arises when two superconductors are separated by a thin layer (barrier), in which is suppressed, but tunneling remains still possible for Cooper pairs, i.e., the condensate wave function leaks by tunneling from one superconductor into the other. In this case, the superconducting current I (supercurrent) of Cooper pairs through the barrier can be written as I = I c sin , = L ­ R, (1)

where is now the phase difference of the left and right superconductors and Ic is the critical current of the d.c. Josephson effect. As long as I < Ic , i.e., the current of an external source does not exceed the critical one, there is no voltage drop across the Josephson junction and the charge is transferred through the junction only by the superconducting current component. In addition to the superconducting current produced by spatial phase nonuniformity, there are cases when the phase difference depends on time. In this case, voltage drop V across the Josephson junction is observed and we can employ the known Josephson relationships d ------ = 2 eV , dt V = 2 eV / . (2) (3)

At time-independent V, formulas (1) and (2) imply that the superconducting current oscillates at Josephson frequency V . This self-sustained Josephson oscillation serves as a basis of all spectroscopic applications of the Josephson effect. Nonlinearly interacting with an external time-dependent excitation, these oscillations cause changes in the d.c. I­V characteristic of the junction or produce a response at intermediate frequencies. An elementary model of the Josephson junction is the so-called model of a resistive­shunted junction (RSJ) or resistive model. In this case, the total junction current is represented as a sum of the superconducting
Vol. 48 No. 10 2003

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


1126

SHUL'MAN et al.

current and the current of normal quasi-particle excitations (electrons): V / R = I c sin = I , where R is the resistance of the junction in the normal (non-superconducting) state. This equation together with the Josephson equation (2) comprises the system which describes all features of small-sized junctions. It is conventional to use dimensionless units which are determined by the relationships: i = I / Ic, V c = Ic R, v = V /V c, = c t . c = 2 eV c / , = / c ,

Let us consider Kanter­Vernon relationship (6) and analyze whether it can generate the Hilbert transformation. At first glance, its right-hand side looks like the kernel of the Hilbert transformation. If we could describe the junction response to the a.c. current with continuous spectrum Si() by the formula


i ( v ) P . V . d S i ( ) i ( v , ) ,
0



Then, it is possible to rewrite the set of equations in the dimensionless form: d / d + sin = i , d/d = v . (4)

where P.V. is the Cauchy principal value of the integral, we might hope to obtain spectrum S() of the radiation incident on the junction via the inverse Hilbert transformation. However, the response becomes linear rather than quadratic with respect to the a.c. amplitude in the vicinity of the resonance voltage v = . Therefore, the situation should be considered in more detail. Theory of Hilbert-Transform Spectroscopy Kanter­Vernon relationship (6) is obtained by neglecting fluctuation currents if () in the junction. In practice, these fluctuations are always present and must be included into the expression for the current. Let us rewrite the first equation of system (4) as d / d + sin = i + i f . (7)

If the d.c. bias voltage on the junction is ensured, the d.c. (time-averaged current) I­V characteristic will be the simple linear relation i = v, since, in this case, the supercurrent is a purely harmonic function of time. With the direct current through the junction (currentdriven mode), the solution to Eqs. (4) gives a nonlinear I­V characteristic, which can be written as i = sgn v v + 1 ,
2

(5) A solution to Eq. (7) yields a finite width for all harmonics of the Josephson oscillation spectrum and eliminates divergence in formula (6) [8]. The analytical expression for the response to a monochromatic signal has the form [1]: ~ 2() i v + v ­ i (v , ) = ----------- -------------------------------- + ------------------------------- . (8) 2 2 2 2 8 i v (v + ) + (v ­ ) + Here, it is assumed that the fluctuation current is -correlated in time and the conditions v , are met, where is the fluctuation-induced linewidth of the Josephson oscillation. Formula (8) for the current response follows from formula (6.49) in [8] for the voltage response (the latter was divided by an expression for the differential junction resistance calculated within the framework of the resistive model). Certainly, if real junctions were exactly described by the resistive model, both of the expressions would be absolutely equivalent in respect to Hilbert-transform spectroscopy. However, formula (8) is less sensitive to deviations of the behavior of I­V characteristics of real junctions from the ideal one described by (5). One reason for this is the insignificant contribution of variations of the quasi-particle current to the total response as compared to that of variations of the supercurrent at the zero frequency at a given d.c. junction voltage V .
Vol. 48 No. 10 2003

where v is the time-averaged bias voltage and the overbar means the time averaging. It should be noted that, in this case, the voltage across the junction is a periodic (but not harmonic) function of time. Its fundamental frequency is determined by Josephson relationship (3) for the time-averaged bias voltage v = v . Now, let an additional monochromatic external curi rent ~ () at an angular frequency be supplied to the junction. Then, the I­V characteristic changes. The corresponding expression describing the current response at a given voltage v = const was obtained by Kanter and Vernon [16] in a quadratic approximation:
2 1~ i -2 i ( v , ) = ---- ----------------- . 4 i v ­ 2

(6)

It is clear that this expression is invalid in the vicinity of the "resonance" bias, where v = . Detailed analysis shows that a current step (Shapiro step) is formed at this point on the I­V characteristic. Its value linearly i depends on the first degree of the AC amplitude ~ (). In this d.c. interval, the frequency of the Josephson oscillation at the step remains synchronized with the external frequency (see, for example, [8]).

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


HILBERT-TRANSFORM SPECTROSCOPY BASED ON THE A.C. JOSEPHSON EFFECT

1127

Now, we can find that the a.c. current with spectral density Si() leads to the response described by the formula [2]


i (v ) =



0

d i ( v , ) (9)

1 ­v -= --------- ­ -- d S i ( ) ------------------------------- . 2 2 8iv ( ­ v ) + ­



It is necessary to note that various types of antennas are usually used to couple radiation with the Josephson junction. Thus, the following relationship between a.c. current spectrum Si() and the spectrum of electromagnetic radiation affecting a Josephson junction holds: Si ( ) = K ( ) S ( ) ,
2

where K() is the transfer function of the antenna. To determine the true signal spectrum, it is necessary to take into account this function. It can be seen that the expression in brackets in Eq. (9) is the Hilbert transform of function Si() as the limit 0 is taken. By introducing the function 8 g ( v ) = -- i v i ( v ) (10)

tribution in the integrand can be replaced by (' ­ ), ^ and obtained function S () becomes equal to sought spectrum Si(). This procedure formally corresponds to the passage to the limit 0 in expressions (9) and (12). At the same time, if the function Si() is more narrowband than the spectrum of the Josephson oscillation and one can assume that it is proportional to ( ­ 0), the latter expression yields intensity of the Josephson oscillations J(0 ­ v , v ) at frequency 0 as a function of bias voltage v . It is important to note that the linewidth of the Josephson oscillation actually depends on the bias voltage. This fact is explicitly shown in Eq. (13) by the second argument of function J(0 ± v , v ). Within the framework of the resistive model, fluctuation current if () originates from equilibrium thermal fluctuations in normal resistance R of the junction and, in accordance with [8], the oscillation linewidth is expressed (in the dimensional form) as follows: 1 ( V ) ( v ) c = -
2 22 Rd ( V ) Ic 2 e -------------- kT 1 + ---------------- . ----- 2 R 2 I ( V )

and applying to it the Hilbert transformation in the form ^ S() = H
v

1 g(v ) [ g ( v ) ] -- P . V . dv ------------- , v ­
­





(11)

we obtain (see also [2]) 1 ^ S ( ) = -- d' S i ( ' ) -------------------------------- . 2 2 ( ' ­ ) + ­




(12)

The spectrum Si() estimate is designated in ^ expressions (11) and (12) as S (). In formula (12), it is easy to recognize the convolution of the spectral density of the induced alternating current with the Lorentz distribution. The latter is the positive-frequency part J( ­ v , v ) of the spectrum of intrinsic Josephson oscillations (see [8, formula (6.63)]) (v ) (v ) S J ( , v ) ---------------------------------------- + ----------------------------------------2 2 2 2 ( ­ v ) + ( v ) ( + v ) + ( v )(13) J ( ­ v , v ) + J ( + v , v ), which was obtained theoretically in the resistive model with -correlated thermal noise. If Si() slowly varies over a frequency interval of order of , the Lorentz dis-

Here, Rd = d V /d I is the differential resistance of the junction. It can be seen that depends on the bias voltage. This circumstance makes the transfer from (9) to (12) via Hilbert transformation (11) questionable even within the scope of the resistive model. Below, we explain under what conditions the dependence of on v does not hinder applying Hilbert-transform spectroscopy. We should also note that the characteristics of real Josephson junctions may differ from those of the ideal resistive model. On the one hand, it is known that the resistive model adequately describes certain Josephson junctions. Its applicability for junctions made of lowtemperature superconductors is discussed in [14]. In the case of HTSC junctions, the application of the RSJ model to description of junction high-frequency and noise characteristics [4, 6, 17, 18] and squids [19] has been successfully tested in experiments. On the other hand, it is important to know to what extent the results obtained by the Hilbert-transform spectroscopy method depend on the fact that the Josephson junction is not ideal. This problem is the subject of current studies, especially investigations of HTSC junctions or more intricate Josephson structures, such as coherent chains of Josephson junctions. In particular, it was only supposed in [2] that expression (12) is suitable for measuring the spectrum of the Josephson oscillation. However, it had been already applied to the quantitative analysis of the Josephson oscillation spectrum [4, 6] without any justification. Solving Eq. (7) for a noisy junction exposed to the monochromatic radiation by the method similar to that
Vol. 48 No. 10 2003

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


1128

SHUL'MAN et al.

used in [20], it is possible to derive the generalization of Eq. (9), which allows one to eliminate limitations of Hilbert-transform spectroscopy by the Lorentz shape of the Josephson oscillation spectrum. It can be easily seen that, taking into account Eq. (13), (9) can be represented as i (v ) 1 J ( ' ­ v , v ) (14) 1 -= --------- -- d S i ( ) -- P . V . d' ----------------------------- . ' ­ 8iv ­ ­






In turn, expression (14) is modified to the form i (v )


(15)

1 1 J ( ­ v ', v ) -= --------- -- d S i ( ) ­ -- P . V . dv ' ----------------------------- . v' ­ v 8iv ­ ­





When ( v ) / v 1 over the entire essential bias voltage range, the second argument of spectral density J in (15) does not prevent us from calculating the Hilbert transformation in accordance with Eq. (11) for an arbitrary (but smooth) dependence of on v . This statement can be proved by using the Bedrosyan theorem about "freezing" slow factors involved in the Hilbert transformation [21] (see also [22, 23]). As a result, we generalize the convolution in expression (12) in the form


1 ^ S ( v ) = -- d S i ( ) J ( ­ v , v ) .
­



(16)

surement procedure. Function g( v ) in Eq. (10) must be formed of measured data. To this end, I ( V ) characteristic and response I ( V ) must be measured as functions of bias voltage V . The bias voltage sweep range must include the interval [0­ V max ] with sufficiently high maximum voltage V max in order to provide for effectively infinite limits of integration with respect to voltage in expression (11). We introduce the term hilbertogram for function g( v ) by analogy with the term interferogram in Fourier spectroscopy. Resolution in Hilbert-transform spectroscopy depends on linewidth of the Josephson oscillation , the length of the voltage sweep interval, and the spacing of voltage sampling points disposed usually equidistantly. The linewidth depends on intrinsic voltage fluctuations in the junction and the external noise level. The latter should be suppressed to values as low as possible over a wide frequency range: from the line frequency to TVstation frequencies. The experimental setup must make it possible to modulate the input radiation and measure the phase I of response using a lock-in amplifier. Function g( v ) determined by Eq. (10) can be obtained from the measured data by a straightforward procedure [26]. The simplest case takes place for the d.c. voltage-driven junction, when two curves (the current and response) should be measured as functions of the bias voltage. For the current-driven junction, it is necessary to calculate the current response from the measured voltage response using the known expression I = ­ V / Rd . (17)

It follows from (16) that Hilbert-transform spectroscopy can be used to measure both the spectrum of the radiation incident on a junction and the spectrum of the Josephson oscillation itself, regardless of the actual spectral distribution of the latter, if one of the spectral distributions involved is considerably narrower than the other. The proof of this statement is similar to those used in optical and noise spectroscopy (see, for example, [24, 25]). Expression (16) allows one to find the conditions such that Hilbert-transform spectroscopy can be used to experimentally measure the Josephson oscillation spectrum or the dependence of H on V and to determine the scheme of these experiments. If the radiation spectrum of an external source is measured and it is necessary to take into account a finite linewidth of the Josephson oscillation, expressions (15) and (16) make it possible to solve the deconvolution problem in the functional space of either measured or Hilberttransformed functions. Measurements in Hilbert-Transform Spectroscopy For clarity of the following data transformations involved in calculations, let us briefly describe the mea-

Differential resistance Rd( V ) can be directly measured or numerically calculated from the measured I­V characteristic. The need for expression (17) was discussed after formula (8). Figure 1 illustrates all necessary input data for Hilbert-transform spectroscopy. The hilbertogram was obtained from these data using expressions (10) and (17). The sample was an HTSC Josephson junction irradiated by a backward-wave oscillator (BWO) at 410 GHz. 2. DATA PROCESSING IN HILBERT-TRANSFORM SPECTROSCOPY Discrete Fourier Transformation and Hilbert-Transform Spectroscopy Due to the simplicity of the Hilbert transformation of trigonometric functions, H H H
xy

[ cos ax ] = ­ sin ay , a > 0; (18)

xy

[ sin a x ] = cos a y ,

xy

[ exp ( i ax ) ] = i sgn ( a ) exp ( i ay ) ,
Vol. 48 No. 10 2003

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


HILBERT-TRANSFORM SPECTROSCOPY BASED ON THE A.C. JOSEPHSON EFFECT

1129

it can be readily analytically applied to functions represented by the Fourier series or integral. A measured (at discrete points) hilbertogram can be represented by a truncated Fourier series whose coefficients are calculated via the DFT. After that, the Hilbert transform of this series can be found. However, for this method could provide for a satisfactory representation of the spectrum of the studied radiation, it proved necessary to continue the bias voltage sweep until measured values of the hilbertogram decrease to the measurement noise level. Because of a slow (as 1/ V ) falloff of the response outside the region where the major intensity of the measured spectrum was concentrated, we had to increase the voltage up to values significantly exceeding Vc , which resulted in a substantial increase in current noise, or a damage of the junction by large currents, and, in any case, reduced the upper boundary of spectroscopy operating range for the given junction. We can process the curve, whose measurement is stopped, when it still has noticeably nonzero values. Due to the periodic continuation of the function by its discrete Fourier series, a discontinuity of the first kind arises at the boundaries of periods, which, as a result of the Hilbert transformation, is converted into the logarithmic divergence slowly decreasing from the point of discontinuity and yielding a large error even far from the period (measurement interval) boundaries. Since the response function is odd with respect to the change of the bias voltage sign, it is possible to reduce this discontinuity effect by odd-extending the response measured on the positive bias voltage section to the negative bias voltages. This method is equivalent to applying the discrete sinus Fourier transformation on the positive bias. However, in this case, the constant component of the recovered spectrum is lost and, in certain regions, the spectral density has inadmissible negative values. The trace of these processing defects often occurs in the literature. A typical illustration can be seen, for example, in [26, Fig. 1c]. Transfer to the Integral Equation with the Hilbert Kernel This work proposes another approach. By applying the inverse Hilbert transformation, we can invert Eq. (11) and rewrite it as a convolution integral of ^ sought spectrum estimate S () with the Hilbert kernel H
­1 v

I, µA; Rd, ; V, arb. units; g( V), arb. units 200

100

1

2 0 3

­ 100

4

0

500

1000 V, µV

1500

Fig. 1. An example of input experimental data as functions of the bias voltage: (1) current I, (2) differential resistance of the junction Rd , (3) voltage response V, and (4) hilbertogram g( V ).

to zero above a certain frequency and the integrals of them are finite. On the contrary, the Hilbert transform of these functions has long tails and slowly falls as 1/ V as v . The integral equation makes it possible to use a relative locality of the sought solution and avoid hilbertogram measurements in the voltage region distant from the spectrum existence domain. The measurements can be restricted to the voltage region where the hilbertogram value still exceeds noise (for further information see Section 4). To solve convolution integral equation (19) with the Hilbert kernel, the RECOVERY software package [9] was used. It is based on the maximum-likelihood method. As was shown in [10], this algorithm allows one to attain superresolution for a nonsingular kernel of the integral equation. This work demonstrates that this is also true for a singular kernel of the Hilbert transformation. Discrete Hilbert Transformation and Hilbert-Transform Spectroscopy In the communications theory and digital signal processing, as a rule, periodic or almost periodic functions with a bounded spectrum are applied. This spectrum can be described with a sufficient accuracy by the DFT coefficients if we correctly select a sampling frequency and avoid aliasing (see, for example, [27]). In this case, the Hilbert transformation can be analytically applied
Vol. 48 No. 10 2003

^ 1 S() ^ [ S ( ) ] ­ -- P . V . d ------------- = g ( v ) , ­v
­





(19)

^ which yields a singular integral equation for S (). Measured radiation spectrum S() and its estimate ^ S () usually belong to the class of functions with a compact support and finite energy; i.e., they are equal

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


1130

SHUL'MAN et al.

to a finite trigonometric series representing a measured function. When the Kramers­Kronig relation is used in optics to find absorption bands from the refractive index spectra or expression (11) is employed in Hilbert-transform spectroscopy, one has to deal with nonperiodic functions represented by the Fourier integral. In these cases, to eliminate contributions of periodic continuation of the measured function, one could try to follow the recommendations proposed in [23] (to complete the set of measured values with zeros until the measured interval doubles). However, if the measurements were stopped at a nonzero response value, it is impossible to avoid by this trick a discontinuity in the original function and, respectively, a logarithmic singularity in its Hilbert transform. However, in the case of Hilbert-transform spectroscopy, it is possible to make the following modification in the described method. If we extend the measured function beyond the right end as an even function, the function obtained in the double interval is continuous and vanishes at both ends, since, at the zero bias voltage v = 0, the function g( v ) vanishes by its definition (10). Then, applying the discrete sinus Fourier transformation and using the obtained coefficients in the inverse cosine Fourier transformation, we find the sought spectrum accurate to the constant component. The latter can be determined via extrapolation of the first several harmonics of the calculated DFT coefficients to the zero frequency. Certainly, this procedure ensures adequate results when the input function is measured in the entire interval where the sought function is nonzero. However, in this case, it is not necessary to measure the input function values up to the noise level. The procedure proposed cannot be applied for the inverse process, namely, for calculating function g( v ) ^ from spectrum S () when solving integral equation (19), since errors of calculated values g( v ) at the interval end are always large and slowly decrease as the interval increases. Therefore, the convolution in (19) is calculated using the discrete Hilbert transformation with the kernel digitized at half-integer argument values, i.e., at the points located in the middle between the points ^ where g( v ) and S () are calculated (see, for example, [28, 29]). In this case, the convolution integral can be calculated using the fast DFT. 3. DEFINITION OF SUPERRESOLUTION IN HILBERT-TRANSFORM SPECTROSCOPY There is a natural connection of Eq. (19) with the Kramers­Kronig relations in electrodynamics and optics [20], and there are algorithms and computer programs which are specially intended for that kind of calculations [31­34]. All methods used in these works can be referred to as linear methods, since they involve a linear data expansion in various sets of basis functions.

According to [10], any linear restoration method can only modify the amplitudes of Fourier harmonics rather than generate new ones, which are absent in initial data or lost in input noise. Therefore, the linear methods cannot guarantee superresolution. Now, let us define superresolution coefficient SR as the ratio of the Fourier spectrum width of the output signal to the Fourier spectral width of the input signal: SR = ( ) out / ( ) in . (20)

Quantity for any signal is determined by the formula


1 = ------------------ S ( ) d , S max ( )


0

(21)

where S() is the spectral signal intensity on the positive frequency semiaxis and Smax() is its maximum value. It should be noted that definition (20) does not formally coincide with that from [10], since the kernel of the Hilbert transformation 1 K ( x ­ y ) = ---------x­y (22)

has an infinite bandwidth both in the direct x-space and in the Fourier-transform -space. However, new definition (22) can be equally used both for the kernel of the Hilbert transform and ordinary Lorentz or Gaussian kernels. As it was shown in [10], the maximum possible superresolution is mainly determined by the experimental signal-to-noise ratio and can be readily predicted using the simple formula SR Es dB 1 = -- log 2 1 + ----- ------ , 10 3 En

max

where Es is the energy of the measured signal, En is the input noise energy, and dB = 10 log ( E s / E n ) is the input signal-to-noise ratio. By virtue of well-known formulas (18), the Hilbert transformation retains the values of the harmonics of any function. Hence, ( )
out

= ( ) in ,

and, according to our definition of superresolution (20), we should always obtain SR = 1 for the Hilbert transformation. This is a very important general conclusion. At first glance, we cannot obtain SR > 1 in Hilbert-transform spectroscopy. In fact, this is not the case. Actually,
Vol. 48 No. 10 2003

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


HILBERT-TRANSFORM SPECTROSCOPY BASED ON THE A.C. JOSEPHSON EFFECT

1131

we always have a finite bandwidth of Josephson oscillation J( ­ v , v ) in Eqs. (14)­(16), so, in addition to the Hilbert transformation in (11), we deal with the convolution integrals of input spectrum Si() either with the Hilbert transform of the spectrum of Josephson oscillations in (14) or with J( ­ v , v ) itself in (16). As a result, the spectral bandwidth of an actual kernel in the integral equations for Si() is always finite, and the concept of superresolution in Hilbert-transform spectroscopy is well defined. This circumstance enables one to increase the bandwidth of the output signal via the recovering procedure. In contrast to methods that involve the direct calculation of the Hilbert transform of the input experimental data, in this work, we propose to solve convolution integral equation (19). In the calculations, we used the Dconv2_n program, which is a modification (for details see [35]) of the Dconv2 program from the RECOVERY package [9] based on the maximum-likelihood method. This is a nonlinear method ensuring superresolution [10]. A detailed description of the recovery procedure is given in [9­11]. Work [12] describes the frequency-domain analysis of superresolution for the mean-maximization algorithm, which virtually coincides with our algorithm, the latter being a modification of the Tarasko algorithm [36]. We should also mention the work [37] devoted to infrared spectroscopy, in which the Kramers­Kronig relation is also reduced to a system of integral equations of the first kind solved by the nonlinear minimization method. 4. EXAMPLES OF DATA PROCESSING IN HILBERT-TRANSFORM SPECTROSCOPY Numerical Simulation Before applying the Dconv2_n program to actual experimental data, we tested the method using numerical simulations. As a rule, we can use only a finite (rather than infinite) bias interval over which hilbertogram g( v ) is measured. This interval has to be matched with the finite frequency band 1 2, 2 < (usually 1 = 0) where the measured spectrum Si() is located. We define this aspect as data limitation. Figure 2 illustrates the effect of such a limitation on the final result for various methods of the Hilbert transformation implementation. In one case, the method is based on the fast DFT combined with optimal filtering [38]. In another one, the nonlinear method realized in the Dconv2_n program is used. Figure 2a shows the original spectrum (curve 1), which is to be obtained using the Hilbert-transform spectroscopy technique, and the input hilbertogram with the 20-dB signal-to-noise ratio calculated from it (curve 2). The hilbertogram is truncated at about half its maximum value. The transformation results are illustrated in Fig. 2b. A drastic difference between the results yielded by both methods at

high-frequency spectrum edge can be seen. The implementation of the Hilbert transformation in accordance with formula (11), where the hilbertogram was represented by its own Fourier series using the fast DFT, led to a large deviation of this result from the correct behavior of the spectrum near the point 2 . At the same time, the solution to integral equation (19) did not exhibit such deviation, and the noise remained approximately at the same level as in the original signal. The false tail of the spectrum obtained by the DFT is the trace of the logarithmic divergence of the Hilbert transform at the interval end point where a discontinuity occurs. This discontinuity is due to the periodic repetition of the hilbertogram represented by the Fourier series. It is useful to note that the application of the even­ odd extension trick to the hilbertogram (see Section 2) gives a spectral curve which virtually coincides with the curve 1 (Fig. 2b). We have failed to obtain, using the same method, hilbertogram 2 from Lorentz distribution 1 (Fig. 2a). The reason is the localizability of the spectrum (the tail of the Lorentz distribution drops as 1/2) and the extension of the hilbertogram obtained from it (decreases as 1/). In addition, due to the oneto-one correspondence of the direct and inverse Hilbert transformations, the repeated application of the Hilbert transformation to g( v ) is equivalent to the convolution of Si() with the delta function. One more example of simulation is the test of Hilbert-transform spectroscopy using the spectrum of frequency-modulated (FM) oscillations. This case gives us an opportunity to simulate the effect of intricate radiation with a continuous spectrum on the junction. The time-dependent function describing oscillations with the harmonic frequency modulation can be written as f ( t ) = cos 0 t + --- sin t . (23)

The Fourier expansion of this almost periodic function can be found in many textbooks (see, for example, [39]), but, in our case, it is more convenient to deal with covariation f(t) and its spectral density. Using the definition of covariation function () and its spectral decomposition, which is suitable for periodic or almost periodic functions
T /2

1 ( ) = lim -T T 2 S ( ) = lim -T T () =


­ T /2 T /2

dtf ( t ) f ( t + ) ,


0

d ( ) cos ( ) ,

(24)

k ( ­ , )



S ( k ) exp ( ­ i k ) ,

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS

Vol. 48

No. 10

2003


1132 ­ S(), g( v ) 6 (a) 4

SHUL'MAN et al.

1 2 0 2 ­2 ­4

S(), S() 1 3 (b)

0 2

­1

0

Fig. 2. Data limitation effect in measurements: (a) (1) input Lorentz peak S() with the half-width D = 60 and (2) inverse Hilbert transformation g( v ) of the Lorentz peak with the normal noise added to the original signal at a signal-to-noise ratio of 20 dB (the ^ curve simulates the measured response) and (b) (1) original Lorentz distribution S(), (2) distribution S () recovered via calculation using formula (11) and the fast DFT with optimum filtering in calculations of the Hilbert transformation, and (3) distribution ^ S () obtained by solving integral equation (19) using the Dconv2_n program (thin noisy curve).

we obtain the following expression in the case of function (23): 1 2 ( ) = -- J 0 ------ sin ------ cos ( 0 ) , 2 2 12 S ( k ) = -- J k --- [ 4
mk

Here, Jk(x) is the Bessel function with an integer index; 0, > 0; (­, ); k = 0, 1, 2, ...; the elements of set {k} are determined as solutions to equations

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS



1

200

400 V, µV

600

800

1000

| ± 0| = k; the integer parameters are k ­ 0 m = -------------------- , k + 0 n = -------------------

and mn is the Kronecker delta: (25)
mn

+ nk ] .

1, = 0,

m=n m n.

We consider the case of a very low modulation frequency, when 0 / / 1.
Vol. 48 No. 10

(26)
2003


HILBERT-TRANSFORM SPECTROSCOPY BASED ON THE A.C. JOSEPHSON EFFECT

1133

Thus, the sum over k in (24) can be approximated by the integral with respect to as follows:
Original curves

() = 2


k > 0



­ J, g( v ) 2 1 0 ­1 ­2 0 200 400 600 800 1000 V, µV 2 1 3 (a)

S ( k ) cos k (27)

d 2 2 2 ------ ------- J µ --- + J --- cos ( ) . 2 2


0

Expression (27) determines the continuous approximation of the spectral density of an FM oscillation: 2 2 S ( ) = ------- J µ --- + J --- , 2 ­ 0 µ = ------------------ , + 0 = ------------------ .

(28)

S(), S() 3 2 1

1'

(b)

Under conditions (26), this formula can be simplified by using Langer's uniform asymptotic representation for Bessel functions ([40, Section 7.4.4]). Passing to the limit 0 in (28) and neglecting the second, exponentially small, term in brackets, we obtain S ( ) = 12 2 -- ( ­ ( ­ 0 ) ) 2
­ 1/2

2'

0

200

,

­ 0 <
2 3/2

400 600 800 (/2) â 0.4835, GHz

1000

2 ­ 0 exp ­ -- ------------------ 1 ­ ---------------------- 2 3 ( ­ 0 ) -----------------------------------------------------------------------------------1/2 2 4 ­ 0 1 ­ ---------------------- 2 ( ­ 0 ) ­ 0 > .

(29)

Fig. 3. Simulation of FM oscillation spectrum measurements: (a) (1) the Josephson oscillation line with the Lorentz profile and half-width D = 50, (2) convolution of FM spectrum S() (shown by curve 2') with the Lorentz distribution 1, and (3) hilbertogram g( v ) obtained using the inverse Hilbert transformation of curve 2 (normal noise is added to obtain the 40-dB signal-to-noise ratio) and ^ (b) (1') spectrum S () of the FM signal recovered using Eq. (16) and (2') original spectrum S() of FM oscillations given by the formula (30).

It can be easily seen that spectrum (29) of FM oscillations is exponentially small beyond the region | ­ 0| < , and below we employ the following expression for the spectral density of FM oscillations: S0 ------------------------------------ , ­ 0 < 2 0 ­ - S ( ) = 1 ­ -------------- 0, ­ . 0

(30)

Expression (30) can be derived by time averaging of the oscillating spectral component of the harmonic signal: S ( ) = ----4
/

Figure 3a shows the Lorentz-profile line of Josephson oscillations, convolution of the Lorentz distribution with FM spectrum (30), and simulated experimental data (hilbertogram) at the 40-dB signal-to-noise ratio. In Fig. 3b, curve 1' is the spectrum of FM oscillations recovered in two stages from hilbertogram 3 (Fig. 3). At the first stage, integral equation (19) with this hilbertogram serving as g( v ) is solved. The result is curve 2 (Fig. 3a), which does not demonstrate any signs of a complicated structure of the FM oscillation spectrum. Then, Eq. (16) is solved for Si(). As the kernel of the integral equation, the Lorentz profile is used (curve 1, Fig. 3a); as the right-hand side, we use the result of the previous stage (curve 2). After 50 iterations, the superresolution SR = 1.955 is achieved. Measurement Results To practically test the integral equation approach, we studied the spectrum of an FM BWO, which was
Vol. 48 No. 10 2003

­ /



dt ( ­ 0 ­ cos t ) .

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


1134 ­ g(V) 1000 (a) 2 0

SHUL'MAN et al.

­ 1000

1

500 S()

600

700

800

900

1000 1100 1200 V, µV

(b) 2000 1' 1000 2'

3

0 250 300 350 400 450 500 550 /2, GHz

Fig. 4. Experiment with an FM BWO: (a) (1) input hilbertogram g( V ) of the BWO spectrum at 410 GHz without frequency modulation; (2) the same with frequency modulation. The signal-to-noise ratio of the original signal is ^ 39 dB. (b) Curve 1' is Hilbert transform S () of curve 1a, ^ curve 2' is Hilbert transform S () of curve 2, and double peak at the center of curve 3 is the FM oscillation spectrum obtained using the Dconv_n program after 200 iterations (superresolution is 1.23).

measured using the Hilbert-transform spectroscopy method. In the measurements, we employed high-temperature Josephson junctions manufactured by the laser ablation of YBaCuO on MgO, Yttrium Stabilized Zirconia (YSZ), and sapphire bicrystal substrates and integrated with log-periodic broadband antennas. The typical characteristics of submicron-wide junctions measured in liquid helium are as follows: the normal resistance is as large as 20 , the critical current is ~100 µA, and the typical voltage is as high as 2 mV. Samples with Josephson junctions and antennas were placed on the flat surface of an extended hypersemispherical sapphire lens in a cryostat with an optical window.

The detector response was measured using a chopper for amplitude modulation of the input signal and a lock-in amplifier for extraction of the detector response. The BWO was used as a source of radiation in a range of 200­550 GHz. The intrinsic linewidth of the BWO radiation does not exceed 1 MHz. Therefore, for producing a broadband radiation, we used the harmonic modulation of the cathode voltage of the BWO. A deviation amplitude of 100 V corresponds to a frequency deviation of about 6 GHz at a central frequency of 410 GHz. The total frequency deviation from the minimum to maximum was about 12 GHz. Figure 4a shows sequential stages of data processing in Hilbert-transform spectroscopy: hilbertograms of current-driven Josephson junction correspond to the incident radiation generated by a BWO without (curve 1) and with (curve 2) frequency modulation. The Hilbert transforms are obtained by solving integral equation (19) whose right-hand side contains hilbertograms 1 and 2 (Fig. 4). At first glance, the results demonstrate nothing interesting. Spectral curves 1' and 2' (Fig. 4b) correspond to the left-hand side of Eq. (16) for the BWO monochromatic and FM spectra. Due to a large linewidth of the Josephson oscillation measured in the modulation-free mode of the BWO, one does not observe the expected complicated dependence of the original spectrum in modulation mode, which is in agreement with the convolution integral in (16). In order to recover an FM spectrum with superresolution, i.e., to eliminate, if possible, the influence of a large linewidth of the Josephson oscillation, one can use the measured spectrum of Josephson oscillations (BWO without frequency modulation, curve 1', Fig. 4b) as the kernel in convolution integral equation (16). However, the result (curve 3, Fig. 4b) proves to be more successful if we use the hilbertogram measured with the monochromatic radiation (curve 1, Fig. 4a) as the kernel in integral equation (15) and substitute the hilbertogram for the FM BWO (curve 2, Fig. 4a) on the left-hand side of Eq. (15) (after its respective transformation from i to g( v )). It makes sense to present the results above in a more conventional way. We can recover an invisible structure (from curve 2, Fig. 4a) using curve 1 as a point spread function and curve 2 as an input signal of the standard Dconv2_n program. The spectrum recovery result is represented by a forked peak at the center of Fig. 4b. The distance between two maximum values is 11 GHz. This is in good agreement with an estimate of 12 GHz from the BWO calibration described above. CONCLUSIONS The problems of Hilbert-transform spectroscopy application are discussed that are caused by using the resistive model of the Josephson junction as its foundaVol. 48 No. 10 2003

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


HILBERT-TRANSFORM SPECTROSCOPY BASED ON THE A.C. JOSEPHSON EFFECT

1135

tion. The generalization of basic relationships is suggested, which makes it possible to study the Josephson oscillation spectrum in structures whose I­V characteristic has features not described by the resistive model. One should bear in mind that the Josephson oscillation linewidth depends on the bias voltage and the most accurate results of Josephson oscillation spectrum measurements can be obtained only with the external radiation frequency scan at a fixed bias voltage (see formula (16) and in [1, formulas (7) and (8)]). The computational methods of the spectrum recovery in Hilbert-transform spectroscopy are analyzed. To improve spectrum measurement resolution and accuracy, it is proposed to apply the integral equation method rather than the DFT. The solution of the integral equation using the maximum-likelihood method is implemented in the RECOVERY program package. The efficiency of the proposed technique is illustrated by the example of the numerical simulation and spectrum measurements of the BWO FM radiation in the 400-GHz range. We attained the superresolution SR = 2 in the simulation and SR = 1.23 in the experiment. It should be noted that a number of fundamental problems arising in practical implementation of Hilbert-transform spectroscopy as a routine measurement method call for further studies, but they do not give grounds for doubts in its potentialities. ACKNOWLEDGMENTS We are grateful to Prof. T. Claeson for useful discussions of the work, T. LindstrÒm and Z. Ivanov (Chalmers University of Technology, Gothenburg, Sweden) for assistance in measurements, and E. Stepantsov (Institute of Crystallography, Russian Academy of Sciences, Moscow) for fabrication of bicrystal substrates and Josephson junctions. This work was supported by the INTAS (project no. 01-686). REFERENCES
1. Divin, Yu.Ya., Polyanskii, O.Yu., and Shul'man, A.Ya., Pis'ma Zh. Tekh. Fiz., 1980, vol. 6, no. 17, p. 1056 (Sov. Tech. Phys. Lett., 1980, vol. 6, p. 454). 2. Tarasov, M.A., Shul'man, A.Ya., Prokopenko, G.V., et al., IEEE Trans. Appl. Supercond., 1995, vol. 6, no. 2, part 3, p. 2686. 3. Hinken, J.H., Superconducting Quantum Electronic, Kose, V., Ed., Springer, 1989, p. 151. 4. Divin, Yu.Ya., Schulz, H., Poppe, U., et al., Appl. Phys. Lett., 1996, vol. 68, no. 11, p. 1561. 5. Tarasov, M., Shul'man, A., Polyansky, O., et al., Proc. ESA Symp. "The Far Infrared and Submillimeter Universe," Eur. Space Agency Symp. SP-401, 1997, p. 445. 6. Kaestner, A., Volk, M., Ludwig, F., et al., Appl. Phys. Lett., 2000, vol. 77, no. 19, p. 3057.

7. Tarasov, M., Stepantsov, E., Ivanov, Z., et al., Inst. Phys. Conf. Ser., 2000, no. 167, p. 611. 8. Likharev, K.K. and Ul'rikh, B.T., Sistemy s dzhozefsonovskimi kontaktami (Systems with Josephson Junctions), Moscow: Mosk. Gos. Univ., 1978. 9. Gelfgat, V.I., Kosarev, E.I., and Podolyak, E.R., Comp. Phys. Commun., 1993, vol. 74, no. 3, p. 335. 10. Kosarev, E.L., Radiotekh. Elektron. (Moscow), 1990, vol. 35, no. 1, p. 68 (Sov. J. of Communications Technology and Electronics (Scripta Technica, Inc.), 1990, vol. 35, no. 12, p. 90). 11. Kosarev, E.L., Inverse Problems, 1990, vol. 6, no. 1, p. 55. 12. Conchello, J.-A., J. Opt. Soc. Am., 1998, vol. 15, no. 10, p. 2609. 13. Rautian, S.G., Usp. Fiz. Nauk, 1958, vol. 66, p. 475 (Sov. Phys.-Usp., 1958, vol. 66(1), p. 245). 14. Likharev, K.K., Vvedenie v dinamiku dzhozefsonovskikh perekhodov (Introduction to Dynamics of Josephson Junctions), Moscow: Nauka, 1985. 15. Barone, A. and Paterno, G., Physics and Applications of the Josephson Effect, New York: Wiley, 1982. Translated under the title Effekt Dzhozefsona. Fizika i primeneniya, Aslamazov, L.G., Bulaevskii, L.N., and Vedeneev, S.I., Eds., Moscow: Mir, 1984. 16. Kanter, H. and Vernon, F.L., J. Appl. Phys., 1972, vol. 43, no. 7, p. 3174. 17. Kautz, R.L., Ono, R.H., and Reintsema, C.D., Appl. Phys. Lett., 1992, vol. 61, no. 3, p. 342. 18. Divin, Yu.Ya., Mygind, J., Pedersen, N.F., and Chaudhari, P., Appl. Phys. Lett., 1992, vol. 61, no. 25, p. 3053. 19. Koelle, D., Kleiner, R., Ludwig, F., et al., Rev. Mod. Phys., 1999, vol. 71, no. 3, p. 631. 20. Volkov, A.F., Radiotekh. Elektron. (Moscow), 1972, vol. 17, no. 12, p. 2581. 21. Bedrosyan, E., Proc. IEEE, 1963, vol. 51, p. 5. 22. Vainshtein, L.A. and Vakman, D.E., Razdelenie chastot v teorii kolebanii (Frequency Separation in Theory of Oscillation), Moscow: Nauka, 1983. 23. Bendat, J.S. and Piersol, A.G., Engineering Application of Correlation and Spectral Analysis, New York: Wiley, 1980. 24. Divin, Yu.Ya., Polyanskii, O.Yu., and Shul'man, A.Ya., Opt. Spektrosk., 1979, vol. 47, no. 1, p. 170. 25. Shul'man, A.Ya., Zh. Eksp. Teor. Fiz., 1981, vol. 81, no. 2, p. 784 (Sov. Phys. JETP, 1981, vol. 54(2), p. 420). 26. Divin, Yu.Ya., Polyanskii, O.Yu., and Shul'man, A.Ya., IEEE Trans. Magn., 1983, vol. 19, no. 3, p. 613. 27. Khurgin, Ya.I. and Yakovlev, V.P., Finitnye funktsii v fizike i tekhnike (Finite Functions in Physics and Engineering), Moscow: Nauka, 1971. 28. Rabiner, L.R. and Gold, B., Theory and Application of Digital Signal Processing, Englewood Cliffs: Prentice Hall, 1975. 29. Soroko, L.M., Gil'bert-optika (Hilbert Optics), Moscow: Nauka, 1981.
Vol. 48 No. 10 2003

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS


1136

SHUL'MAN et al. 36. Tarasko, M.Z., About Solving Technique for Linear System with Stochastic Matrixes, Preprint of Inst. Phys. Power Eng., Obninsk, 1969, no. 156. 37. Kuz'menko, A.B., Tishchenko, E.A., and Krechetov, A.S., Optics and Spectroscopy, 1998, vol. 84, no. 3, p. 402. 38. Kosarev, E.L. and Pantos, E., J. Phys. E. Sci. Instrum., 1983, vol. 16, p. 537. 39. Kharkevich, A.A., Spektry i analiz (Spectra and Analysis), Moscow: Fizmatgiz, 1962. 40. Bateman, H. and Erdelyi, A., Higher Transcendental Functions, New York: McGraw-Hill, 1953, vol. 2.

30. Landau, L.D. and Lifshitz, E.M., Electrodynamics of Continuous Media, Oxford: Pergamon, 1960. 31. Klucker, R. and Nielsen, U., Comput. Phys. Commun., 1973, vol. 6, no. 4, p. 187. 32. Collocott, S.J., Comput. Phys. Commun., 1977, vol. 13, no. 3, p. 203. 33. Collocott, S.J. and Troup, G.J., Comput. Phys. Commun., 1979, vol. 17, no. 4, p. 393. 34. Weideman, J.A.C., Math. Comput., 1995, vol. 64, no. 210, p. 745. 35. Kosarev, E.L., Shul'man, A.Ya., Tarasov, M.A., and Lindstroem, T., Comput. Phys. Commun., 2003, vol. 151, p. 171.

JOURNAL OF COMMUNICATIONS TECHNOLOGY AND ELECTRONICS

Vol. 48

No. 10

2003