1669

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. I, JULY 1994

Delayed Block Transfer Functionin the Frequency Domain
Hiroshi Kanai, Member, IEEE, Toyohiko Hori, Noriyoshi Chubachi, Member, IEEE, and Takahiko Ono

Abstract- When the impulse response of a transfer system islong, it is difficult for the standard cross spectrum method to
obtain an accurate estimate of the transfer function by applying
the fast Fourier transform (FFT) with a finite length window. The
bias error in the estimate is large especially around the resonant
frequency of the transfer system. In this paper, therefore, we
propose an alternative new method to obtain an accurate estimateof the transfer function. The delayed block transfer function is
introduced to detect the components that are correlated to thesignal in the input window but leak from the output window.
Based on these transfer functions, the total characteristics of thetransfer system are estimated accurately. In the latter half of thepaper, we derive the theoretical expressions for the bias errors in
the transfer functions estimated by the proposed and the standard
methods. By thoroughly comparing the resultant expressions,the superiority and the usefulness of the proposed method aretheoretically confirmed. Finally, the simulation experiments show
the advantages of the proposed method.

I. INTRODUCTION

T

HERE has been a dramatic increase in spectrum estima-tion research activities, especially in the past two decades,
since the digital fast Fourier transform (FFT) algorithm wasintroduced about 25 years ago [l]. The FFT has expanded therole of spectral estimation from research novelty to practicaluse. Although there are many disadvantages of such FFT-based
results, the FFT is most commonly applied to unknown signals
in many practical fields before using other parametric techniques because the expanded orthonormal functions employedin the FFT are predetermined and are independent of the signal.
Bingham et al. [4] discussed the computationally fastest
way to estimate the power spectrum of a time series fromthe FFT performed directly on the weighted data set. Welch[5] proposed a method for the application of the FFT tothe estimation of power spectra, which involves sectioning
the records, taking modified periodgrams of these sections,and averaging these modified periodgrams. The weighted
overlapped segment averaging is advocated by Nuttall andCarter [37], [41] to give stability and to minimize the impactof window sidelobes. Other references to the FFT and its
application for power spectrum estimation may be found in
Richards [2], Cochran et al. [3], Jenkins and Watt [6], Bertram
Manuscript received June IO, 1992; revised September 14, 1993. The
associate editor coordinating the review of this paper and approving it for
publication was Prof. Douglas Williams.
H. Kanai and N. Chubahi are with the Department of Electrical Engineering,
Faculty of Engineering, Tohoku University, Sendai, Japan.
T. Hori is with N7T,Japan
T. Ono is with Ono Sokki Company, Ltd., Tokyo, Japan.
IEEE Log Number 9401252.

[8], Glisson er al. [SI, Cooley er al. [lo], [ 111, Oppenheim andShafer [19], Yuen [29], and Amin [56].
In the power spectrum estimation, leakage effects arisingin the frequency domain due to the time domain windowing
can be reduced by the selection of windows with nonuniform
weighting. Bertram [ 121 provided a definitive descriptionof the leakage problem. Harris [30] and Nuttall [43] haveprovided a good summary of the merits of various windows. Nuttall and Carter [46] presented a spectrum estimation
method based on lag weighting. Other references include [16],[ X I , 1391, [4Ol, and W l .
As an application of these FFT-based spectrum estimation
methods, the coherence function has been developed for alinear measure of causality in the transfer function between
a pair of signals. Carter, Knapp, and Nuttall [17], [23] proposed the standard method to estimate the cross spectrum andthe coherence function by partitioning the two signals into
overlapping segments and computing the power spectrum ofeach segment via the FFT. The resultant power spectra arethen averaged to reduce the bias and variance of the resulting
estimates. The coherence function is effective in many fields
as pointed out by Carter and Knapp [20], and it has been
applied to system identification [17], measuring SNR andlinear-to-nonlinear power ratio [ 171, and determining signal
time delay [17], [18], [44]. Knapp and Carter [25] developed
a maximum likelihood estimator (MLE) for determining timedelay between signals based on the cross correlation, which
is identical to one proposed by Hannan and Thomson [15]. A
MLE for coherence is derived by Mohnkern [57]. Chan et al.
[38] proposed a regression approach of the coherence function
by solving the-discrete Wiener-Hopt equation. Youn er al. [45],[47] introduced an adaptive approach based on Widrow’s LMS
algorithm into the estimation of the coherence function for
nonstationary signals. A tutorial review of work in coherence
and time delay estimation was presented by Carter [50], [55].
Other references for the estimation of the coherence functionor the cross spectrum using the FFT and their applications arefound in Jenkins and Watt [6], Carter et al. [18], Talbot [22],Blake [28], Piersol [31], Seybert et al. [34], Barret [35], Chan[48], Cadzow [51], Cusani [58],and Mansour et al. [59]. The
error analyses for their estimates are provided by Benignus[7], Bendat [33], Walker [42], Schmidt [49], Mathews et al.
[52], and Gish et al. [54].
For the multiple input/output cases, Bendat [24], [27] provided methods to estimate the transfer function and the coherence function. Other references for the coherence function,
partial coherence, and system identification are found in Jenk-

1053-587X/94$04.00 0 1994 IEEEIEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. 7 , JULY 1994

1670

ins and Watt [6], Bendat and Piersol [13], [36], Dodd et al.
[21], and Romberg [32].
In the standard method established by these enormousstudies in literature, a transfer function is estimated by theratio of the averaged cross spectrum between the input andoutput to the averaged input power spectrum. If a time window
employed in the FFT is not long enough compared to thelength of the transfer system response, however, a bias error
due to the insufficiency of the time window length is generated.
The cause of the bias error is that the signal in the output
window does not contain the whole response to the signal
in the input window and contains an extraneous response to
the preceding input signal. As the system Q (quality factor)
increases, the bias error in the transfer function estimates getslarge especially around the resonant frequencies if the window
length is kept constant.
A method to determine the window length based on theallowable bias error limit was investigated by Ono et al.
[60]. A coherence function obtained from the signals in the
input window and delayed output window (called a delayed
block coherence function ) is proposed for the detection of a
missing signal in the output that is coherent to the input [60].
Based on the theoretical consideration for the delayed blockcoherence function, a method to estimate the damping factor ofa transfer system was presented in [61]. From the estimate ofthe damping factor, an optimum window length is determinedfrom the allowable bias error limit. It is, however, still difficult
to determine an accurate estimate of the transfer function witha long impulse response by these studies.
In this paper, after pointing out these problems in Section
11, we propose in Section 111 a new method to accurately
estimate the transfer function when the transfer system has along impulse response. By thoroughly comparing the accuracy
between the estimates obtained by the proposed method andthe standard cross spectrum method, we derive theoreticalexpressions for the estimates of the transfer function in Section
IV. From these expressions, the superiority and the usefulnessof the proposed method are theoretically and experimentallyconfirmed in Sections V and VI, respectively.
11. T H E STANDARD ESTIMATIONMETHODAND ITS PROBLEMS

Let us consider a single inputlsingle output system withextraneous noise at the output point only as shown in Fig. 1.
Let us assume that the input measurement r ( t )is essentiallynoise-free, whereas the output measurement y( t ) consists of) the noise
the sum of the ideal linear response z ( t ) to ~ ( tandcomponent n ( t ) .Using the discrete expression obtained bysampling each signal at an interval T , the resultant output
measurement is given byy ( n ) = z(n,)+ n(n)*+ n(n)

Fig. I .

Signal transmission system model employed in this paper.

and mutually uncorrelated. Let the average power of 34.) andn(71) be af and U:, respectively.
By applying an MN-point discrete Fourier transform to the
input signal z ( n ) and the output signal y(n), both of which
are windowed, where each window has MN points in length,(1) is written as+

Y ( k )= H ( k ) X ( k ) N ( k )(2)where k is an integer representing discrete frequencies, andthe discrete spectrum X ( k ) is given byh

The least squares estimate, which we denote by Hall(k),of the transfer function H ( k ) in (2) is obtained by the crossspectrum method [13], [36] such aswhere E [ . ]and * denote the ensemble average and a comp k x conjugate, respectively. For each frequency, the estimate
H,]l(k) minimizes the average power E[lN(k)12] of the noise
component n(n ) .
If the system response is 1ongLompared with the window
length, however, the estimate Hall(k) in (4) does not giveaEurate results. This is due to the following: The estimate
H,ll(lc) of (4) is obtained, assuming that the impulse responseto every input pulse in a block is dropped at the end of theblock as shown in Fig. 2(a). For example, the responses to theimpulses ~ ( o )x,( M N / 2 ) , and J ( M N - 1) are estimated sothat they have, respectively MN-point, MN/2-point, and 1point in length. These relations between ~ ( nand
) y(n) cannot
bzexpressed by a linear system. In the least square estimate ofHall(k) in (4), the nonlinear transfer system is approximatedby a linear system-As
a result, the bias error in the transfer
function estimate Hall( k ) becomes larger around the resonant
frequencies as the impulse response of the transfer system
becomes longer. This bias error is caused even for the caseof SNR = cc and for the case where the average number in(4) is infinite.
111. A NEW METHODOF ESTIMATING
THE TRANSFERFUNCTION

From a different point of view, let us consider the reason the= h,(71) 4 7 1 )(1) bias error happens in the transfer function, which is estimatedby the standard method as follows: Let us divide the input)
the output signal y(n) in a section with M N
where * indicates a convolution operation, and h ( n ) is an signal ~ ( 7 1 andimpulse response of the transfer system, which should be points in length into two blocks (x1(n) and 2 2 ( 7 1 ) ) and (yl(n)estimated below. Let us assume that the input signal ~ ( n )and y 2 ( 7 L ) ) , respectively, where each has M N / 2 points inand the measurement noise n ( 7 t ) are stationary white noise length as shown in Fig. 2(b). Using the spectra X l ( k ) , Y 2 ( k ) ,KANA1 er al.: DELAYED BLOCK TRANSFER FUNCTION IN THE FREQUENCY DOMAIN

M
N
,

I

1671

N POINTSPOINISI

w------I(a)

M N POINTS(b)Fig. 2. (a) Impulse response is dropped at the end of the block in the standard
method. (b) Illustration explaining the cause of the bias error in the estimatesobtained by the standard method.

Y l ( k ) , and Y2(k) of these block signals, the averaged cross] the numerator of (4) is described
spectrum E [ X * ( k ) Y ( k )inby+

E [ X * ( k ) Y ( k )=] E [ X ; ( k ) Y l ( k ) ] -qX;(k)Y2(k)l+(b)
Fig. 3. When the impulse response is long, the signal transmission systemmodel of single input/single output in Figs.1 and 3(a) is also described by themultiple inputlsingle output transfer system as shown in Fig. 3(b).

+(k)YI( k ) ] E[X;(k)Y2(k)l.

From the causality, the third term E [ X ; ( k ) Y l( k ) ] becomeszero. Let Ho ( k ) and H I ( k ) be the transfer functions from one
input block to the output block signal with the same timing
and the one-block delayed signal, respectively. Thus, the crossspectrum E [ X * ( k ) Y ( k )is] given by

E [ X * ( k ) Y ( k )=l E[IX(k)121{Ho(k)+ H o ( k ) + Hl(IC))# ~ [ I X ( k ) l 2 1 { H o (+k Hl(k)))
that is, bias error is caused in the cross spectrum estimated in(4). From this simple example, it is found that it is significant
to divide the input and output signals into short blocks andestimate the cross spectrum separately for each combination
between the input blocks and the output blocks and then sum
up the resultant transfer functions to obtain the total transfer
function. In this section, based on the princip&, a new method
is proposed to decrease the bias error in_H,n(k) of (4) andestimate an alternative transfer function H b k ( k ) ranging fromk = 0 to MN - 1 as described below.
Let us divide the MN-point length signal into M blocks,
where each of consists of N consecutive samples, as shown in
Fig. 3(a). Each block is identified by the block index. Whenthe system response is long, the transmission system in Fig. 1
is represented by a multiple input/output system as shown in
Fig. 3(b). Let an indexed H , ( k ) be the transfer function (called
the delayed block transfer function) from the input signal tothe i-block-delayed-output signal, and let hi ( n )be the impulse
response of H i ( k ) . The output signal ym(n) in the mth blockis the sum of the responses to each input signal in the preceding
blocks and the measurement noise signal n,(n) in mth block

i=O

where X,(k), Y,(k), and "(IC)are, respectively, the spec&(n), andtra of the zero-padded MN-length signals &(.),&(n) of the original N-length signals x,(n), y m ( n ) , andn,(n) in the mth block, where each is obtained as

x k ( n )={;T"(n), if n = 0, 1,.. . , N - 1;
if n = N , N 1,. . . , MN - 1. (7)+From Fig. 3(b), the spectrum Zm-i,m(k) of the response of
H i ( k ) to zm-i(n) is defined as follows:

zm ,m ( k) = Hi (k ) X m--i (k ) .
--i

(8)

The noise term N m ( k ) in (6) can be rewritten as
03

N m ( k ) = Y m ( k )- C ~ i ( k ) x m - i ( k ) .
i=O

The expectations of the product of N m ( k ) and its complex
conjugate N A ( I C ) give the average noise power PN( k ) at the1672IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. 7. JULY 1994

kth frequency such as

Using the relation ( l l ) , V ( k ) is

The term in the {.} brackets of this equation is equal to thespectrum K ( k ) , ( k = 0 , 1 , . . . , MN - l), which is obtainedby applying the MN-point FFT to the MN-point zero-padded
signal U:(.)
defined by

J=0
c c m

+

HZ

c=O

( k)H,*( k ) E [ X T l L - 2( k)xk-j ( k ) ]j=owhere P)v(k) is independent of the block index number
since the signal a ( n ) is assumed to be stationary. The leastsquare estimate of the transfer function is obtained as that) all possible choice of { H I( I C ) } .
which minimizes P N ( ~overBy setting the partial derivatives of Piy(k)with respect to{H,*(k)}equal to zerocxi

E [ X ; , - , ( ~ ) Y m ( k )=]H 1 ( l c ) E [ X , , ~ l ( k ) X ~ , _ , ( k .(9))l1=Owhere j = 0 , l . . ..,ca.Since it is assumed that the input
signals .cm ( n )and xm ( 7 1 ) are mutually uncorrelated for 712 #m', (9) reduces to
I

Using the spectrum b<(k),the total spectrum V ( k ) for the
MN-point signal u ( n ) is described by

In the relation of (13), by replacing K ( k ) and_V(k), respectively, by the delayed K c k transfer function H i ( k ) and thetotal transfer function H b k ( k ) , the estimate is given by

Hbl;( k ) =11-1

1 ( k )exp ( - j 2 r

2)
.

1=0

( k = 0 , 1 , . . . , M N - 1).

E [ X A _ ,(k)Y,(k)] = HJ(k)E[IX,-, ( k )12].

(14)

Thus

E [Xk- ( k )Ym 11
XTIL( k , 12]

Z(k)= E [I

L

(lo)

1

Thus, the estimate z l ; ( k ) of the total transfer function is
obtained from the M spectrum estimates { a ( k ) } of thedelayed block transfer function in (10).

where i = 0 , 1 , . . . , ca. This gives the least square estimate ofI v . THEORETICALDERIVATIONS
FOR
the delayed block transfer function from an input signal x m ( n )THE TRANSFERFUNCTION
ESTIMATES
in a block to the i-block delayed output signal ym+i(n).
To confirm the accuracy improEment in the estimze ofA method to obtain the estimate of the total transfer systemzk( k ) in (14), the estimates of H,ll(lC) in (4) and Hb,+(k)from the resultant delayed block transfer functions is described
below. To begin with, let us consider the relation between the in (14) and their bias errors are theoretically derived in thisspectrum K ( k ) of the ith block signal ui(n),which consists of section. To begin with, the characteristics of the true transfer
N points in length, and the spectrum V ( k )of the signal ~ ( n ) ,function are defined in Section IV-A. Then, after deriving
which consists of M blocks (= MN-point length), which are the theoretical expression for the expectations of the crosss p e c z m and p o x r spectrum in Section IV-B, the estimates
described as follows:
By applying the MN-point FFT to the M-block length of Hall(k) and H b k ( k ) are theoretically derived in Sections
IV-C and D, respectively.
signalW ( I I + N I ) = ~ ~ ( 7 1 ) .71 = 0. 1 . . . . . N - 1:

= O3'. . .

A. Dejnirion of the Characteristics of a Transfer Function

-

'

(l

the resultant spectrum V ( k ) is
AlN-1

v ( ~=c )Po(n)exp( -,jLT$)

n=O

1=1
1
=1

A-1

1N+N-1

1=0

n=zN

=

where { p l } are the poles of the transfer system ((pzI < l), and( C , } are the complex coefficients. By definingu(n)exp(-i2T&)
izI-1 A - 1(

w ( n + N i ) e x p -jar=
i=O

Let us assume that the transfer system is a rational transfer
function of order P and the impulse response h ( n )is described
by

n=O

k(n

+N%)

MN

1

1673

KANA1 et al.: DELAYED BLOCK TRANSFER FUNCTION IN THE FREQUENCY DOMAIN

the MN-point Fourier transform of the impulse response h(n)is obtained by

If the length of the FFT increases to the infinite, that is,
M -+ 00, the true transfer function H , ( k ) of H ( k ) in (17)is given byP

H,(k) =
i=l

ci- 1 - 1piz;.

B. Theoretical Derivations for the Cross Spectrum+

Let Fom(lc;q ) denote the kth spectrum of the signal h( nm N - q ) in the mth block of the response of the transfer
system to the impulse S(n - q ) at a time q in the 0th block
as shown in Fig. 4(a). For the case of m 2 1, by applying
the MN-point FFT to the zero-padded impulse responseh’(n m N - q )++m N - q)h(n + m N - q ) , if n = 0, I , . . . , N - 1;
if n = N , N + 1,.. . , M N - 1
= {O!

’

*0PTPPPd

ZERO-PADDED z’om(n)and substituting (15) into the resultant spectrum, the spectrum
Fo,(IC; q ) is given by1

MN-POINT FFT
ZOm(k)h’(n(b)
Fig. 4. (a) Fo, ( k ;q)denotes the MN-point spectrum of the zero-padded
signal h‘( R + m N - q ) in the mth block of the response of the transfer system
to the impulse 6( n - q ) at a time q in the 0th block; (b) 20, ( k ) denotes the
MiY-point spectrum of the N-point signal r o m ( n ) in the mth block of the
response of the transfer system to the input signal S O ( R ) in the 0th block.

spectrum Zom(IC)of zOm(n)is obtained by
N-1

Zom(k) =For the case of m = 0, by applying the MN-point FFT to thezero-padded ( N - q)-length impulse response h’(n - q )zo(q)Fom(kq )

(21)

q=o

where m = 0 , 1 , 2 , .. .. By substituting (19) and (20) into (21)the spectrum Foo(IC;q) is obtained byP

Foo(k4 ) =i=l

1 - (p&)N--q

cazkq 1 - p ; z t!

(for m = 0) (20)where the derivations of these spectrum Fom(lc;q ) andFoo(k;q ) in (19) and (20) are described in Appendix A.
Let zom(n) (m = 0 , 1 , . . .) denote the N-point signal in themth block of the response of the transfer system to the input
signal ~ ( nin )the 0th block as shown in Fig. 4(b). Using theexpression for the Fom(k;q ) of (19) and (20), the MN-point

Letting Gom(lc)be the expectations of the cross spectrum between X,-,( IC) and 2 0 , ( I C ) in (22) and (23), Gom( I C ) is given byIEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. 7, JULY 1994

1614

Gom(k) = E[X,*(k)Zom(IC)]P= Ne:@)- (pZz,")"c, 1 - p , z $,=l-( P Z ~ $ ) -m~N ~1 - (p,.,",-lpz

of the spectrum X, ( k ) anAK ( k ) of the resultant signals z, ( n )and y,(n), the estimate H,n(k) in (4) is described as follows:
hl-1*( m 2 1)

Hall(k) =M - r n G,(IC)M N2a2(IC)ni=O~.
exp( -jZr$).
(28)

(24)
The term Gm(IC)/N2c%(k)is equal to the z ( k ) in (10).
Thus, using the estimate % ( I C )of the d2ayed block transfer
function, the transfer function estimate Hall(k ) is obtained bywhere these derivations of the cross spectrum Go,(lc) andGoo(k) in (24) and (25) are described in Appendix B.
On the other hand, the expectations of the power spectrum
obtained by the MN-point FFT of the zero-padded input signal z ( n ) ,which consists originally of N points, is calculatedbyAAdescribed in Appendix C, by rearranging after substituting
H , ( k ) in (27), the cross spectra Go,(IC) in (24) and (25),and zo in (16) into (29)
P

=
i=l

c
z
- 1 pzz; [I+-. N 1M-

1 - (pzz,")"N
1 - (pzz,")-l

I(30)

E[IX(k)121

P

Cz*of H,ll(k) is equal toSince the first term
the true transfer functiof3_H, ( k ) in ( 18), the reEaining termshows the bias error AH,ll(k) in the estimate Hall(k).Thus,
the bias error A z ~ ( kis) given by
AZI(k)
Since E[z;(n)zo(q)] = a%6(n- y), the expectations of thepower spectrum is=

E&) H m ( k )
-

P
-

CL
~,

7=1

1 - p,z,"1
N(l - (p,z,")-l}1 - (pzZ,")*tN

M

'
(31)

n=O

= N2c:(k).

(26)

D. Transfer Function Estimate
Methodkz( k )for the Proposed

From the ratio of the average cross spectrum GoTn(k)in(24) or Goo(k) in (25) to the average power spectrum of z ( n )Let us consider the theoretical derivation for the estimatein (26), the estimatesin (10) for the delayed transfer H b k ( k ) of the total transfer system in (14) using the definition
function are obtained byof the impulse response h ( n ) in (15) as follows:
By substituting H;(IC) of (27) into (14) and using thedefinition of zo in (16), the estimate z k ( k ) is obtained byg(k){ G " I ( k ) + hl-1 Gom(k)$"}.

hC. Transfer Function Estimates zl( k )for the Standard
Cross Spectrum Method
As shown previously in (41, %,(IC)denotes the estimateof the transfer function obtained by applying the MN-point
FFT to the M-block signals ~ ( na n) 2 ( n ) , where each block
consists of N points. The estimate Hall(k) is rewritten by

Hbk(k) =~N2.E(k)m=l

By arranging after substituting the cross spectra Goo(k) of(24) and G o m ( k )of (25) into this equation as described inAppendix D, the estimate z k ( k ) of the total transfer system
obtained by the proposed method is given by

As described in Appendix C, the signals ~ ( n and) y(n)(32)
are divided into M-block signals z,(n) and ~ ~ ( 7 1(i) = Let us consider the physical meaning of z k ( k ) . Let
0 , 1 , 2 , . . . , A4 - l), respectively. Using the power spectrum be the transfer function estimate obtained by applying thelXi(k)12 and the cross spectrum Gl-,(k) = E[X:(k)Y,(k)] standard method to the N-point input and output signals. By

g(k)

1675

KANA1 et al.: DELAYED BLOCK TRANSFER FUNCTION IN THE FREQUENCY DOMAIN

STANDARD
CROSS
SPECTRUM
MJTHOD

becomes equal to the number of total points of the signals. In
this case, the spectrum Y,(k) of mth block, which consistsof one-point signal y(m), is given by

BLACKMANTURKEY
METLODPROPOSEDMKTHOD

H.l,(k)r---HHbL(k)---,HRT(k)THE NUMBER ( M )OF BLOCKS

,--

1

M

MN

N-1

1

N

+ m ) e x p (- j 2 r - t N )hY,(IC)=

MN
THE LENGTH ( N )OF EACH BLOCK

n=O

y(n

= Y(m).

Thus, the delayed block transfer function Hi(lc) in (10) is
obtained as*LENGTH (MN) OF
IMPULSE RESPONSE
IS KEPT CONSTANT.

(a)- E[.(m - ~ > * d m > lE[I.(m

-

4l21

h

I

I

I

= Rxy(i)I

MN
N
1
THE LENGTH ( N ) OF EACH BLOCK

where Rzy(i)denotes the estimates of the normalized correlation function between z ( n ) and y(n i ) . E o m (14), the totaltransfer function, which is denoted by H B T ( ~ )is, equal to
the MN-point Fourier transform of the correlation function
such as+(b)

Fig. 5. (a) Relation among the proposed method, the standard cross spectrum
m e g d , and the BlackmazTukey method; (b) theoretical ratio 9 ( k z :i , N )of
A H b k ( k ) in (33) to AHall(k) in (31) is always equal to or less than 1 forvarious lengths N of each block.

Kl(IC)ofsubstituting N into M N of the segment length in(30), the transfer function estimate g ( k ) is given by
Pg(k)= i=l ci 7
1 - p ; z o [1 + N 1 - (p;z;)-l .
1

1 1 - (p&N

1

The second t e r m k t h e right-hand side shows the bias error.
In the estimate Hall(k l o f (30), which is given when thewindow length N of Ho(k> is increased to M N , the bias
k M N
< 1, that is,error of g ( k ) is multiplied by
the bias error is decreasedhinverse proportion to the window
length. I k t h e estimate H b k ( k ) of (32), however, the bias
error of H o ( k ) is multiplied by exponentially decaying term
( p ; ~ t ) ( ~ - <<' ) ~1. Thus, a more accurate estimate of thetransfer function is obtained by the proposed method.
As the same manner in (31), the first term

&lT$;!;i)N

C
i
*

g
k

of
( k ) in (32) is equal to the true transfer function H , ( I C )in 2 8 ) , and the r e m a i n i s term expresses the bias error
AHbk(IC) in the estimate H b k ( k ) . Therefore
AH$-,@) = Z ( k ) - Hm(IC)
P

-

~. Ci
1 - piz;

5

1

N{1 -

x ( p i z ; ) ( M - l ) N (1 - ( P i W ) .

MN-1-

=

(33)

If the number M of the blocks and the length N of each blockin (33) are substituted by 1 and M N , respectively, as shown in
the extreme left of Fig. 5(a), the proposed method coincides
w i L t h e standard cross spectrum m e t h o d , z d in this case,
AHbk(lc)in (33) coincides exactly with AHall(IC) in (31).
Alternatively, when the number M of the blocks and thelength N of each block in (33) are substituted by M N and 1,
respectively, as shown in the extreme right of Fig. 5(a), eachblock consists of only one point, and the number of the blocks

i=O

R X y ( i ) e ~ p ( - jM
2 rNX )

(34)

which corresponds to the spectrum estimates obtained by theBlackman-Tukey(B-T) method [ 191. By substituting M x and1 into M and N of (33), respectively, the bias error A H B T ( ~ )in the estimate Hy~(k)is given by
AH?T(IC) = AH$-,(IC)M-MN,N-l(35)
V. ACCURACY COMPARISON BASEDON
THE THEORETICALLY DERIVED
EQUATIONS%comparing the d i f f e E c e of the bias error between
AHall(k) in (31) and AHbk(IC) in (33), it is clear that
there is difference between them o n e n the third terms(pi.,") ( M - 1 )N {lz(piz;)N}in AHbk(k) and &{l -pi^;)^^) in AHall(k).When the actual length of the signal
involved in the zero-paddeLsigna1 is e q 2 1 to N , let us define
the ratio ~ ( k ; i , Nof) AHbk(IC) to AHalI(k) for each orderi of the impulse response Cip; of the poles in the transfer
system as follows:
A% ( I C )defq ( k ; i , N )=~AZi(k)

1676

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. 7, JULY 1994

-2

(a)

0

15

31
47
63
TIME(P0INTS)

K -3

N=64, M = l

k? -4

20

y -3

N =16,M =4
0 2~/7x
2x
ANGULAR FREQUENCY(RADIAN)

Fig. 6. (a) Examples of the impulse response h ( n ) in (15) employed in
Figs. 7 and 8; (b), (c) frequency characteristics H=(I;) of the signal in
Fig. 6(a) defined by (18). The order P of the rational transfer function isequal to 1, C1 = 1.11111= 0.90, and L p l =The length >llAVis equal
to 64.

F.

The first term~l -~( pN, z( ;P) hlfzN" ) M,Nwhich

we denote by D ( k : i )hereafter, depends not on N but on the total length M N of theimpulse response, which is assumed to be a constant value.
Since the bias error is large especially around the resonant frequency, let us evaluate ~ ( ki ,;N ) at and near the
resonant frequency as follows: At the resonant frequency
IC; = M N L p ; / 2 x of the ith pole, the complex term pi.,"
becomes a real value, which we denote by n (0 < U < 1).
Thus, the ratio ~ ( ki :;N ) and D ( k ;i) become real at k = k;.
In this case, the ratio ~ ( k i ,:N ) at k = IC; is given by(37)
where D ( k ; i ) > 0. Since the partial derivation of a-"' =exp( -N In U ) with respect to N is equal to - u - ~I n a , thepartial derivation of q ( k ; ;i: N ) with respect to N is given by

Since the second term ( a N - lna" - 1) is positive for
0 < uN < 1, the gradientof q ( k Z 1: . N ) is positivefor 0 < a N < 1. Thus, r ) ( k z ;i , N ) increases monotonicallyas N becomes larger in the range of 1 5 N 5 M N as
illustrated in Fig. 5(b). When N is e q u L t o M N , that is,
M = 1. A H d , ( k ) becomes equal to AH,ll(k) as described
previously. In this case, ~ ( ki, ;N = M N ) = 1. Therefore, thebias error A K k ( k ) in the proposed method is always less thanor equal to A Z l ( k ) in the standard cross spectrum method
at the resonant frequency, that is, ,0 2x177r
2
AN G ULA R F R EQ,U,E,N,F,Y.,

-*,,

,

,RH"IHIY,

Fig. 7. Magnitude characteristics of the bias error in the estimate z k ( I ; ) by
the proposed method in (32) for different five combinations of the number 'ifof theblocks and the l e n g t h z o f a block. The dotted line shows the bias error
AHall(k)in the estimate H , l l ( k ) of the standard cross spectrum method in(30): (a) S = 64. A1 = 1; (b) S = 3 2 . .lf = 2; (c) X = 16. A4 = 4; (d)h' = S , d I = 8; (e) A: = ] . ' i f = 64.

Fig. 6 shows an example of the impulse response h(n) in(15) and its frequency characteristics H,(k) in (18) for thecase where the order P of the rational transfer function isequal to 1, C1 = 1, lpll = 0.90, and L p l =The length
M N is equal to 64.
Figs. 7 and 8 show, respectively, E m a g n i t u d e and phasecharacteristics of the bias error A H b k ( k ) of the proposed
method in ( 3 3 ) for five different combinations of M and N
under the condition that the total length M N isAways equal
to 64. The dotted line shows the bias error AH,il(k) of thestandard cross spectrum method in (31). The bias error in
H,ll(k) is large especially at and near the resonant frequency
kl =
x MN.
In Figs. 7(a) and 8(a), the bias error of x k ( k ) coincides
with that of z l ( k ) because the number M of the block isequal to 1 in these figures. As the number M of the blocks isincreased from top to bottom in Figs. 7 and 8, the estimation
accuracy is remarkably improved.
Figs. 9 and 10 show the characteristics of the bias errorA z k ( k ; ) as a function of the value of M , where the total
length M N is 64 and i p l =The magnitude in Fig. 9
and the phase in Fig. 10 of the bias error are represented by
the values at the discrete frequency kl =x M N ~9and k1 = 10, respectively. The dotted line shows the bias

q.

F.

__

KANA1 et al.: DELAYED BLOCK TRANSFER FUNCTION IN THE FREQUENCY DOMAIN

E
3

1611

0

i

t; -4 /

-4,

G. - 2 /
w
0 -

5
0

o1

""

,

,

,

,

,l(b)

-2
-3
-4

i

1 2
4 8 16 32 64
T H E NUMBER ( M ) OF BLOCKS
Fig. 9. Relation between the number hf ofblocks used in the proposed
method and the magnitude of the bias error A H b k ( k ) in (33) at the resonant
frequency k l =x M N N 9 under the condition that the total length
M B is always equal to 64.L p l =The dotted line shows the bias error in
the estimate H T l ( k ) of the standard cross spectrum method: (a) Ipl I = 0.95;(b) 1pi I = 0.9; (c) Ipi I = 0.8.

9-~/32, I,I

0 2 ~ 1 7
T
2lr
ANGULAR FREQUENCY(RADIAN)Kk

Fig. 8. Phase characteristics of the bias error in the estimate( k ) by the
proposed method in (32) for different five combinations of the number Af ofthexocks and the length Afz a block. The dotted line shows the bias error
AHall ( k )in the estimate Hall ( k )of the standard cross spectrum method in(30): (a) = 64.h1 = 1; (b) N = 3 2 . M = 2; (c) N = 16.,!21 = 4; (d)9 = 8. AT = 8; (e) A
Y = 1,211\ = 64.

error in the estimate %](IC)of the standard cross spectrum
method. For three different values of the magnitude of p l(Ipl I = 0.95, Ipl I = 0.9, and Ipl I = 0.5), the results are shownin Figs. 9(a) and 10(a), Figs. 9(b) and 10(b), and Figs. 9(c) and10(c), respectively. It is obvious that the bias error decreases
as the number M of the blocks becomes large. For the case of
M 2 2, the bias error of the proposed method is significantly
less than that of the s s d a r d cross spectrum method, and theresultant bias error AHbk(k) of the proposed method is smallenough for the accurate estimation of the transfer function.

VI. COMPUTERSIMULATIONEXPERIMENTS

-611~

F.

CR.0S.S....

i

SPECTRUM METHOD
K1 - t O

K op- -/

h

w

n

i
(I P I I = IO . ~ , ( C )16 32 64

-71.132,
i
I
1 2
4 8
T H E NUMBER ( M ) OF BLOCKS

In order to illustrate the advantage of the proposed method Fig. IO. Relation between the number A42f blocks used in the proposed
method and the phase of the bias error A H b k ( k ) in (33) near the resonant
to estimate the transfer function H ( k ) or its impulse response frequency El = 10 under the condition that the total signal length A4'V ish(n), we choose the example of the fourth order all-pole always equal to 64. Lp1 = %. The dotted line shows the bias error in thetransfer model H m ( k ) , where poles are 0.996exp(f27r . estimate H ? l ( k ) of the standard cross spectrum method: (a) lpll = 0.95;
65/360) and 0.997exp(f27r .80/360) as shown in Fig. 1 l(1- (b) Ipl( = 0.9; (c) lpil = 0.8.
b). The input signal ~ ( nin)(1) is assumed to be white noise,and the output signal y(n) is contaminated by measured white response h,(n) of H , ( k ) is shown for the same length innoise n ( n ) ,which is uncorrelated with the driving series ~ ( n ) .Fig. ll(1-a). Assuming that the true length of the impulse
The SNR equals 15 dB, and 32 768 points are generated for response is not known a a priori, the impulse response iseach of the signals ~ ( nand
) y(n). The generated two signals estimated for the length M N = 256 points in the following~ ( nand) y(n) are then divided into 128 sections, each of experiments for both the standard method and for the proposed
which has length ( M N )equal to 256 points. The true impulse method.

1678

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42. NO. I, JULY 1994

TIME n (POINTS)

Fig. I I . (1) Magnitude characteristics of the fourth-order all-pole transfer
model H,(k) employed in the computer simulation in Section VI and
its impulse response h , ( n ) ; (2) transfer function estimate H,ll(k) and
its impulse response h , l l ( k ) obtained by the standard method described in
Section 11.
I

Fig. 1 l(2) shows the transfer function estimates z l ( k ) andthe impulse response estimate h a l l ( n ) , which are obtained byapplying the standard method in Section I1 to the signals x;(,n)
and y(n,) in each segment with length MN = 256, which
is divided above. The number of the nonoverlapping average
operation E[.] in (4) is equal to 128 times. As described
Fig. 12. Impulse response estimates h , ( n ) , ( 1 = 0. 1 . . . . 8), which arepreviously in Section 11, the estimate Hall(k) is obtained so obtained by applying DFT to the block transfer function estimates g(k)that the impulse response to every input impulse in a segment multiplied by shift coefficients exp(-j2x . k z / . W ) . In the proposed method
is truncated at the end of the s e E e n t . In the magnitude in Section 111, each segment with length -IfAV = 126 is divided into
nonoverlapping 8 blocks, where each has 32 points in length (*U = 6.
of the transfer function estimate Hall(k) in Fig. 11(2-b), a 'V = 3 2 ) .
spectrum zero appears due& the truncation. The resultant
impulse response estimate h,l~(n)in Fig. 1 1(2-a) converges
more rapidly than the true characteristic h,(n) in Fig. 11( 1- the lengths !of the shortest and longest paths from the inputa). From these results, the true length of the impulse response block to the output block are 0 and-N - 1, respectively, as
shown in Fig. 13(a). However, for h l ( n ) in Fig. 12(b), theh,(n) cannot be estimated by the standard method.
!of the shortest and longest paths are 1 and 2 N - 1,
lengthsFigs. 12 and 14 show the results obtained by the method
respectively,
as shown in Fig. 13(b). Thus, the duration-time
proposed in Section 111. For the results in Figs. 12 and 14(1),each segment with length M N = 128 is divided into nonover- of the resultant impulse response estimate is 2 N - 1 for h l ( n ) ,lapping 8 blocks, each of which has 32 points in length, h ~ ( n .).,. , h8(n.) as shown in Figs. 12(b)-(i). For the estimatesthat is, in (7), (lo), and (14), the number A4 of the blocks hs(7t) in Fig. 12(i), the later half of the estimates is shifted toand the length N of each block are 8 and 32, respectively. the beginning of th5estimates due to the aliasing in the above
Fig 12 shows the impulse response estimates h,i(n).(i = DFT operation of H , ( k ) . e x p ( - j 2 ~. Ici/M).
Thus, by eliminating the beginning part of the estimate
0 , 1 , 2 : . . . , 8 ) , whickare obtained from the discrete Fouriertransform (DFT) of Hi(k).exp( -j27r.ki/M) in the right-hand h s ( n ) and summing up the resultant estimates K ( k ) .
side of (14). The impulse response & ( n ) shows the transfer exp(-j2r . k i / M ) for i = 0, 1?. . . , 8 in (14), the transfer
k(k
and
)
its
,impulse response h b k ( n )characteristics from the zero-padded input signal & ( n ) in the function estimates z
mth block with length N to the zero-padded i-block delayed is obtained as shown in Figs. 142-b) and (1-c), respectively.
output signal yL+;(n).For the estimate ho(71) in Fig. 12(a), The impulse response estimates h b k ( n . ) in Fig. 14(1-a) almost
h

~

h

A

h

KANA1 er al.: DELAYED BLOCK TRANSFER FUNCTION IN THE FREQUENCY DOMAIN+ho(n)

1679

MN is divided into MN blocks, where each has one point inlength. These results almost coincide with those in Fig. 14(1).
VII. CONCLUDINGREMARKS

A

+I

(71)

(c)

Fig. 13. Illustrations for explaining the relations between the length E of theestimated impulse response in Fig. I2 for each block transfer function and theamount of time for which the output block is delayed from the input block.

T I M E 17 (POINTS)*I2

FREQUENCY k(RADIAN)

This paper proposes a new method that accurately estimates
the transfer function of a transfer system with a long impulse
response. To begin with, a delayed block transfer function
was introduced to detect components that are correlated to thesignal in the input window but leaks from the output window.
The delayed block transfer function is defined as the transfer
function that is calculated by use of the input signal and theoutput signal in the a delayed output window. Next, from theresultant delayed block transfer functions, the total transfer
function is estimated, and the estimation accuracy is improvedby the proposed method.
We have also described the derivation of the thoroughtheoretical expressions for the estimates of the transfer functions obtained by the proposed and standard methods. By
comparing these expressions, the accuracy and the usefulnessof the proposed method was confirmed. From the computer
simulation experiments, the advantage of the proposed method
was also confirmed.
Two issues remain for further research as follows: In the theoretical derivations and the simulation experiments, the blocktransfer functions are computed over contiguous nonoverlupping blocks of the input and output signals, where each block
is cut off from the signals by using a a rectangular window. As
described in Section I, however, considerable work has been
done by Carter and Nuttall [37], [41] in applying Welch's
approach [5] to the estimation of the spectra and cross spectra
when the data blocks are overlapping. Moreover, it is also
necessary to determine the optimum window applied to eachblock when it is cut off from the original signal. These
important issues are currently under investigation. It is also
important to apply the proposed method to some practicalexamples.
APPENDIXA
Derivations for the Spectrum Fo,(k; q ) in Section IV-B

Fig. 14. (1) Transfer function estimate Hyk(/i)
and its impulse responseb G ( I ; )obtained by the proposed method described in Section 111 (,!f = 8 ,= 32); (2) transfer function estimate H B T ( I ; ) and its impulse response b y T ( I ; ) obtained by the proposed method described in Section IV-D(-21 = 128. AV = 1).

coincide with the true characteristics H,(k) in Fig. 1 I(1-a)
except for the amplitude Gk(0)for the first point, whichis affected by the contaminated noise components. F x thesame reason, the magnitude of the transfer function Hbk(k)in Fig. 14(1-b) has bias error, especially for the base part ofits resonant characteristics. These errors will decrease withan increase in the average number. From these figures, theestimates obtained by the standard method are improved by
the proposed method.
) (34) and its
Fig. 14(2) shows the estimates H y ~ ( ( k in
impulse response g k ( n ) , where each segment with length+

The kth spectrum Fom(k;(I) of the signal h(n m N - q )in the mth block of the response of the transfer system tothe impulse S(n - q ) at a time q in the 0th block as shownin Fig. 4(a) is obtained as follows: For the case of m 2 1,
by applying the MN-point FFT to the zero-padded impulse
response h'(n m N - q ) defined by

+

h'(n

+ mN

-

q)

h(n + m N - q ) ,= (0,if n = 0,1, . . . , N - 1;
if n = N , N 1,.. . , MN - 1+

the spectrum Fom(k;q ) is given by
N-l

Fom(k;q ) =h(n + m N - q ) z t n ,(for m 2 I) (A.1)

n=O

where zo = e x p ( - j 2 ~ / M N ) as defined in (16). For the caseof m = 0, by applying the MN-point FFT to the zero-padded

1680IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. 7, JULY 1994

( N - y)-length impulse response h'(n - y)whereh'(n - 4 ) =

h(7~
- y),

if n = y , q + 1,.. . , N - 1;otherwise

1, i f n = 0 ;0, otherwise.

lj(7L) =Letting each spectrum component of X ( k ) have the power
the spectrum F ~ o ( ky); is obtained by

0 2.;(IC)N-l

Foo(k;q ) =h (n - y)$.

(A.2)

at the kth frequency

By substituting h (n ) defined in (15) into (A.l) and (A.2), thespectra Fom(k;y) and F~o(k;q ) are, respectively, described byGom(k.1= XciP,

n=q

N-1

(B.2)

= -2

N

N-1

P

- (Piz$)" N a & ( k )1-pix;

i=lz,knnp;n

n=O

P

03.3)

n = O i=l
P

N-1

i=l

For the case of rn = 0, using (23), the cross spectrum Goo(k)is given byn=O

Goo(k)= E[X;(k)Zoo(k)l

i=l
P

2=1

n=O

1 - (p&)N-q
1 - paz;.

(A.4)

Equations (A.3) and (A.4) are used in (19) and (20), respectively, in the text.
APPENDIX B

Derivations f o r the Cross SpectrumGO,( k ) in Section IV-B
The expectations Go, ( k ) of the cross spectrum between the
spectrum X o ( k ) of the input signal in the 0th block and thespectrum Zo,(k) of the N-point signal in the mth block of
the response of the transfer system to the ~ ( nin )(22) and(23) is given by

Go,(k)= E[X;(k)Zom(k)lEquations (B.3) and (B.4) are used in (24) and (25), respectively, in the text.
APPENDIXC
Derivations for the Transfer Function Estimates zl( k )for
the Standard Cross Spectrum Method in Section IV-C
denotes the estimate of the transfer
As shown in (4), %,(IC)function obtained by applying the MN-point FFT to the M block signals ~ ( T L and)y ( ~ &where),
each block consists of N
points. The estimate z l ( k ) is rewritten by

M-1

(B.1)
for the case of m 2 1. Since the variance of the input signalx(n) is assumed to be
as described in Section 11,

02

E[z.(n)xo(Y)l = d 6 ( n - 4 )X ( k )=

Xi(k)exp(-j2a!3
i=O

M-1

Y ( k )=

x(k)enp(-j2n;)
i=O

(C.1)

16x1

KANA1 et ai.:DELAYED BLOCK TRANSFER FUNCTION IN THE FREQUENCY DOMAIN

and Xi(IC) and Y , ( k ) are spectra obtained by applying the
MN-point FFT to the zero-padded MN-point signals xi(.)
and y,!(n), respectively, as defined in (12). Since the signals{xi(.)}are mutually uncorrelated, the denominator of (4) is
calculated from (C.l) such as

The term G m ( k ) / N 2 a ~ ( kis) equal to the z ( k ) in (10).
Thus, using the estimateof the dzayed block transfer
function, the transfer function estimate H,ll(k) is obtained by

g(k)

M-1

M - m - H m ( k )exp( - j 2 7 r $ ) .

(C.8)

m=O
E[IX(k)121

A

By rearranging after substituting H m ( k ) in (27) and zo in

M-1

=E

[{

Xf(IC)ejZ.rrki/M

M-l

E[IXi(WI

=

M-1

1

h

Hall(k) =N2a2 ( I C )

i=O

M-m

~

"

Go, (

m=l

= MN2a2(IC)where the definition of E[1X(lc)12]is different from that ofE[lX(k)12]in (26).
On the other hand, the numerator of (4) is written by+M-l(1 - r N ) ( l- T - " )M - rn--rm"]
M

l.-T-l(c.9)

m=lwhere the cross spectra Gom(k)in (24) and Goo(k) in (25) arenxn =used, and T is defined by p&. Using the relationz(l-zK)
KZK+'
of (C.9)l-Z, the last term E,";is rearranged as follows:E,"=,

E [ X *(IC)Y(k)I

VT"

i=O

M-1
"T

m=l

Using the causality of the transfer system

E [ X f ( k ) K ( k )=] 0 f o r i > 1.

(C.4)

Since the input and output signals are assumed to be stationary,
the cross spectrum Gl-i(k) = E[X,*(IC)K(k)]between X i ( k )and % ( I C ) in the right-hand side of (C.3) is given by

Gi-i ( I C ) = E[X,*(k)K ( I C ) ]= E[X,*(IC)K-i( I C ) ] .

for 1 2 i.

(C.5)

Thus, the numerator of (4) is obtained by
M-1 M-1

E [ X * ( I C ) Y ( I=C ) I~~-;(k)exp

i=o l=i
M-1 M - 1 - i

Gm(IC)exp(-j27r2)

=
i=o

m=O

By substituting (C.2) and (C.3) into (4)
M-1
m=O

M-m

G m ( k ) exp ( q 2 4 - 3 . (C.7)

-

TN

M(1-T

)

{ ( M 1)-- MTN

+ T M N } . (C.10)

Using this result, the term in the square brackets [.] in (C.9)
is rewritten as follows:i1[ . ] = N + M(1 - r - l ) ( l - T N ) M ( 1 - T " ) 2

Since the term in the braces {.} of (C. 11) is equal to (1 r N ) ( l- r N M )the, term [.] in (C.9) or (C.11) is rearranged as[.] = N

(1 + M(1
r-1)

TNM)

-

(C.12)

1682

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. 7, JULY 1994

Thus, the theoretical expression of H,ll(k) is summarized asfollows:

the term in the square brackets [.] of (D.2) is rearranged asfollows:[ ‘ 1 = N ( l - ( p , z , “ ) - l } + 1 - (pz.,”)”+ (1 - ( p ~ z : ) y } ( p ~ z , “ ) N ( l- (p&)(M-l)N= N { 1 - (p,.,”)-’}+ (p,z$)(”-l)N

1(1 - (Paz0”)”I.

(D.4)
Thus, the estimate z k ( k ) of the total transfer system obtainedby the proposed method is given by

hf-1

m=O

P

=

a=1

1

[

I

1 - (p,z(y’”c,7
1 + -.
1- z ) , z ~NM I-(p,z$)-’

1. (C.13)This expression is used in (30) of the text.
APPENDIXDThis expression is used in (32) of the text.

Derivations of the Transfer Function Estimatezk ( k )f o r the Proposed Method in Section IV-D
In this Appendix, the theoretical derivation for the estimatez k ( k ) of the total transfer system in (14) is derived using thedefinition of the impulse responsLin (15) as follows:
By substituting the estimates H ; ( k ) of (27) into (14) andusing the definition of zo in (16), the estimate z k ( k ) is
obtained by

ACKNOWLEDGMENT
The authors would like to thank Dr. K. Kid0 of Chiba
Institute of Technology and Dr. H. Suzuki, Dr. K. Yamaguchi,
H. Imaizumi, and Dr. M. Ohashi of Ono Sokki Co., Ltd. fortheir effective discussion.

Substituting the cross spectra Goo(k) of (24) and Go,(k) of(25) into this equation

Since the last term ~equal to~~~(in the~ lasti equation
z
$
)is ~

REFERENCES
J. W. Cooley and J. W. Tukey, “An algorithm for the machine calculation of complex Fourier series,’’ Math. Computat., vol. 19, no. 90, pp.
297-301, 1965.
P. I. Richards, “Computing reliable power spectra,” IEEE Spectrum,
vol. 4, no. I , pp. 83-90, 1967.
W. T. Cochran et al., “What is the fast Fourier transform?,” IEEE Trans.
Audio Electroacoust., vol. AU-15 no. 2, pp. 45-55, 1967.
C. Bingham, M. D. Godfrey, and J. W. Turkey, “Modem techniquesof power spectrum estimation,” IEEE Trans. Audio Electroacoust., vol.
AU-15, no. 2, pp. 56-66, 1967.
P. D. Welch, “The use of fast Fourier transform for the estimation of
power spectra: A method based on time averaging over short, modified
periodgrams,” IEEE Trans. Audio Electroacoust., vol. AU-15, no. 2, pp.
70-73, 1967.
G. M. Jenkins and D. G. Watts, Specfral Analysis and Its Applications.
San Francisco: Holden-Day, 1968.
V. A. Benignus, “Estimation of the coherence spectrum and its confidence interval using the fast Fourier transform,” IEEE Trans. Audio
Electroacoust., vol. AU-17, no. 2, pp. 145-150, 1969.
S. Bertram, “On the derivation of the fast Fourier transform,” IEEE
Trans. Audio Electroacoust., vol. AU-18, no. I , pp. 55-58, 1970.
T. H. Glisson, C. I. Black, and A. P. Sage, “The digital computationof discrete spectra using the fast Fourier transform,” IEEE Trans. Audio
Electroacoust., vol. AU-18, no. 3 , pp. 271-287, 1970.
J. W. Cooley, P. A. W. Lewis, and P. D. Welch, “The fast Fouriertransform algorithm: Programming considerations in the calculation ofsine, cosine, and Laplace transforms,” J. Sound Vibrafion, vol. 12, no.
3 , pp. 315-337, 19iO.
-,“The applicationof the fast Fourier transform algorithm to the..
estimation of spectra and cross-spectra,’’ J. Sound Vibration, vol. 12,
no. 3 , pp. 339-352, 1970.
S. Bertram, “Frequency analysis using the discrete Fourier transform,”
IEEE Trans. Audio Electroacoust.,vol. AU-18, no. 4, pp. 495-500, 1970.
J. S. Bendat and A. G. Piersol, Random Data: Analysis and Measurement~Procedures. New York: Wiley, ,1971.
P. R. Roth, “Effective measurements using digital signal analysis,” IEEE
Spectrum, vol. 8, pp. 62-70, Apr. 1971.
E. J. Hannan and P. J. Thomson, “The estimation of coherence andgroup delay,” Biometrika, vol. 58, pp. 469431, 1971.
A. Eberhard, “An optimal discrete window for the calculation of power
spectra,” IEEE Trans. Audio Electroacoust., vol. AU-21, no. 1, pp.

I

KANA1 et al.: DELAYED BLOCK TRANSFER FUNCTION IN THE FREQUENCY DOMAIN

37-43, 1973.
G. C. Carter, C. H. Knapp, and A. H. Nuttall, “Estimation of themagnitude-squared coherence function via overlapped fast Fourier transform processing,” IEEE Trans. Audio Elecrroacoust., vol. AU-21, no. 4,
pp. 337-344, 1973.
G. C. Carter, A. H. Nuttall, and P. G. Cable, “The smoothed coherencefunction (SCOT),” Proc. IEEE, vol. 61, no. IO, pp. 149771498. 1973.
A. V. Oppenheim and R. W. Schafer, Digital Signal Processing. Englewood Cliffs, NJ: Rentice-Hall, 1975.
G. C. Carter and C. H. Knapp, “Coherence and its estimation via thepartitioned modified chirp z transform,” IEEE Trans. Acousr. Speech
Signal Processing, vol. ASSP-23, no. 3, pp. 257-264, 1975.
C. J. Dodds and J. D. Robson, “Partial coherence in multivariate random
processes,” J. Sound Vibration, vol. 42, no. 2, pp. 243-249, 1975.
C. R. S. Talbot, “Coherence function effects on phase difference
determination,” J. Sound Vibration, vol. 39, no. 3, pp. 345-358, 1975.
A. H. Nuttall, “Estimation of cross-spectra via overlapped fast Fouriertransform processing,” NUSC Tech. Rep., 4169-S, Naval Underwater
Syst. Cent., New London Laboratory, July 1 I , 1975.
J. S. Bendat, “Solutions for the multiple input/output problem,” J. Sound
Vibration, vol. 44, no. 3, pp. 311-325, 1976.
C. H. Knapp and G. C. Carter, “The generalized correlation method
for estimation of time delay,” IEEE Trans. Acoust. Speech Signal
Processing, vol. ASSP-24, no. 4, pp. 320-327, 1976.
H. Babic and G. C. Temes, “Optimum low-order windows for discrete Fourier transform systems,” IEEE Trans. Acousr. Speech Signal
Processing, vol. ASSP-24, no. 6, pp. 512-517, 1976.
J. S. Bendat, “System identification from multiple input/output data,” J .
Sound Vibration, vol. 49, no. 3, pp. 293-308, 1976.
W. K. Blake and R. V. Waterhouse, “The use of cross-spectral density
measurements in partially reverberant sound fields,” J. Sound Vibrarion,
vol. 54, no. 4, pp. 589-599, 1977.
C. K. Yuen, “A comparison of five methods for computing the power
spectrum of a random process using data segmentation,” Proc. IEEE,
vol. 65, no. 6, pp. 984-986, 1977.
1301 F. J. Harris, “On the use of windows for harmonic analysis with thediscrete Fourier transform,” Proc. IEEE, vol. 66, no. I , pp. 51-83, 1978.
~311A. G. Piersol, “Use of coherence and phase data between two receivers
in evaluation of noise environments,” J. Sound Vibration, vol. 56, no.
2, pp. 215-228, 1978.
1321 T. M. Romberg, “An algorithm for the multivariate spectral analysis oflinear system,” .I.
Sound Vibration. vol. 59, no. 3, pp. 395404, 1978.
[331 J. S. Bendat, “Statistical errors in measurement of coherence functions
and inputloutput quantities,” J. Sound Vibration, vol. 59, no. 3, pp.
405421, 1978.
A. F. Seybert and J. F. Hamilton, “Time delay bias errors in estimating
frequency response and coherence functions,” J. Sound Vibration, vol.
60, no. I , pp. 1-9, 1978.
1351 S. Barrett, “On the use of coherence functions to evaluate sources ofdynamic excitation,” Shock Vibrarion Bull., no. 49, pt. I , pp. 43-58,
1979.
[36] J. S. Bendat and A. G. Piersol, Engineering Applications of Correlationand Spectral Analysis. New York: Wiley, 1980.
1371 A. H. Nuttall and G. C. Carter, “A generalized framework for power
spectral estimation,” IEEE Trans. Acoust. Speech Signal Processing, vol.
ASSP-28, no. 3, pp. 334-335, 1980.
[38] Y. T. Chan, J. M. Riley, and J. B. Plant, “A Wiener filter approach tocoherence estimation,” in Proc. IEEE Inr. Con$ Acoust. Speech Signal
Processing (Denver, CO), Mar. 1980, pp. 646-649.
[39] N. R. Strader, 11, “Effects of subharmonic frequencies on DFT coefficients,” Proc. IEEE, vol. 68, no. 2, pp. 285-286, 1980.
[40] R. J. Webster, “Leakage regulation in the discrete Fourier transform
spectrum,” Proc. IEEE, vol. 68, no. IO, pp. 1339-1341, 1980.
[41] G. C. Carter and A. H. Nuttall, “On the weighted overlapped segment
averaging method for power spectral estimation,” Proc. IEEE, vol. 68,
no. IO, pp. 1352-1354, 1980.
(421 A. W. Walker, “The effect of bandwidth on the accuracy of transfer
function measurements of single degree of freedom response to random
excitation,” J. Sound Vibration, vol. 74, no. 2, pp. 251-263, 1981.
1431 A. H. Nuttall, “Some windows with very good side lobe behavior,”
IEEE Trans. Acousfics Speech Signal Processing, vol. ASSP-29, no. I ,pp. 84-91, 1981.
.
[44] K. Scarbrough, N. Ahmed, and G. C. Carter, “On the simulation of aclass of time delay estimation algorithms,” IEEE Trans. Acoust. Speech
Signal Processing, vol. ASSP-29, no. 3, pp. 534-540, 1981.
[45] D. H. Youn, N. Ahmed, and G. C. Carter, “On using the LMS
algorithm for time delay estimation,” IEEE Trans. Acoust. Speech Signal
Processing, vol. ASSP-30, no. 5, pp. 798-801, 1982.

1683

1461 A. H. Nuttall and G. C. Carter, “Spectral estimation using combined
time and lag weighting,” Proc. IEEE, vol. 70, no. 9, pp. 1115-1125,
1982.
[47] D.H. Youn, N.Ahmed, and G. C. Carter, “Magnitude squared coherencefunction estimation: An adaptive approach,” IEEE Trans. Acoust. Speech
Signal Processing, vol. ASSP-31, no. 1, pp. 137-142, 1983.
[48] Y. T. Chan and R. K. Miskowisz, “Estimation of coherence function andtime delay with ARMA models,” IEEE Trans. Acoust. Speech Signal
Processing, vol. ASSP-32, pp. 295-303, 1984.
[49] H. Schmidt, “Resolution bias errors in spectral density, frequency response and coherence function measurements I-VI,” J. Sound Vibration,
vol. 101, no. 3, pp. 347427, 1985.
[50] G. C. Carter, Toherence and time delay estimation,” in Proc. Ecole
d’CtC de physique fhkorique Traitement du Signal Processing (LesHouches, Haute Savoie, France), 1985, pp. 5 14-57 1.
[51] J. A. Cadzow and 0. M. Solomon, Jr., “Linear modeling and thecoherence function,” IEEE Trans. Acoust. Speech Signal Processing, vol.
ASSP-35, no. I , pp. 19-28, 1987.
[52] V. J. Mathews, D. H. Youn, and N. Ahmed, “A unified approach tononparametric spectrum estimation algorithms,” IEEE Trans. Acousr.
Speech Signal Processing, vol. ASSP-35, no. 3, pp. 228-349, 1987.
[53] C Y . Chi, D. Long, and F-K. Li, “Roundoff Noise analysis for digital
signal power processors using Welch’s power spectrum estimation,”
IEEE Trans. Acoust. Speech Signal Processing, vol. ASSP-35, no. 6,
pp. 784-795, 1987.
[54] H. Gish and D. Ccchran, “Invariance of the magnitude-squared
coherence estimate with respect to second-channel statistics,” IEEE
Trans. Acousr. Speech. Signal Processing, vol. ASSP-35, no. 12, pp.
1774-1779, 1987.
[55] G. C. Carter, “Coherence and time delay estimation,” Proc. IEEE, vol.
75, no. 2, pp. 236-255, 1987.
[56] M. G. Amin, “A new approach to recursive Fourier transform,” Proc.
IEEE, vol. 75, no. 11, pp. 1537-1538, 1987.
[57] G. L. Mohnkem, “Maximum likelihood estimation of magnitudesquared coherence,” IEEE Trans. Acoust. Speech Signal Processing,vol. 36, no. 1, pp. 130-132, 1988.
[58] R. Cusani, “Performance of fast time delay estimators,” IEEE Trans.
Acoust. Speech Signal Processing, vol. 37, no. 5, pp. 757-759, 1989.
[59] D. Mansour and B. H. Juang, “The short-time modified coherence
representation and noisy speech recognition,” IEEE Trans. Acoust.
Speech Signal Processing, vol. 37, no. 6, pp. 795-804, 1989.
[60] T. Ono, M. Ohashi, and H. Suzuki, “A new method for the determinationof an optimum window length in frequency response function estimation,” in Proc. Amer. Soc. Mech. Eng. Pressure Vessels Piping Con$,vol. 177, pp. 233-238, 1989.
[61] T. Ono, M. Ohashi, T. Maeda, and H. Kanai, “Estimation of thedamping factor of a signal transmission system by use of a delayed
block coherence function,” J. Acousr. Soc. Japan, vol. 46, no. 3, pp.
22CL228, 1990 (in Japanese).

Hiroshi Kanai (A’89-M’91) was bom in Matsumoto, Japan, on November 29, 1958 He received
the B E degree from Tohoku University, Sendai,
Japan in 1981 and the M E and the Dr Eng
degrees, also from Tohoku University, in 1983 and1986, both in electncal engineenng
From 1986 to 1988, he was with the Education
Center for Information Processing, Tohoku University, as a Research Associate From 1990 to 1992, he
was a Lecturer at the Department of Electrical Engineering, Faculty of Engineering, Tohoku University
Since 1992, he has been an Associate Professor at the Department of Electncal
Engineenng, Faculty of Engineering, Tohoku University His present interest
is in digital signal processing and its application on the acoustical, ultrasonic,
and electncal measurements
Dr Kanai is a member of the Acoustical Society of Japan, the Instituteof Electronics Information and Communication Engineering of Japan, theJapan Society of Mechanical Engineers, the Japan Society of Ultrasonics
in Medicine, the Japan Society of Medical Electronics and Biological Engineering, the Institute of Electrical Engineers of Japan, and the Japanese
Circulation SocietyIEEE TKANSACTIONS ON SIGNAL PROCESSING, VOL. 42, NO. I , JULY 1994

1684

Toyohiko Hori was born in Toyama, Japan on
March 5 , 1969. He graduated from Tohoku University, Sendai, Japan, in 1992.
Since 1992, he has been with "IT, Japan. Hisinterest is in digital signal processing and its application to acoustical measurements.

Noriyoshi Chubachi (M'83) was born in Kokura,
Japan, on October 5 , 1933. He received the B.S.,
M.S., and Ph.D. degrees in electrical engineeringfrom Tohoku University, Sendai, Japan, in 1956,
1962, and 1965, respectively.
In 1965, he joined the Research Institute of Electrical Communication, Tohoku University, where he
was an Associate Professor from 1966 to 1978.
Since 1979, he has been a Professor at the Department of Electrical Engineering, Tohoku University.
From 1982 to 1983, he was a visiting Professor of
Electrical and Computer Engineering, University of California, Santa Barbara.
He has worked on ultrasonic transducers and delay lines, surface acoustic
devices, acoustoelectronics, piezoelectric materials, acoustic microscopy, andrelated problems.
Dr. Chubachi is a ineinher of the Institute of Electronics and Communication Engineers of Japan, the Society of Japanese Applied Physics, theAcoustical Society of Japan, the Japan Society of Ultrasonics in Medicine,
the Japan Society for Nondestructive Inspection, the Japan Society of Medical
Electronics and Biological Engineering, the Japan Society of MechanicalEngineers, and the Institute of Electrical Engineers of Japan. He served aschairman, Tokyo chapter of IEEE UFFC Society from 1987 to 1988.

Takahiko Ono was born in Yokohama, Japan, in
September, 1951. He received the B.A. degree from
Keio University, Faculty of Business and Commerce, Tokyo, Japan, in 1974 In 1989, be received
the Dr Eng. Degree in electrical engineering fromTohoku University, Sendai, Japan.
In 1976, he joined Ono Sokki Co., Ltd. as the
product planning engineer and was named a memberof the board in 1983. Since 1991, he has been the
President and Chief Executive Officer of Ono Sokh
CO , Ltd Besides his business, he is now a part-timeteacher at Tokyo University of Agriculture and Technology
He is a member of the Japan Society of Mechanical Engineering and the
Acoustical Society of Japan