- Email: [email protected]
Processing noisy line spectrograms as digital pictures
Patrem Nerogaltion
Pergamon Press 1977 . Vol. 9 . pp. 137-146 . Printed in Great Britain
PROCESSING NOISY LINE SPECTROGRAMS AS DIGITAL PICTURES D. G . Niciiot . Department of Defence, Australian Defence Scientific Service, Weapons Research Establishment, G .P.O . Box 2151, Adelaide, S . Australia (Received 22 December 1976 ; in revised form 13 April 1977)
Abstract-Contiguous short-term spectral estimates of non-stationary acoustic signals are often displayed as intensity-modulated frequency-time pictures ('spectrograms') . Some of the problems and advantages in storing and processing certain types of spectrograms as digital pictures are explored . The particular type of spectrogram discussed is from a source producing a time-variant 'line' spectrum consisting of one dominant primary line and several much weaker harmonic or quasi-harmonic lines. This type of spectrogram is commonly produced by rotating machinery . It is supposed that this signal is measured amidst additive broad band noise . The resulting spectrogram consists of a family of related meandering lines, the weaker members of which tend to be lost in the noise . In the proposed picture analysis scheme the raw spectral estimates are stored as a digital spectrogram. This spectrogram is then processed by line enhancement methods to extract the primary line to form a 'family-template' . Using this template the weaker lines can be detected and extracted . By re-assembling the extracted components, a cleaned-up low-entropy frequency-time picture of the signal is produced which shows more structure and which is much more amenable to automatic analysis than the original spectrogram .
Watanabe,IU Meisel1 `o) that high dimensionality of the data is the greatest barrier to successful computer pattern recognition . (e .g.
INTRODUCTION
Overview and summary
The production and analysis of spectrograms (that, is, intensity-modulated frequency-time pictures) is common in many technical studies . Examples of such studies include speech research, vibrational analysis and radio astronomy . The present work proposes a picture-processing technique for the enhancement and data compression of a class of spectrograms containing 'line-structure' as defined in the next section . The type of spectrograms discussed is produced by analysis of acoustic emissions from rotational machinery, and as such is particularly important in the field of automotive research. Here, such spectrograms are used for such purposes as engine diagnosis, tuning and noise suppression . The aim of the picture processing is firstly to detect structure in the spectrogram which is too weak to be visible to the naked eye- This is done by reducing the background noise and enhancing the signal structure . Reducing the noise automatically reduces the entropy of the spectrogram because the noise is `bad information.' (This is possible because we have some knowledge of the signal which is external to the given spectrogram and can use this to distinguish signal from noise .) This partially achieves our second picture processing aim of data reduction . Further data compression can be achieved by removing some of the redundancies in the signal-structure. By such means the noisy high-entropy input spectrogram is transformed into a high-quality low-entropy picture which is much easier to analyse . This data-reducing transformation is particularly important if this analysis is to be done by computer, as it is generally recognised
Line spectrograms
Line spectrograms are produced by sources having 'line spectra', that is, sharp narrow peaks in a onedimensional power-frequency plot such as Fig. I . When plotted as spectrograms, these peaks appear as two-dimensional lines in the frequency-time plane as in Fig. 2. If the signals were stationary then the lines would be straight and parallel to the time axis . Clearly this is not the case in Fig . 2 which are contiguous short term spectra of noise from an idling automobile engine to which band limited white noise (30--800 Hz) has been added . As the engine speed drifts up so do nearly all the lines in Fig. 2. This is because practically all lines in Fig. 2 are multiples of a fundamental frequency corresponding to the cylinder firing rate. We shall use the term 'related lines' to include all lines which are multiples (not necessarily integer multiples) of some fundamentalWe include non-integer harmonics ('quasi-harmonics') because certain components, e .g . noise arising from gearing reduction, can be of this form . In the examples discussed only integer harmonics are present but for generality non-integer lines are still looked for . Now if the background noise is sufficiently high the weaker members of a relatedd family can be lost and, because of the changing frequencies . traditional detection methods (such as time integration) do not retrieve them . Other detection methods such as `order tracking', which is offered as an option on various commercial analysis systems, rely on exact knowledge of the changes in fundamental frequency . 137
1 38
D . G. NtCHOL
LINEAR UNITS
100
200
300
400
500
FREQUENCY (Hz)
Fig . I . Raw spectral estimates from non-stationary source .
This paper describes a scheme for enhancing the weak lines without such knowledge . The only restriction is that one member of the family is reasonably strong compared to the background and to the other members . The other lines can be much weaker and breaks in the dominant line can be tolerated. Implementation A flow chart of the process developed is shown in Fig. 3 . The input signal is sampled by the analogto-digital converters at stage A and the short-term power-spectrum computed by the Fast Fourier Transform (FFT) processor. Succeeding estimates of the spectrum are stored on disc until the spectrogram buffer is filled. Typically the size of this would be
0
100
512 x 300 16-bit words. This buffer is displayed. as in Fig. 2, by a dot matrix plot on a Tektronix 4010 graphics terminal . Only eight grey levels are plotted by the display routine and even this requires sacrificing some spatial resolution to get the requisite number of picture-elements (`pels') on the screen . When the digital spectrogram is completed (stage C), the first picture processing begins at stage D . In this stage the stored spectrogram is processed to select the pels belonging to the dominant line . Further line enhancement is applied to the resulting sparse spectrogram to remove any residual noise and any traces of the weaker lines . By this stage the approximate position of the dominant line has been extracted. In stage E the position is determined more accurately
200
300
FREQUENCY (Hz)
Fig .
2.
Raw digital spectrogram containing
500 x 200
picture elements.
139
Processing noisy line spectrograms as digital pictures
knowledge of the signal, to design a series of matched filters to detect the weaker members of the line family. In this section we firstly consider the production of the raw spectrogram, then discuss the spectral characteristics of the sources of interest It is then shown how this knowledge can be used to design the required filters providing we can extract the dominant line . Spectrogram generation H
The spectral estimates we use are of the 'periodogram' type 1 '-si that is, after appropriate anti-aliasing filtering the signal is sampled at a rate f, until a buffer of size N, is filled. The acquisition time for this is T where
G
REASSEMBLE LINE L
THRESHOLD OUTPUT
PROPOSED PICTURE PROCESSING SCHEME
Fig. 3 . The proposed analysis scheme .
(1) Denoting the digitised signal * by r and the discrete Fourier transform (DFT) of this by R then these are related by
and any gaps in the line are filled in . This results in a template of the dominant line and this can be used (stages F and G) to detect other members of the family by a matched filter process. Thus by stage G we have decided what lines are present but in so doing have lost their temporal structure . As this time variation may be of interest in the analysis of source behaviour it is desirable to restore some of this information to the final `clean' spectrogram. This is done in stage H. The scheme outlined above has been implemented on a Data General Supernova mini-computer with 28K of memory and disc drive using a mixture of FORTRAN and assembly language . At present the FFT and power spectrum computation arc done by software, and for higher sampling rates (> 1 kHz) some data points are lost if the processing is done on-line . As a temporary measure to overcome this, prior to the installation of a hardware FFT processor, the analog signal is first digitised, then stored on magnetic tape and this is used as the input to the processing . As a guide to the total times involved in this processing the measured times for the example given in the text are shown as Table 1 . These times include both terminal and CPU operations .
DIGITAL SPECTROGRAM ANALYSIS
The detection method we are proposing is based on using the extracted dominant line, and our a priori
1 N,-i i V, R, =- ~ r4 exp - A', x-o Ns .
where the frequency has been normalised to units of l/T (`bin-widths') . The time series r is real so only Ns/2 independent points are produced by (2). Thus for each spectral estimate N = NJ2 real and imaginary points are produced by the DFT processor . As the power spectral estimate P is given by P„=R,Rn * ,
* No windowing in the time domain is performed in view of the subsequent two-dimensional filtering of the spectra . Also, non-uniform windowing reduces the effective resolution of the analysis. For certain applications the linewidth is important information . I
Stages
Processing
Time (min)
A-C
Digitise analog signal compute power spectrum, store spectrogram Extract dominant line design template Apply filters Threshold and re-assemble
10
F G H
Total time : 52 min .
(3)
then for each estimate we have N real points to store . If there are M such estimates per spectrogram and we wish to retain L bits of intensity information then LMN bits of storage are needed . It should be noted that in displayed spectrograms, such as Fig . 2, only 3 bits of intensity or 'grey-level' are shown . The full 8 (or 16) hits are available for computation, however . In our presentation we plot time along the y-axis and frequency along the v-axis, the transpose of some other displays (e.g . speech spectrograms).
Table
D-E
())
H
25 9
140
D. G .
Line profiles To design the required family, the spectral 'profiles' of the line structure must be known together with their time variation . The profiles are considered in this section and the time variation in the next. Spectrograms such as Fig . 2 arise from a process of the form
r(t) = s(t) + n(t),
(4)
where n(t) is a white noise process and the signal s(t) is given by K
s(t) _ Z ak exp i2nbka(t)t, k=1
(5)
where the b's are constants, such that b k is the 'order' of the k-th quasi-harmonic and the a's and a (t) are functions of time. The a's correspond to the 'line strengths' and we shall treat these as being constant over one acquisition period T a (t) is the time-varying 'centre frequency' of the fundamental line . The centre frequency of the k-th quasi-harmonic is b k a (t) . For convenience we shall normalise these frequencies to bin units (i .e . IIT) and the time to units of T We suppose that over one acquisition period, i .e. the interval a (t) may be approximated by a straight line. Thus over this interval
a (t) = (fo + fit),
K
Designing the test fine family The template family we are designing are to be applied in the frequency-time plane . In the time domain the signal and noise are related by equation (4), Each 'row' of the spectrogram is thus described by a function such as :
p(j)
=
IS(J)12
+ S(f)N*(f) +
S*(f)N(f)
+ IN(f)h, (10)
where S(f) and N(f) are the short-term estimates of the Fourier transforms of s(t) and n(t) respectively . For the two-dimensional spectrogram this equation becomes of the form
p(ft) =g(f,t)+nAf,t)+h(ft),
#
f
of these 6 harmonics are plotted below each other ; the actual frequency separation of their midpoints depends on fo . The continuous case is shown in the left-hand column and it can be seen that the 'width' of the k-th harmonic profile is proportional to the shift bk whereas the 'height' falls off with increasing shift . In fact, for very high shifts the profile approaches a rectangle!" The other columns show the result of evaluating (8) for the discrete case. The output depends on whether fo is an exact multiple of 1/T or not . In the centre column we assume it is and in the right-hand case assume the other extreme, where fo is one and a half multiples of 1/T, is shown. Actual profiles tend to fall in between these extremes as can be seen in Fig . 1 . In this figure some of the curious 'splitting' which appears for f = 3 in Fig. 4 is also evident .
(6)
where fo and (t are constants. Here fo is the fundamental frequency at t = 0 and /3 is the total shift in frequency over this time . The output of the Fourier processor (at t = }) is S(f)= ~ak k= 1
NiCHOL
` expi2nb k(fo +ft)texp(-i2rft)dt,
(7) assuming for the time being that S(f) is continuous. Now, completing the square and changing the integration variable in (7) leads to : K
S(f)= Y, ak/2(fb k)# k-1 rI x exp(-)nd,/2fbk) J
skllRbk)++(ftbk)* 4100 1 -(IIbk1 1
x exp(irzz''/2)dz,
(8)
where dk
=
NA - f
(9)
The integral in (8) is the Fresnel Integral and is well known from diffraction theory . It must be evaluated numerically, Gersh and Kennedy (6) have solved a similar equation to (8) for a wide range of shifts in describing the spectrum of gliding tones . For our purposes the solutions shown in Fig . 4 provide sufficient illustration . For this case we have K = 6, f3 = 1, b k = k and a, = a 2 . . . a b = 1 . In Fig . 4 the profiles* * Actually we have plotted the estimate of the power spectrum, i .e. p(f) = S(f)S*(f).
(11)
where g(f,t), nAf,t) and h(ft) are the two-dimensional signal power, cross-spectral and noise power terms respectively. One method we might use to detect the signal g(ft) in the spectrogram is by a correlation detector such as template matching which is a form of matched filtering (7) . This seems a reasonable approach as the signal and noise are uncorrelated . In the present case it is known that g(ft) is reps sented by a family of K lines s1, 325 s 3 . . . . sK. Knowing one of these enables us to design templates representing all other (possible) members of the family, as the time variation of the centre frequency of each other member is given by multiplying the known centre frequency, say a'(t), by a constant (the 'order factor) and the profile can be found in the previous Section . Denote the line corresponding to the order factors b by sb . Now the actual line family consists of K lines and is represented by :
Of t)
_ k=1 1: sb,, ,
( 12)
but the order factors b k and the number K are unknown . To find these a template corresponding to each possible family member is designed and applied to
141
Processing noisy line spectrograms as digital pictures B
A CONTINUOUS
C DISCRETE
HARMONIC N,o 1 l
2
2
2
3
3
3
4
4
4
5
5
5
6
6
6
POWER SPECTRA OF NON-STATIONARY HARMONIC FAMILY
Fig . 4 . The power spectra produced by short-term spectral analysis of a linearly sweeping family of harmonics. the noisy picture (equation (11)) . For the line of order b the template matching gives output zs = IT sb (ft) P(P),
(13)
where the summation is performed over the whole picture . The most probable members of the family are those corresponding to high values of (13) . The scheme we shall use to choose the most likely family is thus : (a) Find the dominant line . (b) Design test templates .sb(f,t), for each other possible line component. (c) Evaluate zb for each such line . Select the test template as representing a family member if z b is significantly large. Reject otherwise. DOMINANT TEMPLATE EXTRACTION
Dominant line detection In this section we briefly describe the method used to extract the dominant line from the spectrogram .
This problem is not so trivial as it might seem from an inspection of, say, Fig . 2. Humans are much better at this type of processing than are serial machines . The basic problem is that the number of possible pathways through even a moderate size picture soon becomes far too large for any exhaustive sequential search routine!' ) The method employed in the present scheme is based on two constraints which make the problem more tractable. These are : (1) The position of the line is a (strict) function of time, i .e. for each time-sample there is only one centre frequency . In the previous sections this was denoted by J = u'(t) . (2) The dominant line is several dB above the noise . The first of these constraints means that the line cannot double back or branch and the second means we need only look at the stronger points on the spectrogram when searching for pathways . However we do allow some sudden shifts in the dominant line which cause gaps to form. Many line tracking tech-
142
D . G . NICHOL
0
TIME
(sets)
100
200 0
100
200
300
400
FREQUENCY (Hz)
Fig. 5(a). Thresholded 'sparse' spectrogram .
0 TIME
(sets)
100
I 200 0
I 100
I I
200
300
400
500
400
500
FREQUENCY (Hz)
(b). Line-enhanced sparse spectrogram .
0
TIME
(sets) 100
200 0
100
200
300
FREQUENCY (Hz)
(c). Template of test-family.
Processing noisy line spectrograms as digital pictures niques101 tend to be badly behaved at such discontinuities, e .g . that at t = 130 sec in Fig . 2 . Figure 5 shows the results at three stages in extracting the dominant line from Fig . 2 . The third stage is discussed in the next Section but the first two stages are as follows.
distribution, as many chords have similar slopes, whereas chords for off-line points (e .g . Q) tend to be more randomly distributed . Figure 6 shows the distribution functions for points a and flTo produce these distribution functions from Fig . 5(a) we use only non-zero points in a local region, R, 0, , ( of 20 x 20 pets in this case), surrounding each non-zero point (tit, n). The test we use to distinguish distributions corresponding to noise pets from those corresponding to line pets, is whether in a given distribution of say K, angle bins, one bin has more than d % of N n,,,,, the maximum possible points. If not, the pet is rejected . By the maximum possible, we mean the case where the dominant line is straight and passes diagonally through the region R . ; thus N m ,, x equals the numbers of pets on the diagonal . To produce Fig. 5(b) we put K„ = 6 and i. = 30°x. These numbers, which were found by experiment, have proved satisfactory for all spectrograms from this source. Other sources may require different values, however. Some competing line structure is still visible in Fig . 5(b). Before proceeding to the next stage this can mostly be removed by using the first constraint listed above, i .e . that there is only one centre frequency per picture row . Thus as the competing structure is weaker it usually can be removed simply by accepting only the pet in each row which has the strongest line transform value according to Moore's
Stage 1
Take the stored spectrogram and form a new one by selecting the K x maximum pets from each row . The value of K, selected should increase with the noise level but if it is too high the subsequent processing is unnecessarily long as too many noise pels will have been included . As the dominant line in Fig . 2 is very strong the choice of K,, = 2, as shown in Fig . 5(a), is adequate to retain most dominant line points . There are still gaps and noise points, however . Stage _' To remove the noise points we take the thresholded picture.. such as Fig . 5(a) . and compute the local angular distribution function for each non-zero pet . Mooret 101 and Moore and Parkert i ' 1 discuss extensively the use of this function which is simply the distribution of the slopes 0 (where 0 'C 6 < n) of the chords joining the point in question to all other nonzero points . Basically the idea is that points on the line structure, such as x in Fig. 5(a), will have 'peaky'
1
2
143
3
4
ANGLE BIN NUMBER
20
10
0
Fig . 6 . Angular distribution function for points x and (J of Fig . 5(a) .
144
D. G.
algorithm.") A technique* based on `distance weighting' (Nichol")) is used to remove the remainder . It should be noted that the techniques described above can actually be used for much noisier situations than illustrated here and have been designed with such cases in mind . Figure 5(b) agreeswell with human-eye performance and we can now proceed to `fill in' the gaps . Gap filling and
interpolation
The steps used to produce Fig . 5(c) from 5(b) were : Step I
Find gaps by searching for blank rows in the sparse spectrogram (Fig. 5b) . Find the non-zero end-points of the gap . Denote these endpoints by (m, m,) and (rnv
n2) .
Step 2 As a first approximation connect (m,,n,) and (M2, n2) by a straight line . Find the intersection of this line with rows n, + 1 to n 2 - 1 . Denote the bincoordinates of the intersections by Step 3 Even with additive noise the bin containing the centre point of the dominant line is likely to be a local maximum in its row . Check each intersection . . .,x, 2 _, to make sure it is. If not assign the intersection value to the closest local-maximum bin in the same row . Repeat steps 1-3 until all time gaps are filled. Finally, we attempt to determine more accurately the frequency of the dominant line by the use of interpolation on the raw spectrogram (Fig. 2). A five-point parabolic procedure was found to give satisfactory results. The template Y'(t) is now ready for use . This interpolation is desirable to reduce errors in higher order templates . STRUCTURE DETECTION AND DATA COMPRESSION Line detection
Ideally the family of possible components can now be designed from the stored dominant template a'(t). This procedure, if carried out exactly, would be very time-consuming even if a table look-up scheme for choosing the correct profiles were implemented . In practice it has been found satisfactory to regard all profiles as rectangular, with the width increasing linearly, and the height decreasing linearly, with frequency. This approximation is good for high shiftstsl and as the higher harmonics (and thus higher shifts) tend to be weaker, this is satisfactory from a detection viewpoint. * This method uses a clustering algorithm to associate stronger segments of the dominant line. The weaker segments are accepted or rejected according to their `distance' from the stronger segments .
NicHoL
The test family design is thus made much easier (and quicker) because the width of the component is simply the change in shift of a'(t) multiplied by the order factor (the b,'s) . The order factors are chosen so that (a) the resulting component is entirely in the frequency-time picture . Thus no attempt is made to detect related lines which drift above a frequency of f,/2 ; (b) components are computed for intervals such that the maximum separation of adjacent members is (exactly) one bin . This results in some overlap but this can be overcome by `thinning' as described below . After computing line components thus, we can evaluate equation (16) for each possible line . Figure 7(a) shows a plot of z, vs order factor b (order with respect to the dominant line) . The overlap mentioned above can be seen as a broadening of the peaks in Fig . 7(a) . To detect the structure this plot undergoes a process which : (a) removes a running mean (20 bins wide in this case), (b) suppresses all points which have a stronger `neighbour' within a specified range (± 5 bins) (this is an elementary thinning operationt 13,t4)) (c) thresholds the thinned output of the above . Figure 7(b) shows the plot produced from processing Fig. 7(a) thus . It is clear that all lines detected are harmonics of a sub-harmonic of the dominant line, This is a good test of the process as we have not restricted our search to harmonics but to any related lines . In this case we know that the source is rich in harmonics of the cylinder firing rate and only such lines have been detected . Data
compression
For some purposes an output such as a Fig . 7(b) is all that is required . It is certainly a highly compressed (300 :1) form of Fig . 2. For other purposes we may wish to restore some time information . For example, we may suspect some vibration of commencing when the fundamental reaches some frequency threshold. Figures 8 and 9 show two different levels of time restoration, In Fig . 8 we show the intensity of each line (as detected in Fig . 7b) as the grey level on a time harmonic order plot . No information about frequency variations is retained, however . In Fig . 9 we plot the selected line family out in full and for visual display this is probably the best method . In comparing Fig . 2 and 9 we note that, as mentioned above, there appear to be no 'false alarms' detected whereas the background noise has been eliminated and so have been several non-related lines visible in Fig . 2 . (These may be more easily seen by looking along Fig. 2 from an oblique angle .) This method would thus appear to have potential in removing competing structure as well as random noise . Note that not all time variation is restored in Fig . 9, some smoothing of position and intensity having been carried out (over 3 time samples) . For com-
Processing noisy line spectrograms as digital pictures
145
Fig. 7(a). Output of test matched filter family applied to raw spectrogram .
br (linear units)
I
A
I
I 1 1 2 3 4 5
6
ORDER FACTOR (b). De-meaned and thinned filter output .
Fig. R. Order plot of detected line showing time structure .
7
1 46
D. G . NicuoL 0 TIME (sees)
100
200 0
100
200 300 FREQUENCY (Hz)
400
500
Fig . 9 . Frequency restored plot of detected lines (compare with Fig . 2).
puter analysis there is no need to restore frequency information to each harmonic . We can simply store this for one line and store the intensity values for individual harmonics. If we store Fig . 9 in this fashion, then the compression achieved from Fig . 2 is in the ratio of 75 :1 . Further compression would result if we smoothed over more time samples : the choice dependss on the particular application .
CONCLUDING REMARKS
We have not yet carried out a detailed error analysis for the above process . Under some circumstances, the matched template approach is optimum from a signal-to-noise ratio viewpoint . However errors in measuring the template position, as well as the approximation used in determining the profiles for each test component will tend to reduce the performance . Experimental analysis, using known signals, suggests that the loss from all causes, for low-order harmonics is ca . 2 dB . That is, for the low to medium order harmonics the threshold to be set in the line detection stage is some 2-2.5 dB higher than it should theoretically be if the exact line position and profile were known . As position errors are more important for higher harmonics, this will undoubtedly increase withh order . In spite of these losses it is believed that the proposed scheme produces very useful results and gives a significant improvement on unaided human operator performance,
REFERENCES
S . Watanabe, Automatic feature extraction in pattern recognition, in Automatic Interpretation and Classification of Images, (A. Grasselli, Ed.), pp. 131-136 . Academic Press, New York (1967) . 2 . W. S . Meisel, Computer-Oriented Approaches to Pattern Recognition, pp. 24-25 . Academic Press, New York (1972). 3 . J . B . Thomas, An Introduction to Statistical Communication Theory, pp . 96-97. Wiley, New York (1969). 4 . J . W . Tukey, An introduction to the calculations of numerical spectrum analysis, in Spectral Analysis of Time Series ; (B . Harris, Ed.), pp . 25-26 . Wiley, New York (1969) . 5 . P . J . Daniell, Symposium on autocorrelation in timeseries, J. R . Statist . Soc. B 8, 88-90 (1946) . 6. W . Gersch and J . M . Kennedy, Spectral measurements of sliding tones, IRE International Convention Record, Part 2, pp . 157-170 (1960). 7 . H. C . Andrews, Automatic interpretation and classification of images by use of the Fourier domain, in Automatic Interpretation and Classification of Images, (A. Grasselli, Ed .), pp . 187-198. Academic Press, New York (1972) . 8. U . Montanan, On the optimal detection of curves in noisy pictures, Comm . ACM 14, 335-345 (1971) . 9 . J . R . Williams and G . G . Ricker, Digital line trackers, IEC Technical Note 4800-606 (1971). 10 . D . J . Moore, An approach to the analysis and extraction of pattern features using integral geometry, IEE Trans . Man Syst . Cybern . 2, 97-102 (1972). 11 . D . J . H . Moore and D . J . Parker, Analysis of global pattern features, Pattern Recognition 6, 149-164 (1974). 12 . D . G . Nichol, Separating competing line signatures in digital spectrograms, Tech Note 1624, Australian Defence Scientific Service (1977). 13 . A. Rosenfield and M. Thurston, Edge and curve detection for visual scene analysis, IEEE Trans . Comput. 20, 562-569 (1971). 14. D. G . Nichol, The processing of bathythermograph data : a picture analysis approach, Pattern Recognition 8, 209-218 (1976). . 1
About the Author-DAVID GAwma NidnoL received the B .Sc. (Hons) degree and Ph.D . degree in space physics from the University of Tasmania in 1967 and 1972 respectively . Since 1972 he has worked on picture processing and pattern recognition problems which arise in sonar research at the Weapons Research Establishment of the Australian Department of Defence . Dr . Nichol is a member of the Australian Institute of Physics .
Pergamon Press 1977 . Vol. 9 . pp. 137-146 . Printed in Great Britain
PROCESSING NOISY LINE SPECTROGRAMS AS DIGITAL PICTURES D. G . Niciiot . Department of Defence, Australian Defence Scientific Service, Weapons Research Establishment, G .P.O . Box 2151, Adelaide, S . Australia (Received 22 December 1976 ; in revised form 13 April 1977)
Abstract-Contiguous short-term spectral estimates of non-stationary acoustic signals are often displayed as intensity-modulated frequency-time pictures ('spectrograms') . Some of the problems and advantages in storing and processing certain types of spectrograms as digital pictures are explored . The particular type of spectrogram discussed is from a source producing a time-variant 'line' spectrum consisting of one dominant primary line and several much weaker harmonic or quasi-harmonic lines. This type of spectrogram is commonly produced by rotating machinery . It is supposed that this signal is measured amidst additive broad band noise . The resulting spectrogram consists of a family of related meandering lines, the weaker members of which tend to be lost in the noise . In the proposed picture analysis scheme the raw spectral estimates are stored as a digital spectrogram. This spectrogram is then processed by line enhancement methods to extract the primary line to form a 'family-template' . Using this template the weaker lines can be detected and extracted . By re-assembling the extracted components, a cleaned-up low-entropy frequency-time picture of the signal is produced which shows more structure and which is much more amenable to automatic analysis than the original spectrogram .
Watanabe,IU Meisel1 `o) that high dimensionality of the data is the greatest barrier to successful computer pattern recognition . (e .g.
INTRODUCTION
Overview and summary
The production and analysis of spectrograms (that, is, intensity-modulated frequency-time pictures) is common in many technical studies . Examples of such studies include speech research, vibrational analysis and radio astronomy . The present work proposes a picture-processing technique for the enhancement and data compression of a class of spectrograms containing 'line-structure' as defined in the next section . The type of spectrograms discussed is produced by analysis of acoustic emissions from rotational machinery, and as such is particularly important in the field of automotive research. Here, such spectrograms are used for such purposes as engine diagnosis, tuning and noise suppression . The aim of the picture processing is firstly to detect structure in the spectrogram which is too weak to be visible to the naked eye- This is done by reducing the background noise and enhancing the signal structure . Reducing the noise automatically reduces the entropy of the spectrogram because the noise is `bad information.' (This is possible because we have some knowledge of the signal which is external to the given spectrogram and can use this to distinguish signal from noise .) This partially achieves our second picture processing aim of data reduction . Further data compression can be achieved by removing some of the redundancies in the signal-structure. By such means the noisy high-entropy input spectrogram is transformed into a high-quality low-entropy picture which is much easier to analyse . This data-reducing transformation is particularly important if this analysis is to be done by computer, as it is generally recognised
Line spectrograms
Line spectrograms are produced by sources having 'line spectra', that is, sharp narrow peaks in a onedimensional power-frequency plot such as Fig. I . When plotted as spectrograms, these peaks appear as two-dimensional lines in the frequency-time plane as in Fig. 2. If the signals were stationary then the lines would be straight and parallel to the time axis . Clearly this is not the case in Fig . 2 which are contiguous short term spectra of noise from an idling automobile engine to which band limited white noise (30--800 Hz) has been added . As the engine speed drifts up so do nearly all the lines in Fig. 2. This is because practically all lines in Fig. 2 are multiples of a fundamental frequency corresponding to the cylinder firing rate. We shall use the term 'related lines' to include all lines which are multiples (not necessarily integer multiples) of some fundamentalWe include non-integer harmonics ('quasi-harmonics') because certain components, e .g . noise arising from gearing reduction, can be of this form . In the examples discussed only integer harmonics are present but for generality non-integer lines are still looked for . Now if the background noise is sufficiently high the weaker members of a relatedd family can be lost and, because of the changing frequencies . traditional detection methods (such as time integration) do not retrieve them . Other detection methods such as `order tracking', which is offered as an option on various commercial analysis systems, rely on exact knowledge of the changes in fundamental frequency . 137
1 38
D . G. NtCHOL
LINEAR UNITS
100
200
300
400
500
FREQUENCY (Hz)
Fig . I . Raw spectral estimates from non-stationary source .
This paper describes a scheme for enhancing the weak lines without such knowledge . The only restriction is that one member of the family is reasonably strong compared to the background and to the other members . The other lines can be much weaker and breaks in the dominant line can be tolerated. Implementation A flow chart of the process developed is shown in Fig. 3 . The input signal is sampled by the analogto-digital converters at stage A and the short-term power-spectrum computed by the Fast Fourier Transform (FFT) processor. Succeeding estimates of the spectrum are stored on disc until the spectrogram buffer is filled. Typically the size of this would be
0
100
512 x 300 16-bit words. This buffer is displayed. as in Fig. 2, by a dot matrix plot on a Tektronix 4010 graphics terminal . Only eight grey levels are plotted by the display routine and even this requires sacrificing some spatial resolution to get the requisite number of picture-elements (`pels') on the screen . When the digital spectrogram is completed (stage C), the first picture processing begins at stage D . In this stage the stored spectrogram is processed to select the pels belonging to the dominant line . Further line enhancement is applied to the resulting sparse spectrogram to remove any residual noise and any traces of the weaker lines . By this stage the approximate position of the dominant line has been extracted. In stage E the position is determined more accurately
200
300
FREQUENCY (Hz)
Fig .
2.
Raw digital spectrogram containing
500 x 200
picture elements.
139
Processing noisy line spectrograms as digital pictures
knowledge of the signal, to design a series of matched filters to detect the weaker members of the line family. In this section we firstly consider the production of the raw spectrogram, then discuss the spectral characteristics of the sources of interest It is then shown how this knowledge can be used to design the required filters providing we can extract the dominant line . Spectrogram generation H
The spectral estimates we use are of the 'periodogram' type 1 '-si that is, after appropriate anti-aliasing filtering the signal is sampled at a rate f, until a buffer of size N, is filled. The acquisition time for this is T where
G
REASSEMBLE LINE L
THRESHOLD OUTPUT
PROPOSED PICTURE PROCESSING SCHEME
Fig. 3 . The proposed analysis scheme .
(1) Denoting the digitised signal * by r and the discrete Fourier transform (DFT) of this by R then these are related by
and any gaps in the line are filled in . This results in a template of the dominant line and this can be used (stages F and G) to detect other members of the family by a matched filter process. Thus by stage G we have decided what lines are present but in so doing have lost their temporal structure . As this time variation may be of interest in the analysis of source behaviour it is desirable to restore some of this information to the final `clean' spectrogram. This is done in stage H. The scheme outlined above has been implemented on a Data General Supernova mini-computer with 28K of memory and disc drive using a mixture of FORTRAN and assembly language . At present the FFT and power spectrum computation arc done by software, and for higher sampling rates (> 1 kHz) some data points are lost if the processing is done on-line . As a temporary measure to overcome this, prior to the installation of a hardware FFT processor, the analog signal is first digitised, then stored on magnetic tape and this is used as the input to the processing . As a guide to the total times involved in this processing the measured times for the example given in the text are shown as Table 1 . These times include both terminal and CPU operations .
DIGITAL SPECTROGRAM ANALYSIS
The detection method we are proposing is based on using the extracted dominant line, and our a priori
1 N,-i i V, R, =- ~ r4 exp - A', x-o Ns .
where the frequency has been normalised to units of l/T (`bin-widths') . The time series r is real so only Ns/2 independent points are produced by (2). Thus for each spectral estimate N = NJ2 real and imaginary points are produced by the DFT processor . As the power spectral estimate P is given by P„=R,Rn * ,
* No windowing in the time domain is performed in view of the subsequent two-dimensional filtering of the spectra . Also, non-uniform windowing reduces the effective resolution of the analysis. For certain applications the linewidth is important information . I
Stages
Processing
Time (min)
A-C
Digitise analog signal compute power spectrum, store spectrogram Extract dominant line design template Apply filters Threshold and re-assemble
10
F G H
Total time : 52 min .
(3)
then for each estimate we have N real points to store . If there are M such estimates per spectrogram and we wish to retain L bits of intensity information then LMN bits of storage are needed . It should be noted that in displayed spectrograms, such as Fig . 2, only 3 bits of intensity or 'grey-level' are shown . The full 8 (or 16) hits are available for computation, however . In our presentation we plot time along the y-axis and frequency along the v-axis, the transpose of some other displays (e.g . speech spectrograms).
Table
D-E
())
H
25 9
140
D. G .
Line profiles To design the required family, the spectral 'profiles' of the line structure must be known together with their time variation . The profiles are considered in this section and the time variation in the next. Spectrograms such as Fig . 2 arise from a process of the form
r(t) = s(t) + n(t),
(4)
where n(t) is a white noise process and the signal s(t) is given by K
s(t) _ Z ak exp i2nbka(t)t, k=1
(5)
where the b's are constants, such that b k is the 'order' of the k-th quasi-harmonic and the a's and a (t) are functions of time. The a's correspond to the 'line strengths' and we shall treat these as being constant over one acquisition period T a (t) is the time-varying 'centre frequency' of the fundamental line . The centre frequency of the k-th quasi-harmonic is b k a (t) . For convenience we shall normalise these frequencies to bin units (i .e . IIT) and the time to units of T We suppose that over one acquisition period, i .e. the interval a (t) may be approximated by a straight line. Thus over this interval
a (t) = (fo + fit),
K
Designing the test fine family The template family we are designing are to be applied in the frequency-time plane . In the time domain the signal and noise are related by equation (4), Each 'row' of the spectrogram is thus described by a function such as :
p(j)
=
IS(J)12
+ S(f)N*(f) +
S*(f)N(f)
+ IN(f)h, (10)
where S(f) and N(f) are the short-term estimates of the Fourier transforms of s(t) and n(t) respectively . For the two-dimensional spectrogram this equation becomes of the form
p(ft) =g(f,t)+nAf,t)+h(ft),
#
f
of these 6 harmonics are plotted below each other ; the actual frequency separation of their midpoints depends on fo . The continuous case is shown in the left-hand column and it can be seen that the 'width' of the k-th harmonic profile is proportional to the shift bk whereas the 'height' falls off with increasing shift . In fact, for very high shifts the profile approaches a rectangle!" The other columns show the result of evaluating (8) for the discrete case. The output depends on whether fo is an exact multiple of 1/T or not . In the centre column we assume it is and in the right-hand case assume the other extreme, where fo is one and a half multiples of 1/T, is shown. Actual profiles tend to fall in between these extremes as can be seen in Fig . 1 . In this figure some of the curious 'splitting' which appears for f = 3 in Fig. 4 is also evident .
(6)
where fo and (t are constants. Here fo is the fundamental frequency at t = 0 and /3 is the total shift in frequency over this time . The output of the Fourier processor (at t = }) is S(f)= ~ak k= 1
NiCHOL
` expi2nb k(fo +ft)texp(-i2rft)dt,
(7) assuming for the time being that S(f) is continuous. Now, completing the square and changing the integration variable in (7) leads to : K
S(f)= Y, ak/2(fb k)# k-1 rI x exp(-)nd,/2fbk) J
skllRbk)++(ftbk)* 4100 1 -(IIbk1 1
x exp(irzz''/2)dz,
(8)
where dk
=
NA - f
(9)
The integral in (8) is the Fresnel Integral and is well known from diffraction theory . It must be evaluated numerically, Gersh and Kennedy (6) have solved a similar equation to (8) for a wide range of shifts in describing the spectrum of gliding tones . For our purposes the solutions shown in Fig . 4 provide sufficient illustration . For this case we have K = 6, f3 = 1, b k = k and a, = a 2 . . . a b = 1 . In Fig . 4 the profiles* * Actually we have plotted the estimate of the power spectrum, i .e. p(f) = S(f)S*(f).
(11)
where g(f,t), nAf,t) and h(ft) are the two-dimensional signal power, cross-spectral and noise power terms respectively. One method we might use to detect the signal g(ft) in the spectrogram is by a correlation detector such as template matching which is a form of matched filtering (7) . This seems a reasonable approach as the signal and noise are uncorrelated . In the present case it is known that g(ft) is reps sented by a family of K lines s1, 325 s 3 . . . . sK. Knowing one of these enables us to design templates representing all other (possible) members of the family, as the time variation of the centre frequency of each other member is given by multiplying the known centre frequency, say a'(t), by a constant (the 'order factor) and the profile can be found in the previous Section . Denote the line corresponding to the order factors b by sb . Now the actual line family consists of K lines and is represented by :
Of t)
_ k=1 1: sb,, ,
( 12)
but the order factors b k and the number K are unknown . To find these a template corresponding to each possible family member is designed and applied to
141
Processing noisy line spectrograms as digital pictures B
A CONTINUOUS
C DISCRETE
HARMONIC N,o 1 l
2
2
2
3
3
3
4
4
4
5
5
5
6
6
6
POWER SPECTRA OF NON-STATIONARY HARMONIC FAMILY
Fig . 4 . The power spectra produced by short-term spectral analysis of a linearly sweeping family of harmonics. the noisy picture (equation (11)) . For the line of order b the template matching gives output zs = IT sb (ft) P(P),
(13)
where the summation is performed over the whole picture . The most probable members of the family are those corresponding to high values of (13) . The scheme we shall use to choose the most likely family is thus : (a) Find the dominant line . (b) Design test templates .sb(f,t), for each other possible line component. (c) Evaluate zb for each such line . Select the test template as representing a family member if z b is significantly large. Reject otherwise. DOMINANT TEMPLATE EXTRACTION
Dominant line detection In this section we briefly describe the method used to extract the dominant line from the spectrogram .
This problem is not so trivial as it might seem from an inspection of, say, Fig . 2. Humans are much better at this type of processing than are serial machines . The basic problem is that the number of possible pathways through even a moderate size picture soon becomes far too large for any exhaustive sequential search routine!' ) The method employed in the present scheme is based on two constraints which make the problem more tractable. These are : (1) The position of the line is a (strict) function of time, i .e. for each time-sample there is only one centre frequency . In the previous sections this was denoted by J = u'(t) . (2) The dominant line is several dB above the noise . The first of these constraints means that the line cannot double back or branch and the second means we need only look at the stronger points on the spectrogram when searching for pathways . However we do allow some sudden shifts in the dominant line which cause gaps to form. Many line tracking tech-
142
D . G . NICHOL
0
TIME
(sets)
100
200 0
100
200
300
400
FREQUENCY (Hz)
Fig. 5(a). Thresholded 'sparse' spectrogram .
0 TIME
(sets)
100
I 200 0
I 100
I I
200
300
400
500
400
500
FREQUENCY (Hz)
(b). Line-enhanced sparse spectrogram .
0
TIME
(sets) 100
200 0
100
200
300
FREQUENCY (Hz)
(c). Template of test-family.
Processing noisy line spectrograms as digital pictures niques101 tend to be badly behaved at such discontinuities, e .g . that at t = 130 sec in Fig . 2 . Figure 5 shows the results at three stages in extracting the dominant line from Fig . 2 . The third stage is discussed in the next Section but the first two stages are as follows.
distribution, as many chords have similar slopes, whereas chords for off-line points (e .g . Q) tend to be more randomly distributed . Figure 6 shows the distribution functions for points a and flTo produce these distribution functions from Fig . 5(a) we use only non-zero points in a local region, R, 0, , ( of 20 x 20 pets in this case), surrounding each non-zero point (tit, n). The test we use to distinguish distributions corresponding to noise pets from those corresponding to line pets, is whether in a given distribution of say K, angle bins, one bin has more than d % of N n,,,,, the maximum possible points. If not, the pet is rejected . By the maximum possible, we mean the case where the dominant line is straight and passes diagonally through the region R . ; thus N m ,, x equals the numbers of pets on the diagonal . To produce Fig. 5(b) we put K„ = 6 and i. = 30°x. These numbers, which were found by experiment, have proved satisfactory for all spectrograms from this source. Other sources may require different values, however. Some competing line structure is still visible in Fig . 5(b). Before proceeding to the next stage this can mostly be removed by using the first constraint listed above, i .e . that there is only one centre frequency per picture row . Thus as the competing structure is weaker it usually can be removed simply by accepting only the pet in each row which has the strongest line transform value according to Moore's
Stage 1
Take the stored spectrogram and form a new one by selecting the K x maximum pets from each row . The value of K, selected should increase with the noise level but if it is too high the subsequent processing is unnecessarily long as too many noise pels will have been included . As the dominant line in Fig . 2 is very strong the choice of K,, = 2, as shown in Fig . 5(a), is adequate to retain most dominant line points . There are still gaps and noise points, however . Stage _' To remove the noise points we take the thresholded picture.. such as Fig . 5(a) . and compute the local angular distribution function for each non-zero pet . Mooret 101 and Moore and Parkert i ' 1 discuss extensively the use of this function which is simply the distribution of the slopes 0 (where 0 'C 6 < n) of the chords joining the point in question to all other nonzero points . Basically the idea is that points on the line structure, such as x in Fig. 5(a), will have 'peaky'
1
2
143
3
4
ANGLE BIN NUMBER
20
10
0
Fig . 6 . Angular distribution function for points x and (J of Fig . 5(a) .
144
D. G.
algorithm.") A technique* based on `distance weighting' (Nichol")) is used to remove the remainder . It should be noted that the techniques described above can actually be used for much noisier situations than illustrated here and have been designed with such cases in mind . Figure 5(b) agreeswell with human-eye performance and we can now proceed to `fill in' the gaps . Gap filling and
interpolation
The steps used to produce Fig . 5(c) from 5(b) were : Step I
Find gaps by searching for blank rows in the sparse spectrogram (Fig. 5b) . Find the non-zero end-points of the gap . Denote these endpoints by (m, m,) and (rnv
n2) .
Step 2 As a first approximation connect (m,,n,) and (M2, n2) by a straight line . Find the intersection of this line with rows n, + 1 to n 2 - 1 . Denote the bincoordinates of the intersections by Step 3 Even with additive noise the bin containing the centre point of the dominant line is likely to be a local maximum in its row . Check each intersection . . .,x, 2 _, to make sure it is. If not assign the intersection value to the closest local-maximum bin in the same row . Repeat steps 1-3 until all time gaps are filled. Finally, we attempt to determine more accurately the frequency of the dominant line by the use of interpolation on the raw spectrogram (Fig. 2). A five-point parabolic procedure was found to give satisfactory results. The template Y'(t) is now ready for use . This interpolation is desirable to reduce errors in higher order templates . STRUCTURE DETECTION AND DATA COMPRESSION Line detection
Ideally the family of possible components can now be designed from the stored dominant template a'(t). This procedure, if carried out exactly, would be very time-consuming even if a table look-up scheme for choosing the correct profiles were implemented . In practice it has been found satisfactory to regard all profiles as rectangular, with the width increasing linearly, and the height decreasing linearly, with frequency. This approximation is good for high shiftstsl and as the higher harmonics (and thus higher shifts) tend to be weaker, this is satisfactory from a detection viewpoint. * This method uses a clustering algorithm to associate stronger segments of the dominant line. The weaker segments are accepted or rejected according to their `distance' from the stronger segments .
NicHoL
The test family design is thus made much easier (and quicker) because the width of the component is simply the change in shift of a'(t) multiplied by the order factor (the b,'s) . The order factors are chosen so that (a) the resulting component is entirely in the frequency-time picture . Thus no attempt is made to detect related lines which drift above a frequency of f,/2 ; (b) components are computed for intervals such that the maximum separation of adjacent members is (exactly) one bin . This results in some overlap but this can be overcome by `thinning' as described below . After computing line components thus, we can evaluate equation (16) for each possible line . Figure 7(a) shows a plot of z, vs order factor b (order with respect to the dominant line) . The overlap mentioned above can be seen as a broadening of the peaks in Fig . 7(a) . To detect the structure this plot undergoes a process which : (a) removes a running mean (20 bins wide in this case), (b) suppresses all points which have a stronger `neighbour' within a specified range (± 5 bins) (this is an elementary thinning operationt 13,t4)) (c) thresholds the thinned output of the above . Figure 7(b) shows the plot produced from processing Fig. 7(a) thus . It is clear that all lines detected are harmonics of a sub-harmonic of the dominant line, This is a good test of the process as we have not restricted our search to harmonics but to any related lines . In this case we know that the source is rich in harmonics of the cylinder firing rate and only such lines have been detected . Data
compression
For some purposes an output such as a Fig . 7(b) is all that is required . It is certainly a highly compressed (300 :1) form of Fig . 2. For other purposes we may wish to restore some time information . For example, we may suspect some vibration of commencing when the fundamental reaches some frequency threshold. Figures 8 and 9 show two different levels of time restoration, In Fig . 8 we show the intensity of each line (as detected in Fig . 7b) as the grey level on a time harmonic order plot . No information about frequency variations is retained, however . In Fig . 9 we plot the selected line family out in full and for visual display this is probably the best method . In comparing Fig . 2 and 9 we note that, as mentioned above, there appear to be no 'false alarms' detected whereas the background noise has been eliminated and so have been several non-related lines visible in Fig . 2 . (These may be more easily seen by looking along Fig. 2 from an oblique angle .) This method would thus appear to have potential in removing competing structure as well as random noise . Note that not all time variation is restored in Fig . 9, some smoothing of position and intensity having been carried out (over 3 time samples) . For com-
Processing noisy line spectrograms as digital pictures
145
Fig. 7(a). Output of test matched filter family applied to raw spectrogram .
br (linear units)
I
A
I
I 1 1 2 3 4 5
6
ORDER FACTOR (b). De-meaned and thinned filter output .
Fig. R. Order plot of detected line showing time structure .
7
1 46
D. G . NicuoL 0 TIME (sees)
100
200 0
100
200 300 FREQUENCY (Hz)
400
500
Fig . 9 . Frequency restored plot of detected lines (compare with Fig . 2).
puter analysis there is no need to restore frequency information to each harmonic . We can simply store this for one line and store the intensity values for individual harmonics. If we store Fig . 9 in this fashion, then the compression achieved from Fig . 2 is in the ratio of 75 :1 . Further compression would result if we smoothed over more time samples : the choice dependss on the particular application .
CONCLUDING REMARKS
We have not yet carried out a detailed error analysis for the above process . Under some circumstances, the matched template approach is optimum from a signal-to-noise ratio viewpoint . However errors in measuring the template position, as well as the approximation used in determining the profiles for each test component will tend to reduce the performance . Experimental analysis, using known signals, suggests that the loss from all causes, for low-order harmonics is ca . 2 dB . That is, for the low to medium order harmonics the threshold to be set in the line detection stage is some 2-2.5 dB higher than it should theoretically be if the exact line position and profile were known . As position errors are more important for higher harmonics, this will undoubtedly increase withh order . In spite of these losses it is believed that the proposed scheme produces very useful results and gives a significant improvement on unaided human operator performance,
REFERENCES
S . Watanabe, Automatic feature extraction in pattern recognition, in Automatic Interpretation and Classification of Images, (A. Grasselli, Ed.), pp. 131-136 . Academic Press, New York (1967) . 2 . W. S . Meisel, Computer-Oriented Approaches to Pattern Recognition, pp. 24-25 . Academic Press, New York (1972). 3 . J . B . Thomas, An Introduction to Statistical Communication Theory, pp . 96-97. Wiley, New York (1969). 4 . J . W . Tukey, An introduction to the calculations of numerical spectrum analysis, in Spectral Analysis of Time Series ; (B . Harris, Ed.), pp . 25-26 . Wiley, New York (1969) . 5 . P . J . Daniell, Symposium on autocorrelation in timeseries, J. R . Statist . Soc. B 8, 88-90 (1946) . 6. W . Gersch and J . M . Kennedy, Spectral measurements of sliding tones, IRE International Convention Record, Part 2, pp . 157-170 (1960). 7 . H. C . Andrews, Automatic interpretation and classification of images by use of the Fourier domain, in Automatic Interpretation and Classification of Images, (A. Grasselli, Ed .), pp . 187-198. Academic Press, New York (1972) . 8. U . Montanan, On the optimal detection of curves in noisy pictures, Comm . ACM 14, 335-345 (1971) . 9 . J . R . Williams and G . G . Ricker, Digital line trackers, IEC Technical Note 4800-606 (1971). 10 . D . J . Moore, An approach to the analysis and extraction of pattern features using integral geometry, IEE Trans . Man Syst . Cybern . 2, 97-102 (1972). 11 . D . J . H . Moore and D . J . Parker, Analysis of global pattern features, Pattern Recognition 6, 149-164 (1974). 12 . D . G . Nichol, Separating competing line signatures in digital spectrograms, Tech Note 1624, Australian Defence Scientific Service (1977). 13 . A. Rosenfield and M. Thurston, Edge and curve detection for visual scene analysis, IEEE Trans . Comput. 20, 562-569 (1971). 14. D. G . Nichol, The processing of bathythermograph data : a picture analysis approach, Pattern Recognition 8, 209-218 (1976). . 1
About the Author-DAVID GAwma NidnoL received the B .Sc. (Hons) degree and Ph.D . degree in space physics from the University of Tasmania in 1967 and 1972 respectively . Since 1972 he has worked on picture processing and pattern recognition problems which arise in sonar research at the Weapons Research Establishment of the Australian Department of Defence . Dr . Nichol is a member of the Australian Institute of Physics .