UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
Wavefront Detection By FrequencyWavenumber Analysis of ThreeDimensional Array Data By LEWIS J. PINSON A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1972 Dedicated to Iris because she is dedicated to me ACKNOWLEDGMENTS The author wishes to express sincere appreciation to Dr. D. G. Childers, chairman of his supervisory committee, for his guidance in selecting an interesting area of research and for his promotion of professionalism. He also wishes to thank Dr. R. L. Isaacson and Dr. L. W. Couch, the other members of his supervisory committee, and Dr. L. Jones and Dr. J. Cross for helpful suggestions. A special note of thanks goes to Dr. R. L. Isaacson who provided the financial support which made it possible for the author to return to the University of Florida. Finally, the author would like to express appreciation to the many people who offered encouragement to pursue a doctoral degree. A particular expression of gratitude is extended to his mother and fatherinlaw, Mr. and Mrs. C. B. Mann, because their encouragement was the first. TABLE OF CONTENTS ACKNOWLEDGMENTS............................................................. LIST OF TABLES.................................................. LIST OF FIGURES................................................ ABSTRACT........................................................ CHAPTERS: I INTRODUCTION. .............................. ... II ARRAY PROCESSING..................................... A. Adaptive Processors for SignaltoNoise Ratio Improvement.................................... .. B. Feature Extraction by Digital Filtering........... C. Wavefront Detection.............................. III FREQUENCYWAVENUMBER SPECTRUMTHEORY................. A. Generalized Plane Waves.......................... B. FrequencyWavenumber Representation of Spatial Fields.................................. C. Estimate of the FrequencyWavenumber Spectrum..... D. High Resolution Adaptive Window................... E. ThreeDimensional FrequencyWavenumber Spectrum... F. Finite Bandpass Signals........................... IV FREQUENCYWAVENUMBER SPECTRUMPARAMETRIC ANALYSIS.... A. Comparison of High Resolution and Conventional Estimates ........... ...... ...... .. ..... ....... B. Array Geometry .................................. Page iii vii viii xi 1 4 4 6 7 8 9 12 13 18 22 25 29 30 32 CHAPTERS: C. Temporal Window Functions........................ D. Additive Incoherent Noise....................... E. Statistical Reliability........................ F. Spherical Wavefronts............................ G. Spatial Prefiltering............................ H. Cost and Computation Requirements................ V FREQUENCYWAVENUMBER SPECTRUMSIMULATION ........... A. Array and Signal Example Parameters.............. B. Comparison of High Resolution and Conventional Estimates........... ................... ... . C. Array Geometry .................................. D. Temporal Window Functions...................... E. Additive Incoherent Noise........................ F. Bandpass Signal Spectral Computation............. G. ThreeDimensional Simulation..................... VI INVESTIGATION OF ELECTROPHYSIOLOGICAL DATA BY FREQUENCYWAVENUMBER TECHNIQUE...................... A. Electrophysiological Signal Characteristics...... B. Analysis of Human VisualEvoked Response......... C. PenicillinInduced Focal Epilepsy in Rats....... VII DISCUSSION..... ........... ... ................. APPENDICES: A B C CROSS POWER SPECTRAL DENSITY......................... 107 SPECTRAL REPRESENTATION OF SPATIOTEMPORAL FIELDS.... 111 FREQUENCYWAVENUMBER SPECTRUM IN SPHERICAL COORDINATES ......................................... 116 34 36 39 42 43 45 48 49 51 53 57 63 68 76 81 81 83 90 103 D EXPERIMENTAL METHOD.................................... 118 E PROGRAM LISTINGS........................... .......... 124 F SINGULARITY OF THE CROSS POWER SPECTRAL DENSITY MATRIX. .............. ... .... .... ................. 163 BIBLIOGRAPHY................. ............ .... ...... ................ 165 BIOGRAPHICAL SKETCH............................................... 170 LIST OF TABLES Table Page 1 Cross Power Spectral Density Estimate...................... 12 2 Computation Time and Cost................................. 46 LIST OF FIGURES Figure Page 1 GENERAL TWODIMENSIONAL PLANE WAVE ..................... 11 2 CONVENTIONAL FREQUENCYWAVENUMBER SPECTRAL ESTIMATE...... 17 3 EFFECT OF FINITE DATA RECORD LENGTH ON RESOLUTION....... 19 4 HIGH RESOLUTION FREQUENCYWAVENUMBER SPECTRAL ESTIMATE.. 21 5 THREEDIMENSIONAL PLANE WAVE............................. 23 6 BANDPASS SPECTRAL ESTIMATION........................... 28 7 SAMPLING THEOREM RELATIONSHIPS.......................... 33 8 TEMPORAL WINDOW FUNCTIONS AND CHARACTERISTICS.......... 35 9 TWODIMENSIONAL ARRAY GEOMETRY AND PARAMETERS ........... 50 10 CONVENTIONAL AND HIGH RESOLUTION PERSPECTIVE PLOTS...... 52 11 CONVENTIONAL ESTIMATE OF THE FREQUENCYWAVENUMBER SPECTRUM .............. ............ ....... .............. 54 12 HIGH RESOLUTION ESTIMATE OF THE FREQUENCYWAVENUMBER SPECTRUM ................................................. 55 13 CROSSSECTION OF CONVENTIONAL AND HIGH RESOLUTION SPECTRAL ESTIMATES ...................................... 56 14 RELATIVE EFFECT OF ARRAY APERTURE ON RESOLUTION......... 56 15 HIGH RESOLUTION ATTENUATION PLOT FOR 4ELECTRODE ARRAY.. 58 16 EXAMPLE TEMPORAL WINDOW FUNCTIONS....................... 59 17 EFFECT OF RECTANGULAR WINDOW FUNCTION FOR DISCRETE FREQUENCY COMPUTATION AND R = 0.1....................... 60 18 EFFECT OF COSINE TAPERED WINDOW FUNCTION FOR DISCRETE FREQUENCY COMPUTATION AND R = 0.1.............. ........ 61 viii Figure Page 19 EFFECT OF RECTANGULAR WINDOW FUNCTION FOR DISCRETE FREQUENCY COMPUTATION AND R = 0.0001.................... 62 20 EFFECT OF COSINE TAPERED WINDOW FUNCTION FOR DISCRETE FREQUENCY COMPUTATION AND R 0.0001 ..................... 64 21 EFFECT OF INCOHERENT NOISE ON SINGLE WAVEFRONT........... 65 22 EFFECT OF INCOHERENT NOISE ON MULTIPLE WAVEFRONTS........ 65 23 SINGLEWAVEFRONT SPECTRUM FOR R = 0.001.................. 67 24 TWOWAVEFRONT SPECTRUM FOR R = 0.01 (COSINE WINDOW)...... 69 25 TWOWAVEFRONT SPECTRUM FOR R = 0.001 (COSINE WINDOW)..... 70 26 BANDPASS SPECTRUM (R = 0.01, RECTANGULAR WINDOW, PASSBAND = 620 Hz)................ ................. 71 27 BANDPASS SPECTRUM (R = 0.01, COSINE WINDOW, PASSBAND = 620 Hz)........ ................ ................... 72 28 NARROW BANDPASS SPECTRUM (R = 0.01, RECTANGULAR WINDOW, PASSBAND = 11.516 Hz)................................ 74 29 NARROW BANDPASS SPECTRUM (R = 0.01, COSINE WINDOW, PASSBAND = 11.516 Hz).................................. 75 30 THREEDIMENSIONAL ARRAY GEOMETRY AND SIGNAL PARAMETERS... 77 31 THREEDIMENSIONAL FREQUENCYWAVENUMBER SPECTRUM.......... 79 32 FREQUENCYWAVENUMBER SPECTRUM OF VER DATA (0.5 sec)...... 85 33 FREQUENCYWAVENUMBER SPECTRUM OF FIRST HALF VER DATA (00.25 sec)..... ... .................... .. ....... ... 86 34 FREQUENCYWAVENUMBER SPECTRUM OF LAST HALF VER DATA (0.250.50 sec).......................................... 87 35 FREQUENCYWAVENUMBER SPECTRUM WITH MAXIMUM VALUE BANDS... 89 36 ERRAT1 SPECTRUM (f1 = 2.7; f2 = 4.7 Hz; PHASE COM PENSATED) ... ........ ............................ ..... 93 37 ERRAT2 SPECTRUM (2 mm grid; fl = 3; f2 = 4; f3 = 5 Hz)... 95 Figure Page 38 ERRAT2 SPECTRUM (1 mm grid; fl = 3; f2 = 4; f3 = 5 Hz)... 96 39 ERRAT3 SPECTRUM FOR PRIMARY HEMISPHERE (f = 2 Hz)........ 98 40 ERRAT3 SPECTRUM FOR SECONDARY HEMISPHERE (f = 2 Hz)...... 99 41 ERRAT3 SPECTRUM FOR PRIMARY HEMISPHERE (fl = 7.8; f2 = 13.7 Hz)............ ....... ...................... 100 42 ERRAT3 SPECTRUM FOR SECONDARY HEMISPHERE (fl = 9.8; f2 = 21.5 Hz)................. ........................... 101 D1 SCHEMATIC DIAGRAM FOR MONITORING AND RECORDING ELECTRICAL ACTIVITY OF RAT NEOCORTEX ..................... 119 D2 PREPARATION OF RAT DATA FOR FREQUENCYWAVENUMBER ANALYSIS.......................... ...................*.... 122 Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy WAVEFRONT DETECTION BY FREQUENCYWAVENUMBER ANALYSIS OF THREEDIMENSIONAL ARRAY DATA By Lewis J. Pinson March, 1972 Chairman: Dr. D. G. Childers Major Department: Electrical Engineering A high resolution frequencywavenumber spectral estimation technique is employed to analyze spatiotemporal signals from three dimensional sensor arrays, to estimate the spectrum for finite temporal frequency bands using FFT techniques, and to estimate the spectrum for discrete temporal frequencies. Direction and true vector velocity of a plane wavefront propa gating across a threedimensional arrayare given by the frequency wavenumber spectrum for discrete temporal frequencies. For finite temporal frequency bands the wavefront direction and bounds on the vector velocity are given. A parametric analysis is performed in which the properties of the frequencywavenumber spectral estimate are determined for known signals and parametric values. Array aperture and the ratio of incoherent noise power to total signal plus noise power are found to determine the resolution in the wavenumber domain. Other parameters, including number of temporal samples per sensor, temporal window function, and number of sensors, are found to affect statistical reliability and computation cost. The frequencywavenumber spectral estimation technique is applied to the analysis of visualevoked response (VER) data and to the analysis of penicillininduced epileptiform electrical activity recorded from the rat neocortex. Results indicate that application of the frequencywavenumber spectral analysis method to electro physiological data is possible and that it provides information for increased understanding of brain function. Application of the expanded frequencywavenumber spectral estimation method to spatiotemporal data from the areas of seismology, oceanography, speech, meteorology, radar and sonar would require in most cases only a change in parametric inputs. xii CHAPTER I INTRODUCTION Over the past few years Biomedical Engineering has evolved from isolated applications of electronic instruments as tools for medicine to sophisticated and everexpanding applications of advanced analytical techniques to biological systems. Other areas of applica tion, which include seismology, oceanography, speech, meteorology and sonar, have need for analytical techniques similar to those applicable to the analysis of biological data. These various disciplines have certain common characteristics which makes it advantageous to develop analytical techniques that have general applicability. One such characteristic is that each of the areas is concerned with measuring and interpreting spatiotemporal signals. The purpose of this research was to investigate a technique for analyzing spatiotemporal signals and develop useful extensions of the technique to provide additional information about the signals. Parametric bounds and example computations were chosen to emphasize application of the technique to the study of electrophysiological data. However, application to other areas is implied and in most cases would require only a change in parametric inputs. Chapter II gives a brief review of arrayprocessing techniques which have been applied to other fields such as radar, seismology, acoustics, and image processing. Some of the properties of these techniques are discussed and compared with the expected properties of the frequencywavenumber technique. Digital filtering array processors are indicated as having general applicability as wavenumber limiting prefilters. In Chapter III the theoretical development is given for the frequencywavenumber spectral estimate. It is shown how the frequency wavenumber spectrum provides information for determining direction and velocity for multiple plane wavefronts. An extension of the frequency wavenumber technique to three dimensions is developed and shown to be valid theoretically. Inclusion of frequency bandpass signals as well as discrete frequency signals has been accomplished using fast Fourier transform (FFT) techniques and a spectral integration scheme. Details of the bandpass feature are given in Chapter III. Many parameters affect the resolution, cost and reliability of the frequencywavenumber spectral estimate. Additive incoherent noise was found to have a major effect on resolution; however, the array aperture and the temporal window were found to interact with the incoherent noise level and thus affected the achievable resolution. Parametric variation, statistical reliability and computation require ments are discussed in Chapter IV. For spherically symmetric wavefronts the frequencywavenumber method gives a scalar output which is insufficient for determining wavefront direction. Chapter IV discusses the case for spherical wavefronts as well as the need for spatial prefiltering. In Chapter V examples are shown for simulation runs which verify the theoretical development and parametric analyses. The frequencywavenumber technique was applied to electro physiological data of two types: 1. human visualevoked response data recorded from the scalp over the occipital region of the brain, and 2. penicillininduced focal epileptic discharges recorded from the rat neocortex. Details and results of these experiments are outlined in Chapter VI. Finally in Chapter VII, a discussion of the frequencywavenumber spectral estimation technique is given in relation to the goals accom plished by this research effort. Significance of the experimental results is discussed along with recommended areas for further investiga tions. CHAPTER II ARRAY PROCESSING Included under the category of array processing are several specific techniques which have been investigated by various researchers. The characteristics of some of these techniques are discussed in rela tion to the arrayprocessing technique termed frequencywavenumber spectral estimation. In general, arrayprocessing techniques are used in applications where their advantages outweigh the disadvantage of increased processing complexity. The advantages offered by arrayprocessing techniques include improved signaltonoise ratio for facilitating the detection of weak signals, enhancement or emphasis of specific features or parameters and the determination of spatial characteristics such as wavefront direc tion. Major differences between the various arrayprocessing techniques usually result from: 1. noise assumptions (stationarity, directionality); 2. signal assumptions (known, unknown, uniform or distorted across the array); and 3. method employed (spatial filtering, stochastic approxima tion, steepest descent, conjugate gradients, noise estimation, etc.). Various methods of arrayprocessing are discussed below. A. Adaptive Processors for SignaltoNoise Ratio Improvement Although many techniques have been developed for improving signaltonoise ratio using array data, most of the techniques are similar in approach and make the assumption that wavefront direction is known "a priori." Lacoss [1] has presented an adaptive linear processor with variable coefficients in which the coefficients are adjusted by a rule similar to that used by the projection gradient method of a quadratic form subject to a linear constraint. The purpose of this processor is to improve signal detection and assumes that the wavefront direction is known with 'sufficient' accuracy. Thus the signals from the sensors can be combined after appropriate time delay adjustments with gains adjusted to provide optimum detection of an unknown signal. An optimum processor is one which converges to a minimum variance, unbiased estimate of the signal. Lacoss derived several iterative methods for reducing sensitivity of the processor to data anomalies such as noise bursts and for reducing required computation time and storage requirements. Other adaptive processing techniques for improving signaltonoise ratio which are similar to those of Lacoss have been derived in the literature. Shor [2] has proposed a gradient adaptive method for narrowband hydrophone arrays which requires knowledge of the signals' autocorrelation function. Adams [3] has done work with adaptive online array processors for receiving deepspace probe signals. Kobayashi [4] has proposed two iterative methods for processing array data: the method of steepest descent and the method of conjugate gradients with projection. These methods are very similar to those of Lacoss and rely on essentially the same mathematical development. No intermediate computation of covariance matrix functions is required for the methods used by Kobayashi. As do most methods, however, this method assumes the only difference between the signals at various sensors is a time delay which can be compensated. Burg [5] used a multidimensional Wiener filtering approach to array processing. In this approach the signal and noise are represented as stationary random processes with known crosscorrelation functions. The output of the processor is a minimum meansquarederror estimate of the signal. Again, knowledge of time delays between sensors is required. Another approach to array processing is that proposed by Claerbout [6]. The crosscorrelation matrix is computed by observing the sensor outputs in a "fitting interval." This information is then used to provide a minimum meansquarederror estimate of the noise for a short, future, projected interval. The predicted noise is then subtracted from the signal. This method reduces the noise level but introduces distortion into the signal. Capon et al. [7] have discussed the advantages and disadvantages of timedomain versus frequencydomain array processors. The signalto noise improvement is best with timedomain methods; however, timedomain methods are more sensitive to nonstationary noise and signal anomalies than the frequencydomain methods. B. Feature Extraction by Digital Filtering The fields of feature extraction and digital filtering are by far too extensive to be given a comprehensive survey in this document; however, feature extraction from twodimensional fields by digital filtering techniques is a form of array processing and is included for that reason. In addition to feature extraction, digital filters may be used in array processors to limit wavenumber content and reduce aliasing effects for a given spatial sampling rate. Wavenumber limited spatial filters are discussed further in Chapter IV. C. Wavefront Detection To determine a wavefront propagation direction it is necessary to use methods which give vector wavenumbers as outputs. Standard beamforming techniques use adjustable phase differences between elements of a sensor array to determine wavefront directions. These methods, however, provide no information regarding wave velocities and are restricted in resolution to the natural beam pattern of the array. A method which overcomes these deficiencies and computes the high resolution vector velocity and direction of propagating wavefronts has been proposed by Capon [8] for processing twodimensional data from the largeaperture seismic array (LASA) in Montana. A comparison of this high resolution estimate of the frequencywavenumber spectrum with conventional estimates showed marked improvement for the high resolution method. The theoretical development of the frequency wavenumber spectrum is given in Chapter III. Extension of Capon's method to include data from threedimensional arrays is shown to be valid theoretically. Methods for including finite frequency bands are discussed and one method is developed in detail. CHAPTER III FREQUENCYWAVENUMBER SPECTRUMTHEORY This chapter gives a detailed development of the theoretical background material supporting the validity of the frequencywavenumber technique as a descriptor of array sensed signals. First the properties of plane waves and their mathematical representations are discussed and it is shown how the frequencywavenumber spectrum provides informa tion describing plane waves. The concept of estimation as related to the frequencywavenumber spectrum is introduced and related to theoretical concepts. Window functions are discussed in relation to their effect on accuracy of the estimate. A high resolution adaptive window is introduced which provides a major improvement in resolution over that obtainable with conventional fixed window functions. The use of a twodimensional array for measuring threedimensional signals is inadequate if it is desired to know true vector velocities of incoming wavefronts. Determination of wavenumbers in two dimensions identifies azimuth angle but not elevation angle. This means that for a given azimuth angle, any number of signals with different elevation angles may be present. With the twodimensional analysis these various signals would be indistinguishable and only their phase velocities could be deter mined. A threedimensional analysis identifies wavenumbers in three dimensions and uniquely determines wavefront directions. This also makes it possible to determine true vector velocity. Theory and methods are discussed for extending the frequency wavenumber technique to accept data from a threedimensional array of sensors and to accept finite bandwidth data in addition to discrete frequencies. Capon's technique is valid for only single discrete frequencies. This could be a major disadvantage where the exact spectral content of the signal is unknown. Also in many applications it may be desirable to include a band of frequencies even if spectral content is known. Rather than compute spectral estimates individually for all frequencies of interest a method was developed for accepting multiple frequencies as well as finite continuous frequency bands and efficiently computing the frequencywavenumber spectrum. A. Generalized Plane Waves A general plane wave in two dimensions can be represented mathematically by [9],[10] S(t,x,y) = exp [2iri(ft i *f)] (1) where: S(t,x,y) = the spatiotemporal signal amplitude f = temporal frequency (Hz) k = vector wavenumber (cycles/m) g ? = spatial coordinate vector (m) For nonisotropic media, the vector wavenumber (k) is spatially dependent and is a function of x and y. For an attenuating medium k has both real and imaginary components. Thus the general vector wave number is given by k = Re[k (x,y)] + Im[k (x,y)] (2) g g g where: Re['] denotes the "real part of" Im[*] denotes the "imaginary part of" If the signal is represented in an isotropic, nonattenuating medium then k is a real, vector constant and is denoted by k. The g spatiotemporal signal representation for a twodimensional, isotropic, nonattenuating medium then becomes S(t,x,y) = exp [2Ti(ft kir)] (3) Figure 1 shows the above relationships and how the vector wave number specifies the propagation direction for the plane wavefront. For simplicity and clarity in demonstrating the frequencywavenumber technique an isotropic, nonattenuating medium will be assumed so that for the following analysis equation (3) is representative of a plane wave. Some of the properties of a plane wave or any general spatio temporal signal which is to be measured by an array of sensors are often more easily studied by measuring the cross power spectrum between any two sensors. The cross power spectral density is a K x K matrix for an array of K sensors and is related to the coherency matrix which provides a measure of the linear dependence among the sensors. The cross power spectral density matrix may be obtained in either of two ways: the correlation method or the direct method. Details of these two methods are given in Appendix A. Table 1 provides a comparison of pertinent criteria for the correlation method versus the direct method for estimating the cross power spectral density. Because of reduced computation time and less susceptibility to signal variations the direct method was chosen for the estimate. Statistical reliability is improved by segmenting the data as described by Welch [11]. In some cases however, increased / 1/ / / Ik12 = 2 + 2 E = vector wavenumber (gives wavefront direction) GENERAL: S(t,x,y) k = exp [2ii(ft k *9)] = Re[k (x,y)] + Im[kg(x,y)] 2 2 = x +y ISOTROPIC, NONATTENUATING MEDIA: S(t,x,y) = exp [2wi(ft KZ)] = exp [2ri(ft Ox yy)] WHERE: k = Re[k (x,y)] = constant S= component of in x direction y = component of k in y direction y = component of k in y direction Figure 1. GENERAL TWODIMENSIONAL PLANE WAVE statistical reliability is achieved at the cost of reduced resolving capability. Table 1: CROSS POWER SPECTRAL DENSITY ESTIMATE CORRELATION METHOD DIRECT METHOD S 2 p(k) eikX 1 k=L/2 (k) ) m=1 F F() L = NO. OF DATA POINTS Fjm(X) IS FOURIER TRANSFORM DEFINITION j (k) IS COVARIANCE OF DATA IN mth SEGMENT OF MATRIX jth SENSOR. M SEGMENTS OF N DATA POINTS EACH FOR EACH SENSOR COMPUTATIONALESS MORE LESS TIME REQUIRED STATISTICAL GOOD FOR LARGE L GOOD FOR LARGE M & N RELIABILITY SUSCEPTIBILITY TO SIGNAL VARIA MORE LESS TION B. FrequencyWavenumber Representation of Spatial Fields In Appendix A it is shown how the cross power spectral density function for K sensors is a K x K matrix. The cross power spectral density matrix is a temporal frequency transform of the cross covariance matrix for an array and its measured signals. The cross power spectral density matrix contains information concerning the frequency content of the signal power. In a similar manner it is possible to obtain information about the wavenumber content of the signal by performing a spatial frequency transform on the cross power spectral density function. The result is a frequencywavenumber spectrum and is defined as P(X,k) = J f(A,x) e2i d5 (4) 00 Vector where f(X,k) is the cross power spectral density function and the integral is a vector Fourier transform with respect to the spatial vector, 1. It follows that the cross power spectral density is then the Fourier inverse of the frequencywavenumber spectrum. 00 f(X,x) P(X,k) e dk (5) 00 Vector This development is based on a treatment by Yaglom [12] and details of the development are given in Appendix B. It is also shown in Appendix B that if the signal waveform is a discrete frequency plane wave then the frequencywavenumber spectrum is a multidimensional delta function about the given frequency and vector wavenumber, i.e., P(X,k) = (X o0) 6( ko) (6) Thus it is shown that wavefront detection is possible through the use of the frequencywavenumber spectrum. C. Estimate of the FrequencyWavenumber Spectrum An exact spectral representation for spatiotemporal signals would require infinite summations over time and space. Since this is not practical a spectral representation using a finite number of sensors (spatial representation) and a finite number of time samples (temporal representation) is computed and distinguished from the exact representa tion by being called an estimate of the exact spectrum. Criteria for measuring goodness of an estimate have been developed and are discussed in Chapter IV, Section E on statistical reliability. As was shown previously there are two major methods for estimating power spectral densities, the correlation method and the direct method. Because of its desirable statistical characteristics and computational efficiency the direct segment method or block averaging technique was used for estimating the cross power spectral density. The resulting matrix is denoted by Oj (X)] to distinguish it as an estimate of the exact spectrum f j(X). This estimate of the cross power spectral density is then used to obtain an estimate of the frequencywavenumber spectrum. To illustrate the technique assume that discrete data are avail able from K sensors. Data from each of the sensors are divided into M nonoverlapping segments of N data points each to give L = MN data points per sensor. The cross power spectral density matrix is computed using the direct segment method. The Fourier transform of the data in the mth segment of the jth sensor is F () a S (nT) ein (7) N n=ln j = 1,,K m = 1,,M where: S M(nT) X = 2ffT a n = the data in the mth segment from the jth sensor = normalized frequency = the weighting coefficients describing the temporal window used in estimating the frequency spectrum. For simplicity assume a = 1; n = 1,,N. In practice, however, side lobe leakage can be reduced and the estimate improved in some cases by using windows different from the rectangular window described here. From (7) an estimate for the cross power spectral density is a K x K matrix whose elements are given by M 1 F(X) Fm (N) m=1 = In order to remove the effect of differences in gain and other causes of improper sensor equalization, a normalization is performed on f i(t) A f 1/ f j& ) =  1/2 IA f^ i ) M fzt )]1 The cross power spectral density matrix is closely related to coherency, the major difference being that phase information is retained in the normalized cross power spectrum and lost in the coherency function. The relationship is given by Y 2(X) = I2 (10) where y2 (A) is the coherency function. The coherency function varies it from zero to one and is a measure of linear dependence between sensor elements. Using the normalized cross power spectral density estimate and the coordinates of the sensors, we can obtain the conventional estimate for the frequencywavenumber spectrum: K K ik2*(5t 5t) ,) = IKl wjw* e (11) K j=l A=l1 where the wj are the weights describing the shape of the spatial window used to estimate the wavenumber spectrum. Again for simplicity we assume w = 1; j = 1,,K so that wj represents the aperture of the sensor array. As was indicated for the temporal window, the spatial window can be designed to provide an improved estimate of the spectrum. One such window, proposed by Capon [8], is a function of the particular wavenumber being considered. It is different for each wavenumber k o and passes undistorted any monochromatic plane wave traveling with a velocity corresponding to the wavenumber k It also suppresses in an optimum least squares sense the power of waves traveling at veloc ities corresponding to wavenumbers other than k Figure 2 shows a flow chart describing the estimation procedure for obtaining the frequencywavenumber spectrum. It provides a summary of the procedure discussed up to this point and the resultant output is called the conventional frequencywavenumber spectrum of a homogeneous random field. SENSOR ARRAY K SENSORS SPATIAL WINDOW FUNCTION SPATIAL TRANSFORM P (A,k) Figure 2. CONVENTIONAL FREQUENCYWAVENUMBER SPECTRAL ESTIMATE MULTICHANNEL TEMPORAL SAMPLING FUNCTION t + f f t f MULTICHANNEL SEGMENTING CROSS POWER SPECTRAL DENSITY 1 F F* M jm tAm M(A) fm=l fji() ) =,, 1/2 D. High Resolution Adaptive Window A finite data record of discrete values may be considered as the product of an infinite representation of the data with a rectangular sampling window as shown in Figure 3a. From Fourier transform theory we know that multiplication in the time domain implies convolution in the frequency domain. The result of the convolution is shown in Figure 3b for the time functions in 3a. An ideal situation would be to have a window function whose transform is an impulse. This of course implies a window function of infinite width which implies an infinite data record length. Thus it is desirable to find a window function which more closely approximates the ideal but with a finite data record. More detail on window functions and their characteristics which have been used in this research is given under Section C of Chapter IV on the parametric effects of temporal window design. The choice of which window to use depends on the particular application and criteria selected to determine a desirable output. There is no window function which is generally accepted as optimum. Capon introduces a high resolution estimate for the frequency wavenumber spectrum which utilizes a window function that is dependent on the frequency and wavenumber being considered. The high resolution frequencywavenumber spectral estimate using this window is given by K K i27TT(5 ) P'(,kI) = A(X,k) A (A,k) fJC(A) ek (12) j=l a=l1 where: AJ(A,;) is the high resolution window function. Comparison of equations (11) and (12) illustrates the fact that the j) I'  I 10 0 c * . >0a iyi r~ O z Iia o z I Ii o Ii I 2 H Fz4 W B . W C,, F1 ca t 0c 00 _ _I  0 !tl z H H C/3 f No SE1 E ~ high resolution estimate is the direct result of choosing an adaptive window function. The window function is given by K Aj (,k) =K K (13) I I q j (X,l) j=l =l where: [qj (Xk)] = [[exp {27rit;(i 5)}'[fj( l]]1 (14) that is, [q j] is the inverse of the matrix resulting from the product of the exponential matrix [exp (2lrik*(~j i))] with the cross power spectral density matrix [f j]. Further, if [q ] is defined as the inverse of [fj] [ j] = [ j]1 (15) then it can be shown that an equivalent expression for the high resolu tion frequencywavenumber spectral estimate is given by K K P'(A,k) = ,q(X) exp [27rik(< (6)] (16) j=1 i=1 Thus P' can be computed by simply taking the Hermitian inverse of the cross power spectral density matrix and then computing the spatial transform. Equivalency of the two expressions, (12) and (16), is difficult to show in general; however, computations with low order matrices have verified their validity. The estimate P' is the power output of a maximum likelihood array processor whose characteristics vary with wavenumber in an MATRIX INVERSION [qj = [fj]l HIGH RESOLUTION ESTIMATE P'(A,k) = 1/S.T. Figure 4. HIGH RESOLUTION FREQUENCYWAVENUMBER SPECTRAL ESTIMATE optimum least squares sense. That is, P' provides a minimum variance, unbiased estimate of the power. For each wavenumber, k it passes undistorted any monochromatic plane wave with a velocity corresponding to the wavenumber k and suppresses in an optimum least squares sense all waves traveling at velocities other than ko ,Figure 4 gives a flow chart outlining the procedure for obtaining the high resolution frequencywavenumber spectral estimate. The major difference between it and the conventional technique is the inclusion of the matrix inversion step. Since [?it] is in general complex, additional care must be exercised in implementing the technique. E. ThreeDimensional FrequencyWavenumber Spectrum The theoretical development of the frequencywavenumber spectrum is not limited to a specific number of dimensions. Therefore the equations developed for the twodimensional frequencywavenumber spectrum are directly applicable to detection of threedimensional wavefronts if the vectors k and x are extended to include components in three dimensions. The detection of wavefronts in three dimensions would be particularly useful for the measurement of seismological events via earth body waves and the measurement of EEG signals which can certainly be considered to possess threedimensional characteristics. Although the theoretical development is straightforward, difficulties arise in the implementation of a threedimensional technique. A pictorial representation of three dimensional wavefront characteristics is given below. Then some of the implementation difficulties are discussed further. Figure 5 shows how a plane wave might be represented in three dimensions. Determination of the three wavenumber components kx, ky , z k k / A a " y mA Azimuth Angle x aE Elevation Angle Vector Wavenumber (Direction of Propagation) ISOTROPIC, PLANE WAVE IN NONATTENUATING THREEDIMENSIONAL MEDIA: S(t, ) = exp [2rri(ft *.)] S(t,x,y,z) = exp [2ri(ft k x k y k z)] x y z FREQUENCYWAVENUMBER SPECTRUM: P'(A,k) = I qij(A) exp [2wik*(k, .)] Lj=l X=1 ] WHERE: [q.Z(A)] = inverse of the cross power spectral density matrix Figure 5. THREEDIMENSIONAL PLANE WAVE and k specifies not only the azimuth and phase velocity of the signal z but provides the added advantage of specifying the elevation angle and group velocity of the signal. Determination of the wavenumber com ponents is accomplished through the estimation of the threedimensional frequencywavenumber spectrum. The major problems involved in the extension of the frequency wavenumber spectrum to three dimensions are complexity of implementation, increased cost, how to display the data in a meaningful way and how to design sensor arrays to best sense data in specific applications such as biological preparations. An order of magnitude increase in computa tion time and thus cost results from computation of the spatial transform in three dimensions over that required for computation in two dimensions. Fast Fourier transform (FFT) techniques are not readily applicable to the spatial transform operation. Since most FFT algorithms require that the number of data points be a power of two and since all FFT techniques are based on data points at equal sampling intervals, the use of these techniques places stringent requirements on the sensor arrays to be used. Equal spacing of electrodes for the measurement of biological signals may in some cases be impossible or even undesirable. If this is so then a more general technique for determining the frequencywavenumber spectrum is required, which precludes use of the fast Fourier transform. Extension of the frequencywavenumber technique to three dimensions was accomplished and test cases were run using both known input signals and data recorded from the rat neocortex using a threedimensional electrode array. Details of the results of the threedimensional frequencywavenumber analysis are discussed in Chapter V for the known input signal. Meaningful display methods for the threedimensional spectrum depend to some extent on the wavefront complexity. A method which provides twodimensional spectra for sequential values of the third dimension was used for all threedimensional data runs because of its ease of comparison with the twodimensional data runs. This method worked quite well for simple, known input signals; however, for more complicated signals it proved to be less clear for demonstrating wavefront directions. The particular array geometry to be used for the three dimensional spectral estimate is strongly dependent on the specific application. However, to demonstrate many of the principles and parametric effects, a simple cubic array geometry with known input signals was used. For practical application of the technique to remotely sensed spatiotemporal signals such as the EEG, certain assumptions about the characteristics of the signal were made. These assumptions are discussed in Chapter VI along with experimental results. F. Finite Bandpass Signals Recall from the equations defining the conventional and high resolution frequencywavenumber spectral estimates that both are func tions of X the normalized frequency. This says that for each frequency contained in the signal a different wavenumber spectrum can be computed. Clearly this procedure is not practical for the case where exact frequency content is not known or where a large range of frequency components is known to compose the signal. The question then arises as to how the technique may be modified to accept a narrow finite band of frequencies. Consider the case of a composite signal given by s = i s (t,xy,z) (17) J where: s (t,x,y,z) = exp [2?i(f t k x k y kz z)] (18) In order for the technique to identify the components of the signal, it must compute the cross power spectral density estimate which includes all the frequencies, fj. For a small number of frequencies this may be accomplished by computing each component and then adding to get the total. This particular technique was tested for the case of two frequency components and appeared to introduce very little distortion into the frequencywavenumber spectral estimate. Up to five discrete frequencies have been included in the spectral estimation technique with only a minor increase in computation time. For large numbers of discrete components or where interest in finite bandpass signals exists an integration technique on the output from an FFT algorithm can be used. The inclusion of multiple frequencies introduces some distortion into the estimation of the frequencywavenumber spectrum, as does the finite bandpass technique. The extent of this distortion is determined by a complex interaction of parameters including incoherent noise ratio, data length, and temporal window function. A method for including finite bandpass data was developed which uses the average integrated spectral output over the selected spectral band. The IBM subroutine, HARM, was used in the fast Fourier transform spectral computations. If the output of HARM is as illustrated in Figure 6 where the spectrum shown is for the jth sensor, mth data segment, then an estimate for the bandpass spectral content is HIGH Fjm(A) = Ai*Af (19) i=ILOW where: Ai = spectral value at ith frequency ILOW =FLOW Af FHIGH HIGH = F H+ 1 FLOW = low frequency cuton FHIGH = high frequency cutoff 1 Af = frequency resolution = record length (sec) The Fjm(A) are then used in the direct segment method for computing cross power spectral density. It is important to note that the inclusion of multiple frequencies either as discrete components or finite bands introduces ambiguity into determination of phase velocities or group velocities. However, if wavefront directions are of prime importance then this method of accepting multiple frequencies will provide more information more efficiently. H 41 HI H44 II II H 44 0 CO Cx CHAPTER IV FREQUENCYWAVENUMBER SPECTRUMPARAMETRIC ANALYSIS In general any spectral analysis technique should produce an estimate that is statistically reliable and which has high resolution. It is of interest therefore to identify parameters affecting the estimate and to quantify their effects on resolution, statistical reliability, and other measures of goodness of the estimate. The spatial and temporal sample intervals, Ax and AT, are determined primarily by the expected bandwidths in their respective transform domains. These sample intervals are chosen to reduce aliasing in the transform domain of the periodically repeated spectrum. Temporal resolution is shown to be inversely proportional to record length and spatial resolution is inversely proportional to array aperture. Therefore for fixed sample intervals, Ax and AT, resolution is improved by increasing the number of sample values. The effect of an increased number of sample values on cost is partic ularly important when discrete signals in time and three spatial coordinates are considered such as those for the threedimensional frequencywavenumber spectral estimate. Another parameter which results from choosing the direct segment method for estimating the cross power spectral density matrix is the number of segments M into which each data record is divided. The value of M has an effect on the statistical reliability of the frequency wavenumber spectral estimate and should be large, which implies more data points and greater cost so a tradeoff is indicated. Temporal window design can be utilized to reduce sidelobe leakage and in some cases improve resolution. Various window func tions and their characteristics are discussed. Included under the parametric analysis are discussions of spherical wavefronts in relation to the objective of determining wavefront direction and spatial prefilters in relation to expected maximum wavenumbers. A. Comparison of. High Resolution and Conventional Estimates It was stated that Capon's high resolution spatial window function improves resolution over that obtainable with a conventional window function. It is therefore of interest to give a quantitative comparison of the conventional and high resolution estimates in rela tion to achievable resolution. The following analysis is derived largely from Capon's paper [8]. Consider a single plane wave (wavenumber = k ) propagating across an array of sensors plus a component of noise in each sensor which is incoherent between any pair of sensors. Assume that M, the number of data segments per sensor, and N, the number of data points per segment, are large so that an unbiased, consistent estimate of the frequencywavenumber spectrum is obtained. Then the cross power spectral density matrix elements are estimated as f ,(Xo) = 6j(R) exp [i2'Tok(x xi)] j,A = 1,,K (20) where: A = 2wf T o o f = temporal frequency of incident wave (Hz) lo = vector wavenumber of incident wave (cycles/m) x = vector position coordinates of jth sensor and 6 z(R) = 1 ; J = (21) = lR; j A R incoherent noise power (22) total power of sensor output Using the above expression for the input signal and the equations for conventional and high resolution estimates of the frequencywavenumber spectrum it can be shown that P( =,o) = (1R)IB(Ako)2 + R/K (23) and p( ) R 1R + (R/K) (24) S1R + 2(R/K) Po(c) where: Ar = k (25) o o and K B(k) = exp [27ri*rc ] (26) S=1 and K equals the number of sensors. For K = o P(ZE) = P' (o,Eo) = lR + R/K (27) In the neighborhood of k consider those values of k such that the conventional estimate gives a value close to its peak value of P(Xk). Since R is small (0 : R 1) and K is usually large (about 16) the neighborhood for the conventional estimate can be quite large. For these values of k, the high resolution estimate P'(A ,k) becomes p .(o R 1R + (R/K) P'(A k) = (28) =o K IR + 2(R/K) (lR) or 1 P'(Ao,k) = (1R + R/K) (29) which is already 3 db down from its peak value. Thus P' has a much steeper falloff and gives a higher resolution estimate for the frequencywavenumber spectrum. B. Array Geometry The data to be represented by their frequencywavenumber spectrum consist of finite groups of spatiotemporal samples. That is, spatial sampling is a result of a finite number of discrete sensors located within a finite aperture. Temporal sampling is the result of discrete sampling of a finite time record from each sensor. The relationship between these parameters and the frequencywavenumber spectrum is given by an application of the sampling theorem, modified to include multi dimensional signals. For example the transform pair represented in Figure 7 demonstrates properties which hold for any bandlimited, truncated, discrete signal which includes temporal and spatial signals. IT 2B ^1/T j< 1/AT REQUIRE: AT S 1/2B RESOLUTION IN $ PLANE = 1/T Figure 7. SAMPLING THEOREM RELATIONSHIPS In general, the spectral resolution varies inversely with the record length and the required sample interval varies inversely with the bandwidth of the signal. Sampling theorem concepts can be applied to determine the effects of array geometry on the resolution of the frequencywavenumber spectrum. Application of these concepts is, however, complicated by the additional property of dimensionality for spatial fields. For example, equal sample intervals for a spatial field has meaning only if referenced to some direction in the field. For a threedimensional field, equal sampling intervals can be achieved in directions represented by three perpendicular coordinate axes by constructing a cubic sensor array. However, given no physical limitations on the array size or shape, Petersen and Middleton [13] have shown that the optimum array (i.e., minimum number of sample points per unit hypervolume) is not in general rectangular. The optimum array is represented by the centers of hyperspheres in their closest packed arrangement. For two dimensions the optimum array is rhombic and for three dimensions the tightest packing of spheres is achieved with a hexagonalclose packed configuration. These arrangements not only provide maximum sampling rates for a given number of electrodes but also provide equal sample intervals between any two adjacent sample points. For an equal number of elements and a given spatial sampling rate, arrays other than these have a reduced array aperture and resolution is not as good. For any particular array design the sampling theorem relation ship still holds that resolution in a given wavenumber component is inversely proportional to the aperture of the array in that dimension. In most applications the array design is determined by physical limitations of the recording area. This is particularly true in the case of human EEG recordings in which the array geometry is determined by scalp contour. The rationale for array geometries used in the experimental data examples is given in Section G and in Chapter VI. C. Temporal Window Functions It was shown in Chapter III that the spectrum of a truncated signal is given by the convolution of the signal transform with the transform of the window function. An ideal temporal window function would require an infinite length of time; therefore, it is of interest to investigate various finite temporal window functions relative to 04 04 rI rf H P o HFZ( 00 rz FH4 E > d E4 try M nI C. 8 0 (u' vJ 00 i0 + z c8i U2 8p E j their sidelobe levels, bandwidths, sharpness of cutoff, and statistical reliability. The simplest window function is a rectangle which has a sin x transform of the form In general no window function is optimum. Rather, the best window design for a particular application depends on the requirements for resolution, statistical reliability and sidelobe leakage. Several window functions and their characteristics are shown in Figure 8. For the frequencywavenumber analyses two window functions were used. The rectangular window function was chosen for its simplicity and as a basis for comparison. The cosine tapered window function was also used. Because of its reduced sidelobe level the cosine window is expected to provide a more reliable estimate than the rectangular window. For the frequencywavenumber spectral estimate, the particular temporal window used had little effect on spectral resolution in comparison with other parameters such as additive incoherent noise level and array aperture. Effects of the various parameters on resolution are demon strated in Chapter V. D. Additive Incoherent Noise Since the high resolution estimate requires a matrix inversion of the cross power spectral density matrix [f], it is of importance to determine under what conditions this inverse exists. From matrix theory it is known that an Hermitian matrix has an inverse if the matrix is positive definite [15]. Capon et al. [7] have shown that the cross power spectral density matrix [f] when computed by the block averaging method is nonnegative definite. However, this is not sufficient to guarantee that an inverse exists. In fact, if the number of data blocks M is less than the number of sensors K then the matrix is singular. This comes from the fact that [f] is the sum of M matrices of rank unity. M [2] = [q ] (30) m=l where: [q ] is a 1 x K row matrix whose jth element is given by N inX mJ n=~ jn+N(m1 ) e (31) Since the rank of [f] cannot exceed the sum of the ranks of its constituent parts, then [f] is a K order matrix with rank M. Clearly if M < K the matrix is singular. Thus a necessary but still not sufficient condition for [2] to have an inverse is that M a K. To insure that an inverse exists the matrix [r] is modified to give the matrix f' which can be shown to be positive definite. [f'] = (1 R)[f] + R[I] (32) where: R = ratio of incoherent noise power to total power at a sensor and: [I] is the identity matrix Thus [f'] is a K x K matrix with elements given by f' = (1R) ? +R; m n mn mm m n mm (33) = (lR) fm ; m n Now [f'] is positive definite if for any function aj, the quadratic form Q given by K K Q = a Bafja (34) j=1 =l is positive. Using the definition for [f] and [f'], we have 1R M K N inX K L im Il J ajSjn+N(m.) e + R I jajl j=1 7r S X 5 7T (35) The first term above is always nonnegative for it represents a positive constant (lR) multiplied by the quadratic form for [2] which has previously been demonstrated to be nonnegative definite. Thus the only way for Q to be zero is if K 2 = Slaj12 0 (36) j=l But this implies that aj = 0; j = 1,,K which is a contradiction. Therefore [f'] is positive definite and its inverse exists. The quantitative dependence of the high resolution frequency wavenumber estimate on the magnitude of additive incoherent noise is derived from expressions given in Section A and is given by pA ) R 1R + (R/K) p o' K +K 2Wi() (37) 1K R + (R/K) (1R) I e o j Recall that the peak value for P'(o ,k=k o) is P'(o,(0) = 1R + (R/K) (38) so that P'(Ao' ) R/K _____ 0 _R/K (39) P'(A ,k ~ *K 2Ti(ko 2 (R/K) + (1R) 1 1 e e j=l Effects on resolution of the additive incoherent noise are demonstrated in Chapter V by artificially injecting noise into a known signal. ( See also Appendix F.) E. Statistical Reliability Two desirable properties of any estimator are that it be un biased and consistent. For the random process P it is then required of the estimate P that E[PN] = P (unbiased) (40) lim E[(P FN)2] = 0 (consistent) NKx That is, the mean value of the estimate should equal the mean value of the random process and as the number of observations goes to infinity, the mean square error between the estimate and the random process should go to zero. Capon and Goodman [16] show that the conventional and high resolu tion frequencywavenumber spectral estimates are asymptotically unbiased, that is, in the limit as N (number of data points per segment) gets large the bias approaches zero. The expected value for the conven tional estimate is given by E[P(X )] = f f P(x,) B(i k 12 (A)2 dk (41) 7 o0 Vector K i27trk where: B(k) =i e (42) j=1 is the beamforming response pattern of the array and, 2 WN(x)I2 1 sin (N/2)x (43) sin x is the Bartlett window resulting from the rectangular time window chosen for the temporal transform. Thus P will be an unbiased estimator for P if in the limit the frequencywavenumber window approaches a delta function. f J IWN(xX)22 B(kko I E dk 6(xXo) 6(k) (44) Vector Following Goodman [17] and Capon and Goodman [16], it can be shown using the Complex Wishart distribution that the high resolution estimate is unbiased if it is multiplied by M/(MK+1). That is the expected value of P' is = MK+ P(x,) IB,(X,r,,zo)I 2 IWN(x_X)12 dx d9 7r Vector (45) K where: B'(,k,k ) = 1 A.(A,iK ) exp [27ri(K )* ] (46) j=l is the high resolution beamforming pattern of the array and is a func tion of Z Since in many cases all spectral values are normalized, the question of bias is not so important. Also for M >> K the high resolution estimate gives the same expected value as the conventional estimate. Assuming that the noise received by the array is a multi dimensional Gaussian process, it can be shown for large N that the variance of the conventional estimate is given by VAR [P(A,k)] = {E[P(A, )]}2 ; k = o (47) E[)1 ]}2 Since the variance goes to zero as M goes to infinity the estimate is consistent. Again using the Complex Wishart distribution it can be shown that the high resolution estimate is a multiple of a chisquare variable with 2(MK+1) degrees of freedom and variance given by VAR [P'(A,k)] = M K EP'+(X,)]2 (48) Thus both the conventional and high resolution estimates for the frequencywavenumber spectrum are consistent and can be made to be unbiased. Details of the Complex Wishart distribution and derivations of the means and variances for the estimates are given in [16], [17]. F. Spherical Wavefronts The development of the frequencywavenumber spectrum has been based on the assumption that incoming signals may be represented by sums of plane waves. Since plane waves are represented by vector wavenumbers the frequencywavenumber spectrum determines wavefront direction and, for a given temporal frequency, the velocity of the wavefront. In many applications such as measurement of seismic events from great distances through the earth's crust, the assump tion of plane wavefronts is a good one indeed. For electrophysiological data, where smaller dimensions are found, the assumption of plane wave fronts is probably less valid; however, such an assumption may not invalidate the usefulness of the computed frequencywavenumber spectrum. An analysis of the frequencywavenumber method applied to spher ical wavefronts indicated that velocity and wavenumber magnitude are determined but wavefront direction has little meaning. This comes from the fact that the wavenumber for a spherically symmetric wavefront is a scalar quantity. Thus if the spherical wavefront is given by 2ri(k rf t) S(r) = e r (49) where: k = scalar wavenumber (cycles/cm) f = temporal frequency (Hz) o = 27rf 0 0 then it is shown in Appendix C that the frequencywavenumber technique gives a delta function about k = k and X = . P(X,k) = 6(AXo) 6(kko) (50) where: P(A,k) = frequencywavenumber spectrum From the results given, it is clear that the frequencywavenumber technique does not determine wavefront direction for the spherical wave front. For the case of spherical wavefronts, the objective should be to determine coordinates for wave sources relative to the array of sensors. Accomplishment of this objective requires the use of methods other than the frequencywavenumber method. To apply the frequencywavenumber method to the analysis of electrophysiological data it is assumed that the signals may be represented by approximate plane wavefronts resulting from sources at distances large compared to the array size or by a summation of plane wavefronts. G. Spatial Prefiltering In order to reduce aliasing effects it is important to know the signal bandwidths (spatial and temporal) to establish maximum sample intervals. It has been established that most EEG temporal frequencies are below 50 Hz with most of the energy at 2 Hz or below for averaged visual evoked responses [18]. No limits have been established for spatial frequencies; therefore, some exploratory flexibility is maintained for initial experiments with spatial signals. Spatial prefiltering to limit wavenumber content and reduce errors due to aliasing in the spatial frequency domain was considered for applica tion to the experimental data. The maximum wavenumber content allowable without spectral overlap was determined by electrode spacing for the various experimental data runs. Two types of experimental data were used: 1. human visualevoked response data recorded previously by the University of Florida Visual Sciences Laboratory, and 2. data recorded from the neocortex of rats. Minimum electrode spacing for the visualevoked response (VER) data was two centimeters and had been determined previously by crosstalk considerations for scalp measurements [19]. Minimum spacing for the rat data was determined by available precision for array construction to be one millimeter. Digital spatial filters for prefiltering the experimental data were not used for two reasons: 1. restrictions on array geometry (uniform sampling rates) imposed by available multi dimensional filtering techniques would make application of the filter to the VER data impossible and 2. the current stateoftheart for implementation of digital spatial prefilters is such that their characteristics are not well understood. Rather than introduce ambiguity and possible distortion into the experimental data, an alternate method for checking wavenumber content was used. The method involves a comparison of the wavenumber spectra for a given data set which has been computed for two different spatial sampling rates. Such a comparison was made for data from a rat in which three by three matrices of two and one millimeter grids respectively were used for recording the data. This method allows detection of aliasing errors which would indicate that a higher spatial sampling rate was required. It is expected that incorporation of spatial prefiltering into the analysis technique would eventually be desirable. For further information on digital filtering and spatial filtering the reader is referred to [20] and [21]. H. Cost and Computation Requirements Many factors were involved in the determination of computation time and cost for individual computer runs; however, some factors appeared to dominate over others. Computation time and therefore cost was proportional to the number of electrodes, the total number of wavenumber points in the output spectrum, and the total number of sample points of the spatiotemporal field (K x M x N). For the band pass computation method, increases in N should have less effect because of the use of the FFT algorithm in computing power spectral estimates of the N temporal data samples. This was verified by experimental results. The savings in computation time by use of the bandpass method was low but is expected to increase as the number of time samples is increased. Shown in Table 2 are some of the time and cost figures for various experimental runs. These values could be reduced by the use of a binary program deck which would eliminate compilation time. This was not done for the experimental investigation because of frequent program changes. Also, the requirement that several arrays in the pro gram must be dimensioned exactly make the use of an object deck impractical 00 0O0 po H U U2 g o. O E 0 r< El H 0 *0 COl z 0 Na 5ii H r4 i cn c1 H H r H " 00 *3 0o (N 3 000I * * H H rH  d* .~ c. r %0 n r1* a \DC' \O 47 unless a fixed electrode geometry were to be used in analyzing numerous experimental data trials. The first two lines in Table 2 are representative of the two dimensional simulation experiments. The VER data program included an interpolation routine which added to the required computation time. CHAPTER V FREQUENCYWAVENUMBER SPECTRUMSIMULATION The validity of the concepts of the frequencywavenumber analysis technique was tested via a series of simulation experiments. Characteristics of the frequencywavenumber spectral estimate were observed for known signal and noise conditions. The effects of certain parametric variations were demonstrated by these experiments and theoretical concepts were verified by the results. Because of its relative simplicity and lower cost to implement, a twodimensional analysis was used for those experiments in which the principle to be demonstrated did not depend on the use of three dimensional computations. Relative cost versus performance was dis cussed in Chapter IV. The simulation experiments consisted of software implementation of the frequencywavenumber analysis technique. The computer language used was FORTRAN IV and runs were made on an IBM 360/65 computer. Listings of the various programs used are given in Appendix E. Storage requirements for the program were large (120,000 locations for a twodimensional analysis of a 16electrode array and 256,000 locations for a threedimensional analysis of a 16electrode array); therefore, implementation on a small computer would not be feasible without suitable program segmentation. For fewer electrodes this might be possible and could provide useful results. However, as indicated by the results of simulation experiments, the technique worked best for large numbers of electrodes. A. Array and Signal Example Parameters The array used in the twodimensional simulation experiments is shown in Figure 9. It is a twodimensional rectangular array of sensors with equal spacing in both dimensions. The 4 x 4 (16electrode) array was used for most of the experiments; however, the 2 x 2 (circled electrodes) array was used for some of the experiments to test the resolving power of the array for limited spatial information. Shown also in Figure 9 are the signal parameters used in the twodimensional simulation experiments. Note that the spatial and temporal sampling rates were conservatively high at eight times the signal frequency and wavenumber instead of twice the maximum frequency as required by the sampling theorem. For analysis of practical and nonartificial data, sampling rates slightly higher than twice the Nyquist rate were used to reduce aliasing effects caused by nonideal bandlimited and wavenumber limited signals. Application of the frequencywavenumber technique to electrophysiological data is dis cussed in Chapter VI. It is emphasized that the achievable resolution for a given number of sensors and a given number of data points is improved for sample rates nearer to twice the Nyquist rate than that achieved with a sample rate of eight times the maximum signal frequencies. Computer results consisting of twodimensional plots of the frequencywavenumber spectral estimates for the various simulated condi tions are included in this chapter along with perspective plots and crosssectional views. The directions shown in Figure 9 for the vector wavenumbers were selected so that a vector drawn from the peak value in the frequencywavenumber spectral plot to the wavenumber origin represents y X_ X X X I I I I I I I I I [ye = 2 cm. IX X X  4  .4  xS = 2 cm.  r a) ARRAY GEOMETRY S1 = cos [2lr(flt  S2 = cos [2r(f2t  53 = cos [27T(f3t  Bix yly)] 2x Y2y)] 03x y3y)] COMPOSITE SIGNALS S = Sl + S2 S = S1 + S3 TS = 0.01 sec = sample time fl = 1/8 Ts 12.5 Hz f2,f3 = 12.5 40 Hz variable 81 = .0625 cycle/cm Y, = +.0625 cycle/cm Sl 02 = +.0625 cycle/cm 83 = +.0625 cycle/cm Y2 = +.0625 cycle/cm y3 = .0625 cycle/cm S2 S3 b) SIGNAL PARAMETERS Figure 9. TWODIMENSIONAL ARRAY GEOMETRY AND PARAMETERS the wavefront direction. The particular convention selected is not important except that interpretation of results must be consistent. Array geometry and signal parameters used in the threedimensional simulation experiment are discussed in Section G. Results of the three dimensional simulation are given in the form of a series of twodimensional wavenumber plots for x and y wavenumber components with the z wavenumber component being constant for a given plot. B. Comparison of High Resolution and Conventional Estimates For the first experiment a comparison was made of the obtainable resolution by the high resolution estimation procedure relative to that obtainable with the conventional estimation procedure. A perspective plot on a twodimensional wavenumber grid is shown in Figure 10 for results obtained by the two methods. The peak for both methods occurs at kx = .0625 and k = +.0625 as it should. It is clear x y that the high resolution method provides a sharper peak and thus improves the resolution over the conventional method. One important point to be emphasized is the role of noise in the high resolution estimate. As was discussed in Chapter IV, a small amount of incoherent noise is added to the cross power spectral density matrix to insure that its inverse exists. The level of the injected noise has an effect on the resolution obtained in the high resolution estimate. Thus the results of Figure 10 are noisefree for the conventional estimate and noisy (10% in this case) for the high resolution estimate. It is shown later in the experimental results that the high resolution method can be improved even more when the noise level is reduced. k 0625 x y / I / CONVENTIONAL FREQUENCYWAVENUMBER SPECTRAL ESTIMATE k .0625 Sx .0625 HIGH RESOLUTION FREQUENCYWAVENUMBER SPECTRAL ESTIMATE Figure 10. CONVENTIONAL AND HIGH RESOLUTION PERSPECTIVE PLOTS Shown in Figure 11 is a twodimensional plot of values for the conventional frequencywavenumber spectral estimate. The values shown are in db attenuation so that the peak occurs at a value of zero. Only integer values are printed and the actual spectral values are rounded to their lowest integer values. Figure 12 shows a plot of the high resolution estimate under identical signal conditions. For clarity values of attenuation larger than 20 db are not printed. Figure 13 shows a section view of the results obtained by the two methods. This is a plot of spectral values along the k = .0625 line. C. Array Geometry From the sampling theorem it can be shown that resolution is inversely proportional to the data length. For temporal signals the data length is taken to be the record length in seconds and the resultant effect is on frequency resolution. For spatial transforms the data length is equivalent to array aperture in meters and the resultant effect is on wavenumber resolution. An experimental run was made to test the effect of array aperture on wavenumber resolution. In one case a 16electrode array was used with an aperture of 3xs by 3ys (see Figure 9). In the other case a 4electrode array with an aperture of xs by ys, that is, one third as large in x and y directions as the 16electrode array was used. A section view of the results obtained from the high resolution method for these two arrays is shown in Figure 14. The threedb band widths are indicated and it is observed that for equal noise levels an aperture which is three times larger yields a resolution bandwidth 4 .4 f t 4 4 r4 f .W .4 4 4 M M m af 4 4 4 4 4 r in a .D dD "4 > F. I n a 4 4 4 4 * 4 4 6 * 4 4 6 * 6 * a 4a P. P. 01 4 4 4 4t .4 N N I05 4 4 4 4 40 (a .C 0 in In I IA P. m  10 cP 4 0 a. Ol 4  cl on 0 0 N N * . * 0 * a * a I9 0 Ia r(  N .0 P. 9 6 4 44 00 ar . 0 0 0> @. 0 4 S o * a a * a * 6 * 5 * . * S *~ S *~ a * S * a, * aI  . ~ p C ~ ~C ~ * N 4 o * Nm o 4 N *.4 , < f M 4> M m n N N all 4 4 4 alt 0 0.4 N Itt O (Bo M> O o o r f a o 0 4 N N N alt r4 P O 4 oa Ut U N T N 40 0111 49 C 4", I 0 4 4 m a. Irj > N P 4 4  4 4 4 4 4 4 4 44 4 4 N P. 4 4 ma m t m N 4 4 N r 4 >4 N .V 4 4 N rj r 4 4 N .4 N 4 4 '4 m m 4 4 rI be SA a 00aa a 0 a a 0 l C S S I u I" S w I S I I I I I I I I g l In a 00 0 a 0 a0 o a a 0 a 0 a 0, n an in in a a a m a in cA w '0 a w a 4n 44 4 Im Io em in V o C4 w MI 0 P4 k" P. 0 4 4 4 4 4 4 4* 4 4 4 4 4 4 4 0 a0 o J o00 0 0 0 a I I I N  N .* 4 4 4 4 *' 4' * S S * a * * S N N 4 4 4  a. a a S0 wt in I i *4 1 w4 4 04 4M 'I N> 4 0 I a n *o w e W 0' p. III 4 < i F F N MI F 4 .4 .4 4 a. af m '0 P 4 M < *J l N i m .o 6 4 a me a o .4 .4 .4 .4 .4 I 4L 1 4 0 6 a i .4 .4 .4 .4 .4 . p.. F P. 6 . 9 I zr I H J aa * a a a a a a a a  * 0 a a a a * * a a a a a o a a a a a n U U 6 6 S *W a a * a a a a a a * 0 a a S a a a a c o a a a a * a a a a a a o a a a a a a n N S a a a a S a a a a a a a a a 0 0 a a a a a a N N a a a a a a a a a a 4 a a a a a o a a a a a a a a N a a a a a a a a a a a a a a a a a a 0 a S a a a a a N a a a a a a a 0 a a a a S a 4 a a 44 a 0 a a a a a a C4,' 60 S007 ? ? ? ? o o 00 o a mu o u n o mu 0 mu o urn o urn u mu O w ri a a a a a a a a a a a a a a a m a SN a 0 m Im am o a @ B n * a a J o a a I I a I I a a 1 a. 0 0 0 00 0 C 0 I I I I I t    1.0 k = +.0625 cycles/cm Y 16 electrodes S Conventional .0625 0 +.0625 Figure 13. CROSSSECTION OF CONVENTIONAL AND HIGH RESOLUTION SPECTRAL ESTIMATES 1.0 k = +.0625 cycles/cm 0.707 high resolution only 0.7074 electrodes ^4 electrodes .0625 0 +.0625 RELATIVE EFFECT OF ARRAY APERTURE ON RESOLUTION Figure 14. which is onethird as large. Thus experimental results agree with expected results as implied by the sampling theorem. Additional indication of the resolution dependence on aperture is demonstrated by comparing Figure 12 with Figure 15 which gives the high resolution attenuation plot for the 4electrode example. D. Temporal Window Functions Two temporal window functions were tested to determine their relative effects on the frequencywavenumber spectral estimate. The rectangular and cosine tapered window functions used in the simulation runs are shown in Figure 16. The temporal data segments were multiplied by one of these two window functions before Fourier transformation both for the discrete frequency computations and fast Fourier transform (FFT) computations. For the discrete frequency computation method and a relatively high incoherent noise ratio the resolution of the frequencywavenumber spectrum was the same for both rectangular and cosine tapered temporal window functions. For a ratio of incoherent noise power to total signal plus noise power of R = 0.1 and discrete frequency computation, the frequencywavenumber spectrum of two plane wavefronts was identical for both the rectangular and cosine tapered temporal windows. The spectrum for the rectangular window is shown in Figure 17 and the spectrum for the cosine tapered window is shown in Figure 18. As the incoherent noise ratio is reduced, some smearing is introduced into the twowavefront spectrum for both rectangular and cosine tapered temporal windows. Figure 19 shows the spectrum for a rectangular window with an incoherent noise ratio of R = 104 a rectangular window with an incoherent noise ratio of R = 10 i n in ec n iln i ne ni .4 4 .4 .4 4 4 . * 4. .4 .4 f* 4 .4 .4 n m m n .4 .4 N N fr 0 .4 N N a> a *I 01 * .4.  <* 4. 4. .4 .4" 4 .,4 .4 .4 .4 .4 .4 O> 0> 0 * *4 .4 l P 1 4 4 4 N N * a 4 a i o4 * N U 4 i n N * N 4 tn it . , U rN N, 4 4 .4 n i 4 0 f " in 4. I o 0. U U N U U P @ .4 .4 .4 # 9% L .4 P4 .4 * 4** .4 .4 .4 ' <* .1' .4 .4 ^4 * .4P .4 .4 .4 Al on M Nx iN l  .4 .4 Nu H i .4 4 . 0 .4 N v> 0 .d .4 .4 0. 0 .d N 0. 0 .4 N* .4 .4 .4 0. 0 .4 0 0 N 0 04 .f N .4 .4 * 0 N .4 .4 . O .4N .4 4 . .4 .4N N .4 4 . S0 0 .4 .4 * .4 .4 0 0 000 000 00 000 n 0 inno am aonno a a n in N N n N 0 N U NU R C N nCO N .  .. 4 n i N N n 4 i . * 4 4 4 4 4 4 4 4 4 . 4 N 4 a U. u N 0 L O *' 2 a 0 .4 N 4 > I Ma M 4n~ in m in 4 "4 .4 m in in .4 .4 .4 mm en en mnm .4 .4 .4t * en * .0 4A ol W V .4 .4 .4 * .4. c4 4 . in..~ .4 .4 . .4 IC M 4 0 4 4 M in .. m ^ ** m * w4 ^ ^ ~ ~ ~ ~ ~ II ~ ~ Y u u ~p g~CLY ~ i   I I I I I I I I I 1 4 t RECTANGULAR TEMPORAL WINDOW 0 0.1T 0.9T T t COSINE TAPERED TEMPORAL WINDOW Figure 16. EXAMPLE TEMPORAL WINDOW FUNCTIONS 1.0 1.0 r r r in r( 44 m n .4 .4 M 0 i n .4 .0 *j m N1 45 .4 . In en in 0 *0 a . .4 .4 .W .4 4 e o N . 00 0.4 en .4 .4 .4 .4 .4 .4 .4 .4 .  a a tAc .4 .4 .4 .4 .4 .4 * N .400 'd N . en .4 o 0 0 .4 r .0 .4 .4 e.4 cn O 4 c Mr o N 0 .4 .. 14.40 .4 0  en e n .0 mA .a 9 .4 4 .4 .4 4 . * r enWr 4 4 .4~ .4 .4 .4 .4 .4 .4 p., f f P. I. rl rl C 4 .4 .4 fI rl .4 .4 p. p. U '~ .4  9. p .4 .4 .4 .4 f 9 r. p. i r) .4.4 r  .4 .4 .4 .4 r. p. .4 .4 p. .0 .4 .4  .. s .. a a a v ".44 .4 .4 .4 .4 .4 .4 .4 .4 .4 0 a "s 3S g o o *guul gO IU U il 1 9U IU UI I U UI 11 1U UJ U LU Ui 491 49 UI 19 UJ I Sa a a a I a a a S S I S s I I S w S u 2 ~ ~ f a a 2 2 2 7 7 a? a a A 0 0 m fa W, o :m0N6 .0 in in A F. an o a 0 0 a 0 a 0 a a o 0N p. p. .4 d4 p 4 .4  a0 m .4 .4 .0 .4 d .4 .4 4 M .4 . 6. 4 l .0 re .4 . 4 m ) *4) t  *o m 1 edae ,re .4 .4 4.0 .4 .4 cu H E4 0 Ori z 0 H PI a, L  61 A p 0 0 .0 ft .0 'a 4 P. F P. .0 0 a0 a0 0 0 . .4 .4 t .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 N. .0 .a N I' A 10 a0 oP a a .a a 0 .4 4 4 .4 .4 .4 .4 .4 .4 ps0 d0 .0 P .0 P4 sC_( 9H ~ ~ ~ nP A COO 10 4 m ft a U. N r; P. a a CD do 0 0 .a NO N m P. r 40 40 3 4D 0 cc a SMIf SQO "M 4a4 3 2 2 a 2a2 .4 .4 .4 *J .4 .4 4 .4 .4 .4 .4 .4 .4 .4 .4 .4" ~Z~ fN.. 0 .0 U .0 .0 to o 0 .4 4 .4 ~ ~ ~.4 d .4 0 P4 'd ~ 4 .0 O N X 8 N .0 F .0 .0 0 2 M 4 4 .4 .4 .4 4t .4 .< .4 4< .4 .4 .4 .4 . 4 4 M 0. 4 ON .0 F P. .0 .0 00 .0 .00 S S a > m .4 O .4 .4 .4 .4 B. B a 2 24 2 S 2 2 aF .4 4 .4 .4 .4 .4 .4 .4 1 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 r 4.0 .0.40 0"t o .4 .0 p. p. .0 .0 .0. .0 20 2 S 2 0 ~E'E .4 4 "4 .4 .I .4 .4 .4 .4 .X N ft w r. a 40 u UC L)paI E M n O 4 .0m O p pa .0o w. .2 22 >tf >J >4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4 .4M * 4 6 0 0 M ft 4 in 9 a 4 N go 410 a2 .4 .4 .4 .0 .4 6 4 a / s 4 0 0 .? T 4 .. p 0 0 0 . 4 a 40 N a4 N F 0 .0 0 0 a0 m F 0 a 4 0 to In .0 .A.0 0 a 0 p. 4.m f in 0 ca n m4 in p m 2 0 2 22 8 an. ft inoP~r H a o o o a o 0 0^ 04 .4 .4 .4 .4 .4. .4 .4 .4 .4 .4 ^^ 0 I I I S 8 6 3 0 I 'i r( <( fri l r) rf8E ~ " .0 .0. .0 UA .0 Mi .0 i 0 .0 Mi MI .0. Mi M< i E .000080008 0C090Q 80 0.0 N .4 Op IA N Op IA N N p. ON 20 p. 2 24 N .4 .4 .0 p. .4 .0 4 1 0 < .4 N < .0 0 4 p. 0 < 6 <  5E 3 00800808080 ~~ 0 9 nO S 2 2 S 22 S2 :S 2 2 22 22 o S9 n n~4 2 SSS S22 lS 222 222 3~ 9 Ere 8 g ^? ^S gg oS g~~'~ a s sss gss sss g j^ Sl si s sso ssg ,,oq Y 3 : ". '' 0 ' a ~ ~ uur lu o oooo I*t 62 Is W Pa 0 , a a a a a a a a : : : : : : a a a a a aa a i af a a 0 0 a ta a tA 0 a In a N (A : a : : : a a: a a : a a a a a : : : : a a a a a aa Sa a a a a a a a 4 aa a a a .. a aaa a aa a a a a a a a a a a a a a a a a a 0 ~ a a a a a a a a a a a a : a I a a a a a a a a a a a a a a a a a a a a a a a C a a aa a a aa a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a 2 ~ a a a a 14 a a a a a a a a o ~ a ~ a I ~rC a (~ a B CI a a a a a a a 0 a a 4 ~ a~ a aaaA a a a a a a a a : : 2 a a aa a a aa a a 0 w a a a 0 a a a a p. a a a a l a a a a a a p .. a r ~ a  4 ~ r a a a aaa aa a a a a a a a a a a a a a a a a a a a a a a a a a : a a a~~ a a a a a a a a a a a * a a a ~ a a a a a a a jO a ~ a a a a a a a a a a a ~~ a~ a aa aaa a a ~ a I a aa aaa a a I a a a a a a a a a a a a a a a a a a a a a a a a~ a a a a a a a a a a a a a a ~ a a a a a a a a a a a a a a a~ ~ a aa aaa a a a 1 a 0 ~ 0 00 0 00 00 00 00 I P II SI I II IS lAB LU W ULUW iU w LU4 4 W BA ULUL L U U o~ u 0 0 00 0 00 U ~f Al 0 WI0 A 0 L 0 U W M 0 II I l 4 Pd .4 0 I N l Al N WI7 0 AlL t . 0 P. 1 1 1 Al 7 ~II .0 P 0 a a B a a * 00 00 0 00 II SI I I Figure 20 shows the spectrum under identical conditions except that a cosine tapered window function was used. The figures show that smearing is only slightly higher for the cosine tapered window. The above results apply only to the simulation runs where discrete frequency computations were used to compute the cross power spectral density. For the bandpass method, the results are different and less predictable. It is shown in Section F that for the bandpass spectral computation, the cosine tapered window function causes less error due to sidelobe interaction and produces a frequencywavenumber spectrum with more clearly distinguishable wavefront directions. Further quantitative effects of incoherent noise ratio on resolution are discussed in Section E. E. Additive Incoherent Noise The relative effect of additive incoherent noise on the resolution obtained with the high resolution frequencywavenumber technique was tested by adding known amounts of noise to the cross power spectral density matrix before inversion. For the example signals used, this proved to be a necessary step to insure non singularity of the matrix. Values of the ratio of incoherent noise power to total signal plus noise power of R = 0.0001, 0.001, 0.01 and 0.1 were used to determine the effect of incoherent noise power ratio on resolution. A crosssectional view of spectra obtained for a single plane wavefront and for various values of R is shown in Figure 21. Shon in Figure 22 are the spectral profiles for a signal composed of two plane wavefronts. 64 (4 a P0 0 0 0: E0 W 5 r4 ** a a* a a a a a a 06 0 a a a a a \ i : U0. I us *aa a U a a 0 a a a a I 0 S % 0 Mj 0 t a a a m P. a m P a on P4 : : : o e2 a 0 0 0 { a a a a a a o : : a a : a : a ..a : a :" a: a a a . a a 0 : a : a * a a a.. a a a a a a a a a a a a a a a a a a a. .a. a:a a a a. a a : a a a a a a a a a a a a a* a a a a a 0 a aa M a.a a a a a a a a a a Sa a a a a a a a a a a ... a:aa a a a 0 aa a 0 a 0 a 0 .0 a a a a a a a a a a a a a z a a a a N .a : : : a a a a o a a a a a a a a a a a a a a a a a 0 a a a a a ~. a a a a 1.0 Figure 21. EFFECT OF INCOHERENT NOISE ON SINGLE WAVEFRONT 1.0 .0625 +.0625 Figure 22. EFFECT OF INCOHERENT NOISE ON MULTIPLE WAVEFRONTS Ideally it is desirable to have zero noise level; however, the matrix inversion requires in many cases a finite noise level to insure nonsingularity ( See Appendix F ). This was true of the simple plane wavefront example signals used in the simulation. For the example signals and the matrix inversion subroutine used, a minimum value of R = 6.0 x 105 was found to be necessary for the matrix to be non singular. However, because of the spectral smearing at low noise levels a value of R = 0.001 to 0.01 was found to be best. For data recorded from practical sensor arrays it is expected that any superficial noise addition would not be necessary because of the inherent noise in the data and its general complexity. This was verified by results obtained from application of the frequencywavenumber technique to the analysis of electrophysiological data. As shown in Figures 21 and 22 resolution is improved for both single and multiple component signals as the noise is reduced. However, there is a nonlinear interaction present which causes the resolution to improve less if multiple wavefronts compose the signal. Figure 23 shows the high resolution db attenuation plot for one wavefront incident on the 16electrode array. The incoherent noise ratio is R = 0.001. The improvement in resolution for reduced noise level can be seen by comparing this figure with Figure 12 where the noise ratio was R = 0.1. Figure 23 shows that even the closest wavenumber points to the actual signal values had spectral magnitudes that were greater than 20 db down from the peak value. Further reductions in noise are possible; however, noise ratios less than R = 0.001 were found to introduce smearing into the spectral estimate in the form of sidelobe enhancement when more than one wavefront comprised the signal. 67 a a a a a a a a v a a a 00 1 4i 0 a a 0 a 0 E4a a ot 04 P4 P4~~~ .4 4  a a a a a a a a a a a a 0 a a vi us usw w in in a a a in o in n in 0 w a m a w m m w M ". 0 m i P 0 a. a M f a a : a a .4 : ~~ j j : : r ; : Ir 1 s a a." a a a. a a a a a a a a a s a a" a a a a a a a a a a " ~ a a a a a a a a a a a a a a a a a a a a a a a a a a a a a : : : ::I I : a a a a a a a aT a a a a a a in am an a ai S3 aS 3 S> a S a a a a> a a ar a a a 3t a 3 3 3 a a a a a a a a a a 3 a a a a a a a a a a a a a a 1 a a a a : a: 'a a. a a. a a * a a a 3 a a a a : a a : a a a a : a a : : : :" .. : : a : a a U U a a a a a 3 3 U d U1 a a Spectral plots for twowavefront signals are shown in Figures 24 and 25 for incoherent noise ratios of R = 0.01 and R = 0.001 respectively. Cosine tapered window functions were applied to the temporal data before transformation. As was mentioned in Section D the spectral smearing seen in Figures 24 and 25 was only slightly greater (within 1 db) than that observed with a rectangular temporal window. F. Bandpass Signal Spectral Computation The bandpass spectral computation method described in Chapter IV was simulated by applying the method to signals representing single plane wavefronts. For a signal represented by S1 in Figure 9 the frequencywavenumber spectrum was computed using the bandpass method for a passband from six to twenty Hz. The spectrum was computed for both a rectangular temporal window and a cosine tapered temporal window with respective spectral plots as shown in Figures 26 and 27. Spectral sharpness is best for the rectangular temporal window as it was for the method using discrete computation of the cross power spectral density. However, for signals containing multiple wavefronts and for low incoherent noise levels the bandpass spectrum using a cosine tapered window is more reliable than the bandpass spectrum using a rectangular temporal window. This conclusion is based on simulation experiments for the two window func tions in which the two wavefronts were more distinguishable for the cosine tapered temporal window. This might be expected since cosine tapering is a form of smoothing and should produce a spectrum with less error due to sidelobe effects. An unexplained dependence on specific wavenumber content of the computed spectrum of multiplewavefront 69 a. a a a a a. .a a ; a a a a a a a a a w w w w Ii 1 1.: : *: : : : :, : : : : :a a: : : : t a n a0 n in 0 M a a a a 0 ,. ma M m a a a * a a a a a a a a a S a a a : a a a a u a a a a: o a :0a a q a a a aa a^ a a a a a a o~~~~ a o a a a a a a a a a a a a a a 70 * a a a a a aa0 a a In tti 1i" l a aa a a a a a g a a a S aS a S a a a a a a a a a a a a a a a a a a a a a a a a o a a a a a a a a a a a a a oa a a aao . a a .a. a a a a a a a : a : a .a a a a. a a a a * a : m a a a a a a a a a a : ..aa aa a a a a * a : a a a aa. a a aa * a a a a a a a a a a * a a a a a a a a a a : S a ! : ! S: I ,o o : 2 : : o . a: : : a a a a a : : : : : a a a  a a a a a : a a a o a a a a a a a a a a . a a a a a a a a *a a. a a a a a. . a a a a a a a a a a a a a a a a a a a a a a 0 0 0 .4 .+ a 0 0 o 8* o ? ' o 'o co a .4 .4 P int ; a 4 4; a *M a* a* a a a a a a a a a a a a' * a a a a a a a a : : a a a a a a a * a a a * a a a a a a S a a .4 a a a . a a a a a a a a a a a *~ a a a a a a a a a a aa *~ a aa a a a a a a a a a a * a aa a aa a a a a a a 11 a * o n 0 * * * 0 I a a a a a . . .a a . a .a .a a a . a a . a .4 4 .4 . a a a a a 000 00 0 in o n us II o iM Wn a a n a a * f a U I I I I a a a * * * * * a a a a a * a * a a a a a a a a * a * a a a a a a a a a * a * a * a1 a a a a a a a a a a a a a * al 0 0 a oo a 0i 0 "1 n .10 IAn F o ** N * * c a I S .I I  An An n An An An An A m m m : .4 .4 .4.4 4 d n An An An .4 .4 .4 .4 An W% In In .4 .4 .4 4 in in n A .4 .4 .4 .4 An In AA In .48 .4 .4 d4 .4 .4 d An An 4 .4 4 .4 .4 .4 .4 .4 .4 m .4 O . .4 .4 .4 .4 01 N .,4  ,. 4 A l 0 0 .4 4 . .,4 .4 .4 w 4 *4 N 4 @ 4 9 i01 N An 4 401" An An An An .4.4 .4 .4 Mn M m in m M< m .4 .4 4 .4 in An m in .q .4 .4 .4 An An An An An Ai An An .4 .4 .4 .4f An An An An .4 .4 .4 .4 An An An An in An An in !1 4 W. In 4 $A .4 a a a 4 in In & in in i in m N S m in An ,4 .4 n In & ,4 4 .4 .4 * in __~.4 .4 .4_ __ N nj N N 4 4 4 .24  f. a a* cmx .B P. 0l .4 N < C4 N.N A < 0 0 N .4)~ .4 An m An 0o 0 44 P .4 .4Y .4 .4  0 0 0 .4 N* 01 M 41 41 N 01 4 4 4 4 .4 .4 .4 . 01 01 01 01 .4 L .4 .4 .4 .4 An m An A .4 .4 .4 4 w m in m .4 .4 . n An An in .4 .4 .4 .4 In in in An An An An In .4 .4 .4 .4 In inm in An in An In An An An .4 .4 .4 .4 %n In in An An An An An An An An An .4 .4 .4 .4 40 in %n In 4 .4 .4 .4 in An in An .4 .4 .4 .4 An An An A I4 In .4 .4 n An An An An An ml An An An An An .4 4 .4 .4 An A An An An n M An 1 .4 .4 .4 in An m L .4 .4 .4 .4 an An An A .4 .4 .4 .4 m 9% n in in m A A 14 .4 .4 .4 m m An A .4 .4 "4 .4 iA m A A An An An n in An In m .4 An .4 .4 An An An An i n in In An .4 4 .4 > A n Al jt 4n in An 9% oA An %n o on An in In An n An An 4 .4 I4 .4 An An An 41 An An An An .4 .4 .4 . 9 0 0 a 0 4 4 4 4 4 4 0 0 0C 0 000 0 00 a Q 0 0 0 0 0 000 1 tawI iiiI I 1 uS 6 6 S S S S S S I I S lU LU lU LJ LU lU UJ UJ LU J LU LU LU LU W LU LU LU lU LU o An 0 o o o a a o 0 a a a a 0 0 Ao M N o A 0 on an oin O in O A 0 N A n N ^ 0 41 An N o 41 An NM N' An S 0 N An P. 0 .4 N! 4 .4 .4 4 P. W6 I C a I 0 z U& 0 P. 4 s _ _ signals was also observed. For example a signal composed of wavefronts S1 and S3 (see Figure 9) was successfully resolved by the frequency wavenumber method; whereas, a signal composed of wavefronts Sl and S2 produced an ambiguous spectrum for the same frequency passband. It was found, however, that judicious selection of the passband effectively reduced ambiguity in the frequencywavenunber spectrum. For example, the composite signal of S1 (12.5 Hz) and S2 (15 Hz) was poorly resolved by the frequencywavenumber spectrum for a passband of six to twenty Hz regardless of the value of incoherent noise injected into the signal or the choice of temporal window. The same example was run with a different passband (11.5 to 16 Hz) which produced an improved estimate of the frequencywavenumber spectrum over that obtained with the larger pass band. Figure 28 shows the narrow bandpass spectrum for a rectangular window and Figure 29 shows the narrow bandpass spectrum for a cosine tapered window. A sharper spectrum is achieved with the rectangular window; however, a more accurate spectrum is obtained by using the cosine tapered window. The results shown in Figures 28 and 29 are consistent with the expected result of cosine tapering, which is a form of smoothing. In general, smoothing produces a broader spectral estimate with higher statistical reliability and less error due to sidelobes. Thus the results of experiments using the bandpass computation method indicate that the selected pEssband should be narrow and carefully selected to improve the frequencywavenumber spectral estimate; this may call for spectrum analysis prior to data processing. Also for improved reliability of the spectral estimate, a cosine tapered window was better than a rectangular window when using the bandpass method. 4 a IL I o a a a a a a a I ~ ~ I * a a a a a a a a a a a a a a * 8 a ai a * a a i a * a a a a a a a a a a a a a a S* a a a * a a a a BC i B e e e A * ao a * 04 Na a a e a a a a a ai a a a a *a a a a aO o a a a N a a a a a a a a I a a a a a a * a a a* a l a a a a * a1 a aE. * i a a t * a* a a4 a a .t 0*, o .4 a ~a a a a a a a a 0 0 0 0 0 a 0 m % a w 0 UJ LU UJ LU LU Ii SM 0 0 0 0 0 0 0 a I I 0 IA 0 In 0 0 f IA N 0 F LA 0 0 0 0 0 0 N a a 0 a a u4 a 0 in N 0  * a 0 Sa m = m a E  * . * a * a * a * a * a * a * a * a * a * a * * * a a a * a * a * a * a o a N * a * a * * * a * a a a * a a a a 0 * a * a a a * a * a * a * a * a a a a a * a * a . .. .l .4 4 6 0 0 o 0 0 0 0 0 0 0 I I I I UJ J U UJ UJ lU UJ IU 0 0 0 0 0 0 0 " o 0 ;0 i n N LA r 0 N LA P' 0 ' i* I I I *S I* I *m i z A I i il  i * a a a a * * a a a a a : : : : : : : : : : * o a M* N .4 .4 .4 : a a . o  o * o 0 a 4 a 0" i4 0 o a 0 a S a a S. .4 4 = ** a  0 0 in 4 i a t t a C A4 .I n< r 't P aa a a q a a awt a a m a  f .4 d 0 4 1 a a a 0 0 0 0 0 0 a a a a0 a a a a0 0 a 0 a 0* a a a 0 a a 0 0 a an a m 0 * 0n 4 P in a4 *m an a a a % A a a a C 'a In ." .m .4 .4 a a a % a 4 a a a a a a a a a a* 4 0 a a a a a a a a . pa 0 0 rN in 4; O C 0 0 8 a a a a a a a aA Recall from Section D that for discrete computation of the cross power spectrum, both rectangular and cosine tapered windows produced similar spectra,with the cosine window producing slightly (1 db) more smearing for a low incoherent noise ratio. G. ThreeDimensional Simulation Because of the complexity and projected cost of threedimensional frequencywavenumber computations and since most parametric interactions were sufficiently demonstrated by a twodimensional analysis only one simulation experiment was designed to test the threedimensional technique. The signal consisted of a single threedimensional plane wave front as shown in Figure 30. The 3ightelectrode cubic array shown in Figure 30 was used in the analysis. Since the signal was a single plane wavefront and by reducing the incoherent noise ratio to R = 104 for improved resolution, signif icant savings in computation time and cost were possible by selecting an elevenpoint matrix of wavenumber values in three dimensions instead of the twentyonepoint matrices used in the twodimensional analyses. Results of the threedimensional frequencywavenumber simulation experiment are shown in Figure 31. Twodimensional attenuation plots were generated for an x and y wavenumber grid for specific values of the zdirection wavenumber. In Figure 31a is shown schematically how the spectral plots were generated. Figures 31b, c and d give the attenuation plots for zdirection wavenumbers corresponding to the signal wavenumber (center plot) and one wavenumber point each side of the signal wavenumber. Because of the high resolution obtained for the threedimensional simulation all spectral values were greater than 77 Electrode spacing = 2 cm.  / I / I / I ^/I I I / / I // / I / / / 7 a) ARRAY GEOMETRY S = cos [2w(ft k x ky k z)] f = 12.5 Hz k = .075 cycle/cm k = +.075 cycle/cm k = .075 cycle/cm b) SIGNAL PARAMETERS Figure 30. THREEDIMENSIONAL ARRAY GEOMETRY AND SIGNAL PARAMETERS p, in I M rz4 P4 * a a a a * a a a a a a a a a a a a a a a a a a a . a a a a a a a a a a a a a a a9 * a a a a a a a a a * , a a a d a a a a a , a a a a a a a a a a a a a a a : = a a a * a * a * a Sal * a * a * a * a * a * a * a 41 a a a * a i a 4 a 2 . a + 0 S uj iu iu S t i S00 0 0 a 0 NI 0000o a a 00 III~ 0 o m m a an a a a a a a o a a o c a o I* s* i Sa a ea a a a a a al a a a a a I tl a a a a a1 a a * a a a1 a a, a a a a a a a a a a a a a a a a a a a a a a al l a a a a a a a S S I 0 0 00 00C 0 ft ajm em m omu 0 o 0 0 0 i 0 0 0 S j j a 0 J i* 9* a a o a ai a a1 a a0 a a a* a a* a W T WZ T T T T T 0 0 00 I nI 8 00 0 0 4 a a 4 *mauit a 2 0 5M S3 2 a *4 08 S. 0, ,0 0 0 sU 0 m :x *  I twenty db down from the peak value as indicated in Figure 31. For this reason, only the three attenuation plots shown are included. Thus the threedimensional frequencywavenumber spectral computa tion is demonstrated. For the application of this technique to multiple wavefront signals and bandpass computation the parametric interactions discussed earlier should have analogous effects. Increases in array aperture, number of data segments, and number of data points which are required for improved resolution and statistical reliability produce a greater increase in computation requirements and cost for the three dimensional analysis method. CHAPTER VI INVESTIGATION OF ELECTROPHYSIOLOGICAL DATA BY FREQUENCYWAVENUMBER TECHNIQUE The frequencywavenumber technique was applied to two types of electrophysiological data to determine if the technique might provide useful information about the data. Human visualevoked response (VER) EEG data and penicillininduced focal epileptic activity in the rat neocortex were analyzed via the frequency wavenumber method to demonstrate the application of the technique to these data. Succint interpretation of the results appears premature at this time; however, the results indicated that certain characteristics of the data might be categorized by application of the frequencywavenumber spectral technique. A. Electrophysiological Signal Characteristics For the data analyzed, a power spectral analysis was first performed to determine the temporal frequency content. Recall from Chapter V that spectral smearing and ambiguity in the wavenumber domain were less for those examples in which the selected temporal frequency band most closely approximated the true frequency content of the signal. The expected wavenumber content of EEG signals has not yet been determined; however, physical constraints of electrode size and array spacing determined the maximum wavenumbers which pro vided nonaliased spectra. Ideally the data should be subjected to a wavenumberlimited spatial prefiltering process prior to application of the frequencywavenumber spectrum. However, a more expedient method for checking the wavenumber content was designed and tested using the rat data. The method consists of computing the frequencywavenumber spectrum for electrode arrays of different spatial sampling rates and then comparing the resultant spectra for differences due to aliasing. The VER data had been previously recorded by the University of Florida Visual Sciences Laboratory with an array geometry designed for determining characteristics other than the expected wavenumber content. For a description of the experimental method for recording the VER and related discussions of its implications see [18], [22] and [23]. Section B gives the frequencywavenumber spectral plot for VER data obtained from a subject with intermittent exotropy (subject DLH). For data recorded from the rat neocortex the spatial sampling rate was restricted by the techniques available for the electrode array manufacture and by the electrode diameter. The stainless steel wire electrodes had a diameter of approximately 0.4 mm. Thus an array of these electrodes in their closest spaced arrangement could efficiently sample signals with wavenumbers no higher than 12.5 cycles/cm. To keep tolerance errors low a threeelectrode array with 2,0 mm spacing was fabricated. Sampling intervals of 1.0 mm were obtained by micro positioning of the array. For the 1.0 mm sampling interval signals with wavenumbers as high as 5.0 cycles/cm could be sampled without aliasing errors. Spectra produced by the 2.0 mm and 1.0 mm sampling grids were compared to determine relative errors due to aliasing. Section C describes the results of the application of the frequency wavenumber technique to the analysis of focal epileptic discharges from rat neocortex. Details of the experimental method employed for collection and analysis of these data are presented in Appendix D. B. Analysis of Human VisualEvoked Response Visual evoked response data recorded from an array of sixteen electrodes were analyzed via the frequencywavenumber method. Individual signals represented averaged EEG responses to visual stimulation of the subject and were recorded from scalp electrodes positioned over the occipital region of the brain. The electrode array consisted of sixteen electrodes located on two concentric circles of radii 2.0 cm and 4.0 cm respectively. Eight electrodes were uniformly spaced on each circle. An electrode at the center of the circles was used as a reference. Electrical response activity was monitored for a period of 500 msec after stimulus presentation and the average of fifty such responses was used as the VER data for a given electrode position. All VER responses for the various electrode positions were referenced to the center electrode prior to analysis using the frequencywavenumber method. Analysis of the VER data indicated a need for higher temporal sampling rates than required by the sampling theorem for the expected maximum temporal frequency content of the signals. This conclusion was based on the requirement that each time record be divided into a number of segments greater than or equal to the number of electrodes (16 for the VER data) and on the observation that short time lapse transient events comprised the visualevoked response. Thus a large number of data points for short observation times required an increased sampling rate. Additional temporal sample points for the existing VER data were obtained by a linear interpolation between the true sample values. An initial analysis of the VER data using the frequency wavenumber method was performed for temporal frequencies of 2.0, 5.0 and 10.0 Hz. The maximum wavenumber included in the analysis was given by the inverse of four times the electrode spacing or 0.125 cycles/cm. The resultant frequencywavenumber spectral estimate for the 500 msec visualevoked response is shown in Figure 32. The fact that all spectral values were at either zero or greater than twenty db attenuation levels is partially explained by subsequent analyses in which the 500 msec response time is divided into two 250 msec time periods. The frequencywavenumber spectrum for the first 250 msec is shown in Figure 33 and the spectrum for the last 250 msec is shown in Figure 34. These computations were for temporal frequencies of 8.0, 9.0, and 10.0 Hz to include observed power spectral peaks near these frequencies. Inclusion of the 2.0 Hz temporal fre quency in the spectral computation for the 500 msec record length probably contributed to the ambiguity observed in that spectrum. Spectra for the subdivided response indicate that spatial signal sources changed with time for the 500 msec observation period so that the frequencywavenumber spectrum for the total observation time exhibited a smearing effect of different transient spatial sources. This observation in conjunction with the required large number of 85 S0 0 0 0 0 0 0 0 0 w 0 co 0 c Sa a S * SQ O O O O O 0 0 0 0 a 0 00 0 0 a 0 0 000 a0 0 0 a 0000 0 0 o a o a0 0 0c 0 0 0 a0 0 0 0 a e o a a o a aoo o o o o oA 0 o ooo o : : 0 0 : :a a 0 a 0 SO 00 0 0 0 0 0 0 0 o S a at a o 0 oa o o S0 0 0 0 00 0 0 0 0 C * : : : : aoo 0 o ao o So o o o o o o o a a o a a * S *3 So o o o a %o o o < o o o o o o 0 o tT 0 O W O O i t 0 N 0 0 0 0  o o o oo 0 0 o o o 0 1 a * o oo a o 0 o o a o oo a o o a o * o o o o o o p o o o o a o c o o a o o o a o o 0000 0 000 0 0 0 00 0 0 0 0 I I 86 0 0 a a 0 0 4 2 0 E4 p UN I' a a a a a P a a a * m c. r, u. 0.. a d a a. 0 a a 0 ~ ~ CD a a a a in 0 a a a a a 0 a 0 a a a 0 a an 0 l 0 U% an a an a% 0 in 5m 6 a in a0 ain m a a a m a a P a P a a P : : a a.a.a.a C.a a a a a a a a a a a a a + a aa a a a aa Aa a a a a a a a a + a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a~ a a a a a a a a a a a a a a ~ a a a 01 a a a a a a N .4 a .~ N N a a* a a a a a a a a a o a> a a a' a> a a a a a a 0 0 F U" E" 0 N I 4 0 @ r* ** a* a* a a >>~I. t4 5A I *. in ~ t r ~ m & ai, a a a. .4. . . ^ a. a a a a a a. a' a ~ o3< n o 0 a af a 0 0 0 o 14 1 0 N I 0 Ea a F a N .4 a .4 .4 N ao a a aa a a .4 o . .4 .4 a a a a a a a e *. 0 a a a a 0 a a a a N N a a a a a ~ a a aa a a a a aQ a a ; a I Z a a a a a a a a a a a a a a a* US a a a a a a a a a a a a a an a a a a a a a a a a a a a a a g a a a a a a a a a a a a a a a a a a a ! a! 11 al ai l a a a a a a a a a a a a a a a a. a. a a a a a a a a a a a a a a a a a a a a a a a 0 0 0 0 0 * a S0 0 0 0 0 a 0 0 o o o a S a o o o o a a o o o a a : a a a a a " a a a a a a a a : a a a a a a a a a a a a a a a a a " a a a a * a a a a * a a o a a a a 4< IL 0 I a J 0. 4 5 4 S O 0 ' .0 .4 .4 a ?a 0 0 0 .4 .4 .4 t .0 bi S oa o a o 000 000 00 '. 4 <    a a S* * o* * 0 * 0 a 0 0 0 ; o0 o a S an a * a U S 0 0 0) 0 0 S000 0 0 0 0 0 0 S o o a a Q a : coon: : 0* 0 ** * 0 0 0 S 0 * S* o* a * 0 a a a a o Ma a vi (.4 0 .4 * . o  .4 .4 0 0 0 a 0 0o i a 0 Ma mu A MA MA Ma ma o a a a a in 0 a oi a oi N vI 0 x mi F a .4 Pd * S pa a a* * r4 S I i i Q   n B I ! data points for input to the frequencywavenumber spectral estimation technique formed the basis for the conclusion that higher temporal sampling rates are required. As indicated by Figure 33 the frequencywavenumber spectrum for the first 250 msec time period following stimulus presentation gave a peak value at k = k = 0. This may be the result of a plane wave x y propagating in the km direction or may be the result of a synchronous active discharge exhibited by neuronal populations under each electrode and triggered by a common source. The latter explanation is considered most likely since similar results were obtained with signals recorded from the rat neocortex. Details of the result of application of the frequencywavenumber spectral estimation technique to analysis of artificially induced focal epilepsy are discussed in Section C. The frequencywavenumber spectrum for the last 250 msec of the VER exhibits bands of uniform maximum and minimum values. Again this may be partially explained by attributing it to a temporal smearing of numerous disparate spatial sources. If the problem is analyzed in reverse by taking the inverse spatial Fourier transform of a flat spectral band, a different possible explanation is obtained. Given a frequencywavenumber spectral plot as shown in Figure 35, the inverse spatial Fourier transform is defined by 00 f(xX) = f P(k,X) e2ITix dk W W (51) = e 2i(kxx+k y dydx oW 2 W o 2 