Dual Channel Fft Analysis - Part1 676c15

PREVIOUSLY ISSUED NUMBERS OF BRÜEL & KJÆR TECHNICAL REVIEW

4-1983 Sound Level Meters - The Atlantic Divide Design principles for Integrating Sound Level Meters 3-1983 Fourier Analysis of Surface Roughness 2-1983 System Analysis and Time Delay Spectrometry (Part II) 1-1983 System Analysis and Time Delay Spectrometry (Part I) 4-1982 Sound Intensity (Part II Instrumentation and Applications) Flutter Compensation of Tape Recorded Signals for Narrow Band Analysis 3-1982 Sound Intensity (Part I Theory). 2-1982 Thermal Comfort. 1-1982 Human Body Vibration Exposure and its Measurement. 4-1981 Low Frequency Calibration of Acoustical Measurement Systems. Calibration and Standards. Vibration and Shock Measurements. 3-1981 Cepstrum Analysis. 2-1981 Acoustic Emission Source Location in Theory and in Practice. 1-1981 The Fundamentals of Industrial Balancing Machines and their Applications. 4-1980 Selection and Use of Microphones for Engine and Aircraft Noise Measurements. 3-1980 Power Based Measurements of Sound Insulation. Acoustical Measurement of Auditory Tube Opening. 2-1980 Zoom-FFT. 1-1980 Luminance Contrast Measurement. 4-1979 Prepolarized Condenser Microphones for Measurement Purposes. Impulse Analysis using a Real-Time Digital Filter Analyzer. 3-1979 The Rationale of Dynamic Balancing by Vibration Measurements. Interfacing Level Recorder Type 2306 to a Digital Computer. 2-1979 Acoustic Emission. 1-1979 The Discrete Fourier Transform and FFT Analyzers. 4-1978 Reverberation Process at Low Frequencies. 3-1978 The Enigma of Sound Power Measurements at Low Frequencies. 2-1978 The Application of the Narrow Band Spectrum Analyzer Type 2031 to the Analysis of Transient and Cyclic Phenomena. Measurement of Effective Bandwidth of Filters. 1-1978 Digital Filters and FFT Technique in Real-time Analysis. 4-1977 General Accuracy of Sound Level Meter Measurements. Low Impedance Microphone Calibrator and its Advantages. 3-1977 Condenser Microphones used as Sound Sources. (Continued on cover page 3)

TECHNICAL REVIEW No. 1 — 1984

Contents Dual Channel FFT Analysis (Part I) H. Herlufsen............................................................................................... 3

DUAL CHANNEL FFT ANALYSIS (PART I) by H. Herlufsen, (M.Sc.)

ABSTRACT The first part of this article introduces basic dual channel FFT measurements. The physical interpretation of the Cross Spectrum, which is the fundamental function in these measurements, and the Coherence Function are dealt with in some detail. Two different methods for estimating the complex Frequency Response Function of a system, from the input and the output signals, are derived, and it is shown which of the two estimates, H1(f) and H2(f), should be used in different practical measurement situations. Various excitation techniques for system analysis are described and their advantages and disadvantages for specific applications outlined. A number of practical measurements, using the Brüel & Kjær Dual Channel Signal Analyzer Type 2032/2034, are presented to illustrate the function estimates obtained with the different techniques. Part 2 of this article deals with the applications of the time domain functions, Hilbert Transform and Sound Intensity. SOMMAIRE La première partie de cet article intoduit les mesures de base en analyse FFT deux voies. L'interprétation physique de I'interspectre, qui est une fonction fondamentale dans ces mesures, et la fonction de coherence sont traités avec une certaine profondeur. Deux méthodes différentes pour I'estimation de la réponse en fréquence complexe d'un système, à partir des signaux d'entrée, sont dérivées, et I'on montre quelle estimation de H1(f) ou H2(f) doit être utilisée dans les différentes situations pratiques de mesure. Diverses techniques d'excitation pour I'analyse des systèmes sont décrites en soulignant leurs avantages et inconvénients dans des applications particulières. Plusieurs mesures pratiques, effectuées à I'aide de I'Analyseur FFT à deux voies Type 2032 ou 2034, sont présentées pour illustrer les estimations de fonctions obtenues avec les différentes techniques. La deuxième partie de cet article traite des applications des fonctions temporelles, de la transformée de Hilbert et de l'intensité acoustique. 3

ZUSAMMENFASSUNG Der erste Teil dieses Artikels führt in die Grundlagen von 2-Kanal-FFT-Messungen ein. Die physikalische Interpretation des Kreuzspektrums und der Koherenzfunktion werden in einigen Einzelheiten beschrieben. Zwei verschiedene Methoden zur Bestimmung des komplexen Frequenzgangs aus den Ein- und Ausgangssignalen eines Systems werden genannt und welcher — H1(f) and H2(f) in bestimmten praktischen Meßsituationen angewendet werden sollte. Verschiedene Anregungstechniken für die Systemanalyse sowie ihre Vor- und Nachteile für spezielle Anwendungen werden diskutiert. Es werden mehrere Messungen unter Verwendung des Brüel&Kjær Zweikanal-Signalanalysators 2032 und 2034 dargelegt, um die erreichten Näherungen mit den verschiedenen Techniken zu illustrieren. Der zweite Teil wird Anwendungen von Zeitbereichfunktionen, der Hilbert-Transformation und der Schallintensität umfassen.

1. Introduction The main objective in system analysis is to measure input - output relationships. A dual channel FFT analysis of the input and output of a system (Fig.1) permits calculation of a function which describes its dynamic behaviour, assuming the system is linear. Conversely, this function characterizes the system independent of the signals and can be used to predict the output due to a known input or to calculate the input which will cause a given output. The validity of these results will of course depend upon the quality of this model. An essential application of two channel analysis will therefore be to obtain a measure of the validity of the linear model which forms the basis of the system analysis. The two channel function which can be used for this, "linearity check", is the so-called Coherence Function. The Coherence Function is therefore of fundamental importance for a two channel measurement.

Fig. 1. System with input signal a(t)) and output signal b(t). The Fourier Transform of a(t)) and b(t) are A(f) and B(f) respectively 4

The function which describes the system when presented in the time domain is called the Impulse Response Function h(τ), and when presented in the frequency domain is called the Frequency Response Function H(f). H(f) and h(τ) are related via the Fourier Transform and contain the same information about the system, although in two different domains. For the linear system in Fig.1 the relation between the input a(t) and output b(t) is mathematically described by:

(1.1)

where b(t) is the convolution of h(t) and a(t). In the frequency domain the relation between input and output is given by B(f) = H(f) - A(f) (1.2) where A(f) and B(f) are the Fourier Transforms of a(t) and b(t). Performing system analysis with a Dual Channel FFT Analyzer gives a number of advantages. Firstly, it is possible to measure the Frequency Response Function even though the input signal or the output signal is contaminated with extraneous noise i.e. even in situations where other inputs than the measured are present. Secondly, the estimated Frequency Response Function(s) will represent the best linear fit (in the least squares sense) to the system. This is very important in analysing nonlinear systems where a best linear fit is needed for mathematical modelling. The estimated Frequency Response Function will naturally depend upon the type of signals involved in the analysis and their level. The concepts of Frequency Response Function estimations are dealt with in Sections 4 and 5. In several applications the Dual Channel Analyzer is used to measure two outputs of a system. The time delay or the phase lag between the two signals can be determined from the Cross Correlation Function or the Cross Spectrum which are the basic two channel functions to be defined later. If the two output signals are sound pressure signals measured by two closely spaced microphones, the two channel FFT measurements can be used for calculation of the sound intensity. Transmission path identification and source location are examples of these applications. Another application is analysis of operational mode shapes

5

of structures. This can be done by performing a number of dual output vibration measurements with the operational inputs acting on the system. This application will however not be dealt with in this article. It is the intention in this article to introduce the concepts of two channel FFT analysis in a fairly pictorial form without using a strict mathematical formulation which can be found in the references. 2. Dual Channel Measurements A dual channel measurement can be considered to be a simultaneous measurement of the Cross Spectrum between the signals in the two channels, and a measurement of the Autospectrum (often called Power Spectrum) of each of the two channels. Although several new functions are available with dual channel analysis compared to single channel analysis the only additional calculation which is done in the basic measurement is the calculation of Cross Spectrum between the two signals. All the other functions such as Frequency Response Functions, Coherence, Cross Correlation, Impulse Response etc. are functions which are derived from the three spectra, namely the Autospectrum of the signal in channel A, the Autospectrum of the signal in channel B and the Cross Spectrum between the signals in channel A and B.

Fig. 2. Simplified block diagram of the analyzer in the dual channel spectrum averaging mode 6

Fig.2 shows a simplified block-diagram of the analyzer in the dual channel spectrum averaging mode. Although the analyzer might have other measurement modes such as signal enhancement or amplitude probability, only the measurement mode relevant for normal two channel analysis will be discussed here. The measurement consists of recording the signals in the two channels followed by a Fourier Transform of each of these signals. This will (after multiplication and averaging) result in the two Autospectra and the Cross Spectrum. All the other functions can be computed by post processing and although several of these are needed in the analysis only one measurement has to be performed. This is very important in practical applications. For instance, one estimate of the Frequency Response Function of a system (Fig.1) is computed by dividing the Cross Spectrum between the input and the output signals by the input Autospectrum. The corresponding Impulse Response Function can then be computed from this Frequency Response Function via the inverse Fourier Transform denoted by F-1 in Fig.2. The Coherence Function can be found by dividing the numerical square of the Cross Spectrum by the product of the Autospectra, etc. The function relating the two signals together in the time domain (or delay domain) is the Cross Correlation Function Rab(τ). This function gives a measure of how much the two signals a(t) and b(t) are "alike" with a certain delay τ between them. It can be estimated from the inverse Fourier Transform of the Cross Spectrum. If the signals in the two channels are the sound pressure signals measured by two closely spaced microphones, the sound intensity (not shown in Fig.2) can be calculated from the imaginary part of the Cross Spectrum. All these functions will be dealt with later separately. 2.1. The Autospectra and the Cross Spectrum The two channel analysis is based on the concepts of Fourier Transform, which converts information from the time domain into the frequency domain or vice versa. No information is gained or lost in transforming from one domain to the other. The idea of Fourier analysis is to present the information in such a way that it is easy to interpret and facilitate solutions of the problems. The Fourier Transform of a time signal a(t) defines the complex spectrum A(f) and is given by:

7

(2.1) Likewise we have the spectrum B(f) of the time signal b(t) (2.2) The Autospectra of a(t) and b(t) are then defined by SAA(f) = A* (f) · A(f)

and

SBB(f) = B* (f) · B(f)

respectively,

where * indicates complex conjugation. The Fourier Transform which is performed in practice is the so-called Discrete Fourier Transform (DFT). The Fast Fourier Transform (FFT) is just an algorithm which computes the DFT with a greatly reduced number of arithmetical operations compared to a direct computation. The DFT is a transform that works on finite time records of length T. The time signal is sampled at discrete points in time n · ∆ t, where ∆ t is the sampling time and n is an integer. The DFT results in the Fourier spectrum (given by (2.1) and (2.2)) of the finite time record sampled at discrete frequencies k · ∆ f. k is an integer and ∆ f is the spacing between the calculated lines in the frequency domain. ∆ f is called the resolution and given by the record length T, as ∆ f = 1/T. For a detailed discussion of this subject see Ref. [1], [2] and [3]. Each time history block (i) thus gives an estimate Âi(f) of the Fourier Transform A(f),

(2.3) at the frequencies k · ∆f = k · 1/T. A time weighting function w(t) is often applied to the time signal prior to the calculation of the FFT, but is here omitted to simplify the notation. The expected value of Â*i(f) · Âi(f) is the Autospectrum SAA(f) i.e.

(2.4) 8

In practice averaging is performed over only a finite number of records nd. This will of course introduce an error or an uncertainty in the estimate, which will be discussed in Section 9. Likewise the Autospectrum of the time signal b(t) is given by (2.5) where the individual estimates of B(f) are (2.6) computed at the discrete frequencies f = k · ∆ f = k · 1/T. The Autospectrum is real and gives the distribution of power (or energy) in the signal as a function of frequency. The word power is here used for the mean value of a squared quantity. If the signal is periodic the Autospectrum should be scaled in power i.e. Volt2 or Unit2, (for signals calibrated in a physical unit), as it contains power at some discrete frequencies. Continuous random signals are characterised by a power spectral density (PSD) and the Autospectrum should be scaled in Volt2/Hz or Unit2/Hz. Transient types of signals have energy distributed continuously in frequency and the Autospectrum should thus be scaled in energy spectral density (ESD) Volt2 · s/Hz or Unit2 · s/Hz. The Cross Spectrum SAB(f) between the two signals a(t) and b(t) is defined by A*(f) · B(f), where * indicates complex conjugation. When the DFT is used the Cross Spectrum is found by averaging the individual estimates Â*i(f) · Bi(f) giving (2.7) The comments about the scaling of the Autospectrum will be the same for the Cross Spectrum. It can be shown that Âi(f) and Bi(f) are conjugate even i.e. Âi(-f) = Â*i(f) and Bi(-f) = B*i(f) (see for instance Ref. [2] or [3]). The Autospectra SAA(f) and SBB(f) are therefore real and even and the Cross Spectrum SAB(f) is conjugate even. The information at the negative frequencies is therefore the same as the information at the positive frequencies. This makes it often more convenient to work with the corresponding onesided spectra GAA(f), GBB(f) and GAB(f) defined by 9

The power or energy is then only distributed at non-negative frequencies, see Fig.3. The Cross Spectrum is the fundamental function relating the signals in the two channels. The individual estimates of the Cross Spectrum Â*i(f) · Bi(f) can be written as | Âi(f) | · | Bi(f) ej∆Φi(f) = | Âi(f) | · | Bi(f) | ej(ΦB(f) - ΦA(f)) where | | means numerical value, and ΦA(f) and ΦB(f) are the phase of Âi(f) and Bi(f) respectively. Thus the amplitude of the individual Cross Spectrum estimates is the product of the individual amplitudes and the phase is the phase difference between Bi(f) and Âi(f). The averaging of these individual estimates Â*i(f) · Bi(f) will result in a spectrum SAB(f) having a phase which is a weighted average of the individual estimates ∆Φi(f) (the weighting depending upon the amplitudes of Âi(f) and Bi(f)). The amplitude of SAB(f) will have a value between zero and

depending upon the amount of fluctuation of the phase difference ∆Φi(f) from record to record. This is illustrated in Fig.4. The amplitude of SAB(f) will thus not only depend upon the amplitudes of SAA(f) and SBB(f) but also how much correlation there is between Âi(f) and Bi(f).

Fig. 3. One-sided spectrum G(f) formed from the two-sided spectrum S(f) 10

Fig. 4. Averaging of individual estimates of the Cross Spectrum SAB(f) = A*(f) B(f). a) Phase of each estimate Â*i(f) Bi(f) is the same (b) Some fluctuation in the phase of the estimates Â*i(f) Bi(f) (c) Random phase of the estimates Â*i(f) Bi(f) In the figure the dependency of frequency is left out for convenience

The concept of correlation will be defined and discussed in section 3. Perfect correlation at frequency f will cause ∆Φi(f) to be the same for each estimate (Fig.4.a) whereas if the signals are uncorrelated at frequency f, the estimates ∆Φi(f) will be random and between 0 and 2π causing SAB(f) to be zero (Fig.4.c). Fig.4.b shows a situation where there is some correlation between A(f) and B(f) at the selected frequency. The amplitude of SAB(f) is thus quite difficult to interpret and in itself is not used very often. The phase of SAB(f), however, is the phase difference ("average of") between b(t) and a(t) at frequency f. In an input-output analysis as in Fig.1 the phase of SAB(f) will give the phase response of the system, or in general the phase lag from signal a(t) to signal b(t) as a function of frequency. The Cross Spectrum can also be written as GAB(f) = CAB(f) + jQAB(f), using the one-sided notation. CAB(f) is the part where A(f) and B(f) are in-phase and is called the coincident spectrum or just the cospectrum. The imaginary part QAB(f) is the part where A(f) and B(f) are 90° out-ofphase and is thus called the quadrature spectrum or the quad spectrum. As already mentioned the Cross Spectrum is one of the key functions for calculation of many other two channel functions rather than a function which is used just by itself.

11

2.2. Documentation for Measurement - Measurement Setup Apart from the three basic spectra GAA(f), GBB(f) and GAB(f), the parameter settings on the Analyzer are a vital part of the documentation of a measurement. Fig.5 shows a block diagram of the Brüel & Kjær Dual Channel Signal Analyzers Types 2032 and 2034 for clarification of some of the functions. For the input function it is necessary to know which input channels are used, polarity and sensitivity settings of these, selection of tow- and high- filters, together with calibration settings in Volts/Unit or Units/Volt and which Unit(s) is used for the measurement. For the recording it is relevant to have information about the triggering, delay settings, and selection of frequency span, whether baseband or zoom analysis is used, and for zoom analysis what the centre frequency is. The trigger can be set either on free run or on any of the sources: ch.A, ch.B, external, manual or generator (clock in the Analyzer). The delay between the trigger and recording in channel A must be documented. For systems with inherent delay, (for example acoustic systems), the delay between ch.A and ch.B should be set and documented. Otherwise, systematic bias errors will occur in the analysis. In the digital zoom-processor (see Fig.5) it is possible to select a frequency span (sampling frequency) and a centre frequency, if zoom analysis is performed (as opposed to normal baseband operation). The analysis is performed using the FFT algorithm of the DFT which works on blocks of time data. For several reasons different types of weighting functions can and should be applied to the time history records before the FFT is calculated. Some of these weighting functions will be discussed later. The type of weighting function used in the

Fig. 5. Block diagram of the Brüel & Kjær Dual Channel Signal Analyzer type 2032 or type 2034 12

Fig. 6. Example of a measurement setup from the Brüel & Kjær Dual Channel Signal Analyzer type 2032

analysis is part of the documentation for the measurement. The socalled zero pad mode used for analysis of correlation functions can be considered as being a mode where a special time weighting function, which nullifies the last half part of the record, is used. Also the type of averaging (linear or exponential) and the number of averages should be reported, together with information of how much overlap between the analyzed time records has been used in the analysis, if any overlap conditions have been set. All the above information can be stored as a measurement setup and used as documentation for the measurement. Fig.6 shows an example of such a measurement setup from the Brüel & Kjær Analyzer Type 2032. This measurement setup together with the Autospectra and the Cross Spectrum and perhaps also a time record from each of the channels will make up the total documentation.

2.3. Presentation of Processed Data - Display Setup In the Brüel & Kjær Analyzers Types 2032 and 2034 (see Fig.5) there is a separate processor called the Display Processor or the "Master", which not only controls the Signal Processor to perform the FFT, averaging etc. but also controls all the required post-processing, either after the measurement or during the measurement (live display). This can be division or multiplication of spectra, performing a Hilbert Transform (discussed later) or an inverse FFT (using the Signal processor) etc. Most of the functions are complex and can therefore be presented in a number of formats or coordinates: Real part. Imaginary part, Magnitude, Phase, Imaginary versus Real part (Nyquist plot) or Magnitude 13

versus phase (Nichols Plot). Fig.7 illustrates these different formats for a Frequency Response Function being the transfer accelerance of a mechanical system measured by use of an Impact Hammer, an Accelerometer, and the Brüel & Kjær Analyzer Type 2032. Depending upon the application and what kind of information is required from the measurement, the function, the format, as well as the X and Y scales have to be selected. The scales can be linear or logarithmic and integration or differentiation of the function can also be obtained (if, for instance, mobility or compliance instead of accelerance is desired). This is the documentation for the displayed graph and can be seen as the text in the illustrations on Fig.7. Apart from this so-called display setup there is of course a need for some cursors for reading out of data in different ways. The cursor information is found on the right side of the illustrations on Fig.7.

Fig. 7. Different presentations of a complex function, here the Frequency Response Function of a mechanical system (transfer acceelerance)

14

3. Correlation and Coherence In system analysis where an input and an output is measured it is of great importance, as part of the analysis, to assess the degree of linear relationship between the input and the output. Having two stochastic variables x and y which could be measurements of input and output of a system, the correlation between these variables, is described by the so-called Correlation Coefficient ρ xy defined by:

σ xy ρ xy = σ σ x y

(3.1)

where σ xy is the covariance of x and y given by σ

xy

= E[(x - µx)(y - µy)]

(3.2)

and σ x and σ y are the standard deviations of x and y defined by σx=

and

√E [(x - µx)2 ]

σ y = √E [(y - µy)2 ]

(3.3)

(3.4)

E is the expected value and is found by averaging, while µx and µy are the mean values of x and y respectively, i.e.

and

µx = E [x]

(3.5)

µy = E [y]

(3.6)

ρ xy will have a value between 0 and 1. If there is a perfect linear relationship between x and y i.e. y = αx + β, for all samples of x and y, α and β being constants, ρ xy will be 1. This situation is illustrated in Fig.8a. If the x and y samples are contaminated with some random noise (scatter) as shown in Fig.8.b, ρ xy will be less than 1. Notice the underlying linear relationship between x and . This could be the situation where the system under investigation is perfectly linear but there is some extraneous noise in the measurements (samples) of the input x and/or the output. In a situation where the relation between the samples of x and is well defined but non-linear, the value of ρ xy will also be less than 1, even though the x and y are free of random noise. Fig.8.c shows an example of this situation. 15

Fig. 8. The Correlation Coefficient σxy indicates the amount of linearity between variables x and y If the variables x and y are not at all related to each other, the samples will be randomly scattered as depicted in Fig.8.d and the Correlation Coefficient ρ xy will be zero. Thus, the value of the Correlation Coefficient ρxy is a measure of the degree of linear dependance between the variables (samples) of x and . One of the other functions which can be computed in the two channel 2 FFT Analyzer is the Coherence Function γ (f) defined by

γ2(f) =

GAB(f)2 GAA(f) GBB(f)

=

SAB(f)2 SAA(f) SBB(f)

At each given frequency, f, the Coherence Function corresponds to the Correlation Coefficient Function squared given by

16

ρ2xy =

σ2xy σ2xσ2y

(3.8)

Let x be the complex spectral component at frequency f, A(f), and y be the complex spectral component at the same frequency f, B(f). By comparing (2.4) with (3.3),(2.5) with (3.4) and (2.7) with (3.2) it is seen that the Coherence and the Correlation Coefficient squared are similar (apart from the subtraction of the mean values µx and µy from x and y for σx, σy and σxy and the complex conjugation due to the fact that the spectral components A(f) and B(f) are complex). All the properties of the Correlation Coefficient (squared) will therefore also apply for the Coherence Function. Similar to the Correlation Coefficient, the Coherence γ2(f) of the signals a(t) and b(t) is a function which on a scale from 0 to 1 measures the degree of linear relationship between the two signals at any given frequency f. Without going into a detailed discussion of this subject, only the most common reasons for having Coherence less than one will be discussed very briefly here. Coherence less than one can be due to one or more of the following situations. a) Uncorrelated noise in the measurements of a(t) and/or b(t). b) Non-linearity of the system under investigation. c) Leakage in the analysis (resolution bias error). d) Delays in the system not compensated for in the analysis. a) If the measurements are contaminated with uncorrelated extraneous noise, the individual estimates of the Cross Spectrum Â*i(f)Bi(f) will add as depicted in Fig.9.a, while the Autospectra estimates Â*i(f) · Âi(f) and B*i(f) · Bi(f) being all real are added as in Fig.9.a. The linear related parts of Âi(f) and Bi(f) will add in the Cross Spectrum estimation as shown earlier in Fig.4.a while the uncorrelated noise will gradually average out as illustrated in Fig.4.c. The Autospectra GAA(f) and GBB(f) will include the extraneous noise and the Coherence will be less than one. The more extraneous noise there is in the measurements the lower will the Coherence be. 17

Fig. 9.

Effects on the Autospectral and the Cross Spectrum estimates of: a) extraneous noise in the measurement of a(t) and/or b(t) b) non-linearity in the system The dependency of frequency is omitted for convenience

b) Fig.9.b shows a situation where the gain H(f) (amplitude and phase) of the system depends upon the input level A(f) (related to Fig.1). The sample Â*2(f) · B2(f) has a different phase than the sample Â*1(f) · B1 (f) due to the fact that the input amplitude sample | Â2(f) is greater than the input amplitude sample | Â1(f) |. The argument for low Coherence in such a situation is similar to that where uncorrelated noise was present. In a situation where the phase of H(f) is independent of A(f) but the amplitude of the gain | H(f) | is dependent on the level of A(f), it is easily seen as well that the Coherence will be less than one. Let us for instance assume | Bi(fo) = |Âi(fo)2 for each estimate at frequency fo. It is then found that

18

As a simple example let nd = 3 and |Â1(fo)| = 1, |Â2(fo) = 2, and Â3(fo) = 3 gives γ2(fo) =

(1 + 8 + 27)2 324 = < 1 (1 + 4 + 9) · (1 + 16 + 81) 343

Notice however that the input estimates Âi(f) are assumed to vary in amplitude. Also a non-linear system when excited at one frequency can give rise to response signals at other frequencies. This content of power at other frequencies will be analyzed as extraneous noise at these other frequencies in the output b(t). It is here assumed that the input signal is random causing Âi(f1) to be uncorrelated with Âi(f2), f1 ≠ f2. This will be discussed further in section 5. c) Before discussing why leakage can lead to low Coherence, a brief description of leakage is given. Leakage is a phenomenon which may arise in the frequency domain due to the time limitation of the signal before the FFT calculation is performed. Only a finite time record length of the time signal can be analysed i.e. the signal which is analyzed is the original time history a(t) multiplied by a weighting function w(t). As multiplication in one domain corresponds to a convolution in the other domain (Convolution Theorem for the Fourier Transform, see Ref. [1, 2 and 3]), the spectrum A(f) of the original time signal a(t) will be convolved with the Fourier Transform W(f) of the weighting function w(t). In other words:

giving

19

Fig. 10. Convolution of spectrum A(f), of the original time signal, with the Fourier Transform W(f) of the time weighting function giving the estimated spectrum A(f)

Fig.10 shows how a time weighting w(t) might influence the estimate A(f) of a spectrum A(f). Power in one frequency region leaks into adjacent frequency regions causing the peak amplitude to drop and the amplitudes in valleys to rise as illustrated in Fig.10. The amount of leakage will depend upon the type of weighting function used in the analysis. A more detailed description of leakage and convolution is found in Ref. [1] and [2]. The so-called Hanning weighting, which is one period of a cosine matching the record length of the Analyzer and lifted so it starts and stops at zero, is the normally used weighting function when analyzing continuous random or sinusoidal signals. Use of different weighting functions for different types of signals will be dealt with in Section 5. In explaining why leakage in the analysis might reduce the Coherence, let us consider two situations: 1) The complex Frequency Response Function H(f) = B(f)/A(f) is a constant H in a frequency range wider than the bandwidth of W(f). 2) The complex Frequency Response Function H(f) = B(f)/A(f) undergoes rapid changes in amplitude and phase with frequency. Rapid is here meant relative to the resolution in the analysis (bandwidth of W(f)). This could be the case in the frequency region around a system resonance or anti-resonance. It is assumed that there is no extraneous noise, that the system is linear and that the input signal is a random signal. In the first situation 1) there will be no drop in the Coherence Function as the leakage will be the same in both spectra A(f) and

20

B(f), provided that the same weighting function is used in the two channels. Each of the samples of the spectra Âi(f) and Bi(f) at frequency f can be considered as a sum of the true spectrum A(f) and B(f) and a leakage term Aileakage(f) and Bileakage(f) where B(f) = H(f) · A(f) and Bileakage(f) = H(f) · Aileakage(f) giving

Bi(f) = B(f) + Bileakage(f) = H(f) (A(f) + Aileakage(f)) = H(f) Âi(f)

and Âi*(f) · Bi(f)) = H(f) Âi*(f) · Âi(f) for all the samples in the averaging. Thus | GAB(f)  2 = H(f)2 G2AA = GBB(f) GAA(f) which means that γ2(f) = 1 in the frequency range under consideration. This situation is shown in Fig.11.a. In the other situation 2) where the Frequency Response Function of the system changes rapidly in amplitude and phase (relative to the

Fig. 11

Leakage in the estimates of the spectra A(f) and B(f) a) Frequency Response Function H(f) = B(f)/A(f) is constant (= H) in the frequency range under consideration. The ratio of the leakage Bileakage (f)/Aileakage (f) will be equal to H(f) = H b) Frequency Response Function H(f) = B(f)/A(f) has rapid changes, relative to the resolution in the analysis, in amplitude and phase. Aileakage(f) and Bileakage(f) can vary independant of each other.

21

frequency resolution in the analysis) the leakage will not be related to each other in any simple fashion, see Fig.11.b. The leakage Aileakage(f) and Bileakage(f) can vary to a certain extent independently of each other. This situation thus corresponds to the situation where the input and/or the output signal is contaminated with extraneous noise (Fig.9.a) and there will be a drop in the Coherence Function in this frequency region. How much the Coherence will drop depends of course upon how much the amplitude and phase of the Frequency Response Function changes within the resolution bandwidth in the analysis. The drop in Coherence can also be explained from the time domain description. Having a resonance peak in the Frequency Response Function which is narrower than the analysis bandwidth means that the Impulse Response Function h(τ) of the system is longer than the time record length T. According to the convolution integral (1.1) some of the signal in the time record in channel B will be due to the input signal before the time record in channel A and some of the response signal due to the input signal in the time record in ch. A will be in output signal b(t) after the time record in ch.B. The Coherence therefore drops in the measurement. Decreasing the analysis bandwidth will increase the correlation between the leakage (increase the time record lengh T relative to the length of the Impulse Response Function) and the Coherence will increase. This is illustrated in Fig.12 which shows the magnitude of a Frequency Response Function estimate and Coherence with a resolution of a) 16 Hz and b) 1 Hz. The drop in Coherence at 1040 Hz in a) is thus due to leakage. This phenomenon is also often referred to as resolution bias error. It is important to note that these arguments are valid only if the input signal is random (random from data record to data record). If the signal is deterministic and repeats itself for each data record, leakage will not be detected by the Coherence Function. This can be explained as follows. Having an input signal which repeats itself for each data record will give Âi(f) = Âi(f) and thus Bj(f) = Bj(f), as the leakage is the same for each sample. Each sample of the Cross Spectrum Âi*(f) · Bi(f) will also be identical and the Coherence will be one. Even if the input signal varies in amplitude i.e. Âi(f) = c Âj(f), where c is real, and assuming that the system is linear, we will get Bi(f) = c Bj(f) and thus Âi*(f) · Bi(f) = c2 Âi*(f) · Bj(f). The samples of the Cross Spectrum will add as shown in Fig.4.a and the Coherence will

22

Fig. 12. Magnitude of Frequency Response Function and Coherence Function measured with a resolution of a) 16 Hz and b) 1 Hz 23

be one. This situation occurs in practice when working with an impact hammer as an input force generator for a (linear) mechanical structure. An illustration of this will be given in Section 5 where excitation techniques are discussed. d) Coherence measurement between input and output signals of a system with a physical delay from input to output or with reverberation requires some special attention. Having a situation as described in Fig.13, where there is a delay of τ sec from input to output and a record length (time window) in the analysis of Tsec, the estimated Coherence γ2(f) will be biased by a factor of (1 - τ/T)2 if no delay is set between the data blocks in the analysis. The bias error is derived on the assumption that the signals are white noise (Ref. [4]). In order to get an unbiased estimate of the Coherence Function a delay of #sec has to be set between the time data block in ch.A and the time data block in ch.B. It is thus essential to have a delay setting facility between the channels, when analysis is performed on systems with physical delays. When reverberation (or multiple reflections) is present at the output of the system, the record length T in the analysis should be longer than the reverberation time Trev in the system in order to avoid significant bias errors (Ref. [4]). The reverberation time is defined as the time it takes the output signal to decrease 60 dB after an instantaneous termination of the input signal. A more detailed discussion of some of the subjects concerning Correlation and Coherence is found in Ref. [4].

Fig. 13. Measurements on systems with delays require introduction of this delay between the time data blocks in the analysis

24

4. Frequency Response Function Estimates Probably the most common application of two channel analysis is to measure Frequency Response Functions of physical systems. Referring to the ideal system in Fig.1 the Frequency Response Function defined by B(f) H(f) = , describes the system in the frequency domain. The system A(f) can just as well be described in the time domain by the Impulse Response Function h(τ) defined by the inverse Fourier Transform of H(f) : h(τ) = F-1 {H(f)}. The Impulse Response Function will be dealt with separately in Section 6. A number of assumptions about the system have to be made before the system can be described in of a Frequency Response Function (or an Impulse Response Function): 1. The system must be physically realisable, i.e. it cannot respond to an input before it is applied, or h(τ) = 0 for τ < 0. 2. The system must be time invariant. Its properties may not change with time i.e. h(τ) and H(f) are independent of time. h(τ, t) = h(τ) & H(f, t) = H(f), -∞ < t < ∞. 3. The system must be stable, i.e. it can only respond with a limited amount of energy when excited with a finite amount of energy at the input. This is also true if

4. The system must be linear. This means that if the inputs a1(t) and a2(t) produce the outputs b1(t) and b2(t) respectively, then the input a1(t) + a2(t) must produce the output b1(t) + b2(t), and the input c · a1(t) must give the output c · b1(t), where c is an arbitrary constant. The essence of this is that the functions h(τ) and H(f) characterize the system itself independent of the signals a(t) and b(t) involved. The assumption of linearity is probably the requirement which is most often violated in practical applications. The input signal may have such high amplitude levels that the system will be excited beyond its range of linear behaviour. Also some systems are inherently non-linear and the system function description will not be valid except for very limited ranges of input signal levels. 25

The estimates of Frequency Response Functions using dual channel FFT will give the optimum calculation of H(f) in the least-squares sense (see Ref. [4]). If the input signal is random the best linear approximation to the system is obtained and the Frequency Response Function can therefore be very useful and give meaningful results unless the system is strongly non-linear. Time variations can also cause severe problems in practice. The system characteristics might depend upon temperature, pressure and other parameters giving problems in the linear modelling of systems. However, in some applications, such as speech analysis and synthesis, it is relevant to work with a linear system varying in time i.e., describing the system, being the vocal tract for speech applications, in of a H(f, t). This is of course nonsense in the true Fourier sense (equation (1.2), (2.1) and (2.2)), since the Frequency Domain Function cannot change with time. Having a well defined time window of a given length T the transforms (2.3) and (2.6) could be used to define a Frequency Response Function H(f) from (1.2) for that time window. The next time window (of length T) will then define a new Frequency Response Function for the system via (2.3), (2.6) and (1.2) etc. A typical time base T for speech applications is 25 msec. Analysis of time varying systems will not be dealt with in this article. Ref. [10] and [11] gives some examples of this for speech analysis. The Frequency Response Function of the ideal system of Fig.1, fulfilling the requirements 1) - 4) will now be estimated using the Two Channel FFT Analyzer. The fundamental equation relating the input spectrum and output spectrum is B(f) = H(f) · A(f)

(4.1)

Multiplying by A*(f) on both sides of this equation we get

or

A*(f) · B(f) = H (f) A* (f) A (f)

(4.2)

SAB(f) = H (f) SAA (f)

(4.3)

In practice it is the estimates

of SAB(f) and SAA(f) or GAB(f) and GAA(f) which are measured with the Dual Channel Analyzer. The ratio of these estimates therefore gives a 26

measure of the Frequency Response Function and is called H1(f). Using the one-sided spectra we get from (4.3) H(f) =

GAB(f) ≡ H1(f) GAA(f)

(4.4)

Multiplying (4.1) by B* (f) instead of A*(f) we get

or

B*(f) B(f) = H(f) B*(f) A(f)

(4.5)

SBB(f) = H(f) SBA(f)

(4.6)

Using this equation we see that the ratio of the estimates of SBB(f) and SBA(f) or GBB(f) and GBA(f) also gives a measure of the frequency response function H(f). This estimate is called H2(f) and will in many practical cases give a different result than H1(f). In of the onesided spectra estimates we therefore have G (f) (4.7) H(f) = BB ≡ H2(f) GBA(f) Note that the phase of H1(f) is the same as the phase of H2(f) since

G*BA(f) = GAB(f).

Also note that

H1(f) G (f) GBA(f) = γ2(f) = AB H2(f) GAA(f) GBB(f)

(4.8)

A third method of estimating the Frequency Response Function is found by taking the numerical square of (4.1): B*(f) B(f) = H*(f) H(f) A*(f) A(f) which is

SBB(f) = H(f)2 · SAA(f)

(4.9)

or

GBB(f) = H(f)2 · GAA(f)

(4.10)

The ratio of the output and the input Autospectra GBB(f) and GAA(f) is called the Autospectra estimate and is denoted by | Ha(f) |. As it is real, it contains only information about the amplitude of H(f), unlike H1(f) and H2(f) which contain information of both amplitude and phase. Thus:

| H(f) |2 =

GBB(f) = |Ha(f)|2 GAA(f)

(4.11)

27

The advantages and disadvantages of using the three different methods of determining the Frequency Response Functions in different situations will now be described in of bias errors. In all situations an ideal system will be assumed, and that sufficient amount of averaging is performed in order to remove the random errors. Formulae for random errors will be dealt with in Section 9. In the situations described in Sections 4.1 to 4.5 it is also assumed that sufficiently narrow bandwidth is used in the analysis in order to resolve system resonances (and antiresonances) without any leakage problems (resolution bias). Bias errors due to leakage in the analysis are discussed separately in Section 4.6.

4.1. Ideal Situation The ideal situation is one in which the measured signals a(t) and b(t) are not contaminated by any kind of extraneous noise, and where a(t) is the only input to the system giving the output b(t). This is shown in Fig.14.

Fig. 14. Ideal system in ideal situation

For each individual estimate: Bi(f) = H(f) Âi(f) and therefore from the above section it can be seen that H1(f) = H2(f) = H(f) and Ha(f) = H(f) The Coherence for the measurement will be 1,0, as a(t) and b(t) are ideal and linearly related. In this situation only a few averages are needed in the analysis depending upon the types of signals (type of excitation) involved.

4.2. Noise at output If the measured output signal b(t) is contaminated by extraneous noise n(t) the situation is as shown in Fig.15. The output signal v(t) is linearly related to the measured input a(t) as in the ideal situation. The noise term n(t) is assumed to be uncorrelated to a(t) and therefore not correlated to v(t) either i.e. GAN(f) = GVN(f) = 0. 28

Fig. 15. Ideal system with extraneous noise in measured output signal b(t)

The noise signal n(t) could in a practical situation be noise from transducers, instrumentation or computational noise in the analysis. From Fig.15 we have Bi(f) = Vi(f) + Ni(f) = H(f) Âi(f) + Ni(f) giving

GBB(f) = GVV(f) + GNN(f) = H(f)2 GAA(f) + GNN(f)

and

GAB(f) = GAV + GAN = GAV(f) = H(f) GAA(f)

since

GAN(f) = 0.

The uncorrelated noise is averaged out in the measured GAB(f) as discussed earlier. H1(f) will in this situation therefore give an unbiased estimate of the Frequency Response Function. GAB(f) H1(f) = G (f) = H(f) AA

(4.12)

It can be shown (Ref. [4]) that H1(f) in this situation relates as much of the output b(t) to the input a(t) as possible, minimizing the amount of noise n(t). H1(f) is therefore called the optimum estimate of H(f). For the estimate H2(f) however we get H2(f) =

GBB(f) H(f)2 GAA(f) + GNN(f) = GBA(f) H*(f) GAA(f)

= H(f)

(1 + GNN(f)) GVV(f)

(4.13)

which is an overestimate of H(f) in magnitude due to the noise GNN(f)) measured in GBB(f). The phase of H2(f) is determined by the Cross Spectrum and is therefore correct.

29

The Autospectrum method gives: Ha(f)2 =

GBB(f) GVV(f) + GNN(f) = GAA(f) GAA (f) = |H(f)|2 (1+ GNN(f)) GVV(f)

(4.14)

Ha(f) overestimates H(f) similar to H2(f) by the noise term

The Coherence in this measurement situation is γ2(f) =

or

GAB(f)2 GAA(f) GVV(f) = = GAA(f) GBB(f) GAA(f) GBB(f)

γ2 (f) =

H1 (f) = H2 (f)

GVV(f) GBB(f)

(4.15)

1 GNN(f) (1 + ) GVV(f)

GVV(f) is the part of the measured output signal which is linearly related (through H(f)) to the input signal. From (4.15) it is seen that it is given by GVV(f) = γ2(f) GBB(f). This leads to the definition of the Coherent Power = γ2(f) GBB(f) which is one of the functions which can be computed from the two channel measurement (see Section 2, Fig.2). The uncorrelated noise is given by GNN(f) = (1 - γ2(f)) GBB(f) and the function Non-Coherent Power is therefore defined by Non-Coherent Power = (1 - γ2(f)) GBB(f) With the two channel measurement it is thus possible to calculate the noise spectrum without making a measurement with the signal a(t) "turned off".

30

The signal to noise ratio at the output,

2 GVV(f) is given by γ2 (f)· GNN(f) 1-γ (f)

Another function computed by post processing is therefore

Signal to Noise Ratio =

γ2(f) 1 - γ2 (f)

4.3. Noise at Input Let us consider the situation shown in Fig.16. The noise m(t) is assumed to be uncorrelated with u(t) i.e. GMU(f) = GMB(f) = 0. In practice the noise m(t) could be noise from transducer, instrumentation or computational noise as in the previous situation of Section 4.2. The fundamental relations between the individual estimates are in this situation Bi(f) = H(f) Ûi(f) and Âi(f) = Ûi(f) + Mi(f) giving

GAA(f) = GMM(f) + GUU(f) GBB(f) = | H(f) |2 GUU(f) = | H(f) |2 (GAA(f) - GMM(f))

and

GAB(f) = GUB(f) = H(f) GUU(f)

since m(t) is uncorrelated to u(t). For the H1(f) estimate of H(f) we get

H1(f) =

GAB (f) GAA(f)

=

H(f) GUU(f) GNN(f) + GUU(f)

1 = H(f) 1 + GMM(f) GUU(f)

(4.16)

Fig. 16. Ideal system with extraneous noise in measured input signal a(t) 31

which is a biased estimate of H(f). The more noise GMM(f) there is relative to input signal GUU(f) the lower will the amplitude of H1(f) be. The phase will be correct since it is determined by the Cross Spectrum. The estimate H2(f)) is given by: H2(f) =

GBB(f) |H(f)|2 GUU(f )= = GBA(f) H*(f) GUU(f) H(f)

(4.17)

It is seen that H2(f) is insensitive to uncorrelated noise in the input measurement (apart from a random error due to the insufficient number of averages). This was first pointed out in Ref. [5]. For the Autospectrum method we obtain:

|Ha(f)|2 =

GUU(f) GBB(f) = |H(f)|2 GAA(f) GMM(f) + GUU(f)

=|H(f)|2

1 1 + GMM(f) G (f) UU

(4.18)

The bias error in | Ha(f) | is thus an underestimation by a factor of

This situation with input noise can occur in practice when a shaker is used to excite a mechanical system with stationary force excitation. At the resonance frequencies of the structure the vibration response to a certain force input will be very large (very high mobility). The power ampifier driving the shaker will thus see a very high electrical impedance at the resonance frequencies, causing the current and thereby the input force to drop perhaps as low as the background noise level at those frequencies.

32

The Coherence is given by: |GAB(f)|2 γ2(f) = GAA(f) GBB(f)

=

H 1(f) 1 H2(f) = 1 + GMM(f) GUU(f)

(4.19)

and the input noise is detected by the Coherence Function as expected. Notice that the Coherent Power at the output in this situation is GBB(f). The previously defined Coherent Power = γ2(f) GBB(f) (Section 4.2) has therefore no physical meaning in this situation. Likewise will the previously defined Non-Coherent Power = (1 - γ2(f)) GBB(f) (section 4.2) have no meaning in this situation. The Signal to Noise Ratio,

GUU(f) γ2(f) as in is however given by GMM(f) 1 - γ 2(f)

the situation with noise at the output. As H1(f) gives the optimum estimate of H(f) for the situation of noise at output, H2(f) will similarly in the situation of noise at input give the optimum estimate of H(f) relating as much of the input signal GAA(f) to the output signal GBB(f) minimizing the amount of noise GMM(f) (see Appendix A).

4.4. Other Inputs If there are other inputs to the system, which are not measured in a(t), but which contribute to the measured output b(t) we have the situation as depicted in Fig.17. This is encountered in practice when a controlled and measurable input a(t) is applied to a system which is also excited by some operational and unmeasurable inputs m(t) or if the system is excited by more operational inputs and only one of these is measured. Let us first assume that the inputs a(t) and m(t) are uncorrelated i.e. GAM(f) = 0.

Fig. 17. Ideal system with another input signal ing through the system 33

We thus have

Bi(f) = H(f) (Âi(f) + Mi(f))

giving

GBB(f) = H(f)2 (GAA(f)+ GMM(f))

and

GAB(f) = H(f) GAA(f) since GAM(f) = 0

From this we get and

GAB(f) = H(f) GAA(f) H(f)2 (GAA(f) + GMM(f)) G (f) H2(f) = BB = GBA(f) H*(f) GAA(f) H1(f) =

(4.20) (4.21)

= H(f) (1 + GMM(f)) GAA(f) If the inputs a(t) and m(t) are uncorrelated this situation is equivalent to the situation with uncorrelated noise at the output. H1(f) has no bias error while H2(f) overestimates the amplitude with the "noise to GMM(f) signal ratio" term (1 + G (f) ) = AA As expected the Autospectrum method gives Ha(f)2 =

GMM(f) GBB(f) ) similar to H2(f). = H(f)2 (1 + GAA(f) GAA(f)

(4.22)

As in the previous situations the phase of H1(f) and H2(f) are the same and both correct. If the input signals a(t) (measured) and m(t) (not measured) are assumed to be correlated the formulae will be changed as follows: Bi(f) = H(f)(Ai(f) + Mi(f)) giving

GBB(f) = | H(f) |2 ( GAA(f) + GMM(f)+ GAM(f) + GMA(f) )

and

GAB(f) = H(f) ( GAA(f) + GAM(f)) GAM(f) ) GAA(f)

(4.23)

GMM(f) + GAM(f) ) GAA(f) + GMA(f)

(4.24)

which can be written as H1(f) = H(f) (1 + and

H2(f) = H(f) (1 +

and

 Ha(f) 2 =  H(f) 2 (1 +

34

GMM(f) + GAM(f) + GMA(f) ) GAA(f)

(4.25)

All the estimates have a bias error in amplitude and note that H1(f) and H2(f)) can be biased in phase as well since GAM(f) can be complex. For a detailed discussion of analysis of systems with correlated inputs see Ref. [4].

4.5. Noise at both Input and Output This general situation is shown in Fig.18. The noise signals m{t) and n(t) are assumed to be uncorrelated with each other and with u(t) and v(t). From the fundamental equations: GAA(f) = GUU(f) + GMM(f) 2

GBB(f) =  H(f)  GUU(f)+GNN(f) and

and

GAB(f) = GUV(f) = H(f) GUU(f) 1 1 + GMM(f) GUU(f)

H1(f)

= H(f)

H2(f)

= H(f) (1 + GNN(f)) GVV(f)

Ha(f)

2

1 + GNN(f) GVV(f) = H(f) ( 1 + GMM(f) GUU(f) 2

(4.26)

(4.27)

)

(4.28)

It can be seen that  H1(f)  underestimates  H(f)  due to the noise at input and that  H2(f)  overestimates  H(f)  due to the noise at output.

Fig. 18. Ideal system with extraneous noise in both measured input signal and measured output signal 35

Thus by combining  H1(f)  and  H2(f)  we have a lower and an upper bound of the true  H(f)  since  H1(f)  ≤  H(f)  ≤  H2(f) 

(4.29)

The phase of H1(f) and H2(f) are the same and correct, given by the Cross Spectrum.  Ha(f)  will either underestimate, overestimate or be equal to  H(f)  depending upon the relation between the noise to signal ratios GMM(f) GNN(f) and GVV(f) GUU(f)

The Coherence is

γ2(f) = =

GVV(f) GUU(f) ( GUU(f) + GMM(f) ( GVV(f) + GNN(f)) 1 1 + GMM(f) 1 + GNN(f) ) ) ( ( GVV(f) GUU(f)

(4.30)

Presence of uncorrelated noise signals is detected by the Coherence Function, but it cannot distinguish between input and output noise. As a consequence of this the previously defined Coherent Power ≡ γ2(f) GBB(f) only gives the true coherent power GVV(f), if GMM(f) = 0. Likewise the Non-Coherent Power = (1 - γ2(f)) GBB(f) will only give the correct result GNN(f), if GMM(f) = 0.

4.6. Leakage in the Analysis As discussed in Section 3, leakage can cause deformation of the estimated spectra (Fig.10) and therefore lead to bias errors in the Frequency Response Function estimates. If the resolution in the analysis is too coarse compared to the bandwidth of the system resonances and the input signal is random, the Coherence Function will detect this leakage by having a value less than one at the resonance frequencies (see Figs.11 and 12). The Coherence can therefore give a warning of potential bias errors in the Frequency Response Function estimates. Fig.19 gives an example where a Frequency Re36

Fig. 19. Magnitude of Frequency Response Function estimate H1(f) and Coherence Function from a baseband measurement on a mechanical system. Random noise excitation

sponse Function is estimated in the frequency range 0 - 3,2 kHz (baseband analysis). The system is a mechanical structure and the Frequency Response Function is the point accelerance (acceleration over force) estimated by H1(f). The structure is excited with a random force signal from a shaker, driven by a generator and a power amplifier. The force signal is measured with a force transducer while the response is measured with an accelerometer. A number of resonances (modes of vibration) in the structure are detected and their approximate frequencies can be found from this measurement. However, the low Coherence around the resonance peaks are caused by the lack of resolution in the analysis and the leakage effect will introduce bias errors in the estimates of the peak amplitudes. For the second resonance, for example, the maximum peak level at 956 Hz is 31,1 dB and the Coherence is 0,4 (see the cursor readout in the upper and lower right corners of the plot). Increasing the resolution by performing a so-called zoom analysis should remedy this leakage problem. The result of a zoom analysis with a frequency span of

37

Fig. 20. Magnitude of Frequency Response Function estimate H1(f) and Coherence Function from a zoom measurement (Frequency span of 200 Hz) on the same system as in Fig.19. Random noise excitation.

200 Hz around the second resonance is shown in Fig.20. The Coherence in this analysis is very close to one around the resonance peak indicating elimination of leakage effects and also proving that the low Coherence in the baseband measurement was due to leakage and not extraneous noise or non-linearities. The exact resonance frequency is 953,5 Hz and the true level is 42,2 dB (see cursor read out in Fig.20). The peak level in the baseband analysis was 31,1 dB, indicating severe leakage in the baseband analysis. Let us instead of H1(f) use H2(f) as the baseband estimate of the driving point accelerance. The result is shown in Fig.21, which is based on the same measurement as the one used for calculating H1(f) in Fig.19. The peak amplitude in H2(f) for the second resonance is 39,0 dB (at

38

Fig. 21. Magnitude of Frequency Response Function estimate H2(f) and input Autospectrum for the same measurement as in Fig. 19.

956 Hz, see cursor read-out) i.e. there is only a bias error of 3,2 dB in H2(f)  compared to an error of 11,1 dB in H1(f). H2(f) is thus much less sensitive to the leakage effect at the resonances and a better estimate of peak amplitudes can be found using H2(f) instead of H1(f) if the resolution is too coarse. If the resolution is high enough and the Coherence is one, then H1(f) = H2(f). A mathematical formulation and proof of this is found in Ref. [6]. Only an intuitive explanation shall be given here. Suppose that the input spectrum is flat and therefore has no bias error (this is not true for the example in Fig.19 and 21). The leakage (due to the convolution of the spectrum with the Fourier Transform of the time weighting function, as shown in Fig.10) will cause the levels in GBB(f) and GAB(f) to be underestimated around the resonances. This is illustrated in Fig.22. The bias error in GAB(f) will give the bias error in H1(f) as

39

Fig. 22. Bias errors in GBB(f) and | GAB(f) due to leakage. GAA(f) is as sumed flat

H1(f) =

GAB(f) GAA(f)

and GAA(f) is assumed flat.

The bias errors in GBB(f) and GBA(f) ( = GAB(f) ) tend to be of the same order of magnitude and therefore tend to cancel in the H2(f) GBB(f) estimate, H2(f) = GBA(f) It is shown in Ref. [6] that also in the situation where the input spectrum drops at the resonance frequency of the system, H2(f) is less susceptible to bias errors than H1(f). Such situations very often occur in practice when a shaker via a force transducer is used for excitation of a mechanical system, as mentioned earlier. For example, this is the case for the first and second resonances in Figs.19 and 21. The force spectrum is shown in the lower graph of Fig.21. Furthermore, it could also give rise to problems with background noise for the input force signal at the resonance frequencies (see Section 4.3) and could therefore be another reason for using H2(f) as the estimate of the Frequency Response Function at the resonances. In this discussion it has been assumed that there is no extraneous noise at the output. If there are other uncorrelated inputs to the system which excite the resonances in question as well, H2(f) might overestimate the Frequency Response Function. As argued in section 3 the Coherence Function will only detect leakage when the signals are random in nature. If the signal repeats itself for each data record, the Coherence will be one (assuming no extraneous

40

noise in the measurement) and H1(f) = H2(f), according to equation (4.8). This can be the situation when an impact hammer is used for excitation of a mechanical structure (see Section 5).

4.7. Summary From a Dual Channel FFT measurement two different estimates of the complex Frequency Response Function can be calculated: H1(f)=

GBA(f) G (f) and H2(f) = BB GBA(f) GAA(f)

To summarize, 1. When there is extraneous noise at output, or several independent inputs to the system H1(f) should be used, 2. When there is extraneous noise at input, H2(f) should be used. 3. If there are problems with leakage at resonance peaks (resolution bias) H2(f) gives a better estimate than H1(f). In a practical measurement, different situations occur at different frequencies. In order to get an optimal Frequency Response Function estimate from one measurement, H1(f) should be used at some frequencies and H2(f) at other frequencies. Usually, H1(f) should be used for the notches, where output noise tend to dominate, and H2(f) at the peaks, where input noise or leakage tend to cause the problems. In situations with both input- and output-noise and when leakage does not cause any bias errors H1(f) and H2(f) can be used to give a lower and upper bound respectively, for the true | H(f) |.

5. Excitation Techniques In most applications of system analysis or testing, where input-output relationships have to be measured, it is necessary to excite the system with a well controlled and measurable input. A number of different types of excitation signals are available for the analysis each having its own advantages and disadvantages. For struc41

Fig. 23. Different excitation techniques for structural analysis or testing

tural testing for instance, it is possible to use either an impact hammer or a shaker as shown in Fig.23. If the shaker is selected there are a number of generator signals which can be used such as random, pseudo random, sine etc. The choice of the signal depends among other things upon the test application, non-linear behaviour of the system and time available for the analysis. In this section the most common types of excitation techniques will be described and discussed briefly in of advantages and disadvantages. The types described here will be 1) Random, 2) Pseudo Random, 3) Periodic Impulse, 4) Periodic Random, 5) Sine, 6) Impact.

5.1. Random Noise Excitation A random signal such as shown in Fig.24 is a continuous type of signal which never repeats itself and whose amplitude can only be predicted in of statistical parameters. It is usually described in of its power spectral density (Autospectrum), its Autocorrelation Function (containing the same information as the Autospectrum) and the amplitude probability density. The Autocorrelation Function is related to the 42

Fig. 24. Random noise signal

Autospectrum via a Fourier Transform and will be discussed in Section 7. The amplitude probability density is defined as the probability of having amplitude values within a certain amplitude interval ∆x, divided by the size of that interval ∆x. This function gives a description of how on the average, the instantaneous amplitude of the signal is distributed as a function of amplitude level. It has however no information about the time history or frequency content. Often the random signals found in practice have a gaussian amplitude probability density distribution (see Ref. [4]). The random signal should preferably have a constant spectral density in the frequency range of interest i.e. the Autospectrum should be flat in this frequency range and the signal is called a band limited white noise signal. The main characterisitic of the random excitation signal is that the spectral estimates Âi(f) for each recorded data block will have random amplitude and random phase. At each frequency the system can thus be considered as being excited by different amplitude and phase in each data block analysis. Considering the effects of non-linearities as noise at the output (see Section 4.2), the H1(f) estimate will give the best linear fit to the system or the optimum Frequency Response Function minimizing the effects of the non-linearity (also called the least squares estimate, see Ref. [4]). This is a very important advantage of random excitation. Another advantage of using random noise excitation is that it can be fairly easily shaped to fit the frequency range of interest by filtering and modulating the original broad band white noise signal. Thus the system is not excited by frequencies outside the analysis bandwidth, giving a better dynamic range in the analysis. Fig.25 gives an illustration for a baseband and a zoom measurement. The random signal being continuous does not fit the block length in the analysis. A smooth weighting function (such as the Hanning Weighting) therefore has to be applied which causes leakage in the spectral esti-

43

Fig. 25. Filtering of the random excitation signal giving better dynamic range in the analysis

mates. This is exemplified in Fig.19 and is one of the disadvantages of the random noise excitation. The H2(f) estimate is less susceptible to effects of leakage (section 4.6, Fig.21) and could thus be used, or a zoom analysis would have to be performed (Fig.20). Zooming will, however, increase the analysis time and may therefore not be preferable. It should also be kept in mind that for structural testing (Fig.23) although the shaker is being driven by a signal, which has a flat spectral density in the frequency range of the analysis (Fig.25), the impedance mismatch between the structure and the shaker will cause the input force signal to drop at the resonance frequencies of the structure. This is however usually not a serious problem as H2(f), which is insensitive to extraneous noise at the input (section 4.3), could be used instead (assuming no output noise). 5.2. Pseudo-Random Excitation The pseudo-random signal is specially designed for the DFT analysis which works on blocks of data. It is made up of a segment of a "random" signal of length T which is repeated after every period of time T, see Fig.26. It is periodic and therefore has energy only at discrete frequencies f = k1/T, where T is the period length and k is an integer. The period length T is matched to the record length of the Analyzer, so the frequency components of the pseudo-random signals coincide with the computed frequency lines in the analyzer. One record length of the pseudo-random signal therefore contains all the information in the signal. Rectangular weighting should be used and there will be no leakage in the spectral estimates. This is probably the main advantage of using pseudo-random excitation. The signal is designed in such a way that each frequency component has the same amplitude in the frequency

44

Fig. 26. Example of a pseudo-random signal

range of interest. The phase angle between the different components, however, will be random. The pseudo-random signal can therefore also be considered as a number of sinewaves, having the same amplitude in the analysis frequency range but a random phase, and where the time record in the analyzer contains an integer number of periods for each sinewave. The spectrum can be shaped in the same way as the random signal. This is shown in Fig.27. In the Brüel & Kjær Dual Channel Signal Analyzer Type 2032 or 2034 which is a 801 line analyzer, the pseudo-random signal contains 801 components (sinewaves) in the frequency span which is selected. As the signals involved are deterministic and periodic only a few averages are needed in the analysis, assuming that there is no extraneous noise at the input or output. This can be an advantage specially for low frequency work and zoom analysis, compared to random noise excitation where some averaging always is needed in practice. The main advantages in using pseudo-random excitation are: no leakage in the

Fig. 27. Spectrum of pseudo-random signal for baseband analysis and for zoom analysis 45

analysis, the spectrum can be shaped to only excite frequencies in the range of interest and only a few averages are needed. For a linear system the Frequency Response Function will be calculated at the discrete frequencies fk = k ∆ f = k1/T without any bias error, as there is no leakage in the analysis. This is illustrated in Fig.28 where the same structure as used in Figs.19, 20 and 21 is analysed using pseudorandom noise. The Frequency Response Function estimate at 952 Hz is 38,8 dB (compared to 30,6 dB in  H1(f)  with random noise and the same resolution of ∆f = 4 Hz). Notice that the coherence is one as there is no leakage and no extraneous noise in the measurement. Zooming in on the resonance (Fig.29) using random noise and a resolution of ∆f = 0,25 Hz (no effects from leakage) shows that the true optimum linear estimate of the Frequency Response Function is 38,7 dB at 952,00 Hz as computed in the baseband pseudo-random test (within 0.1 dB). This is advantageous when the data are used for analysis and modelling of linear systems. The effect of sampling the Frequency Response Function

Fig. 28. Magnitude of Frequency Response Function estimate H1(f) and Coherence Function from a baseband measurement using pseudorandom excitation. Same structure as used in Fig. 19, 20 and 21 46

Fig. 29 Zoom analysis around the second resonance (same measurement as in Fig.20). Notice that the magnitude of the Frequency Response Function estimate at 952 Hz is the same as the estimate at 952 Hz from the baseband measurement using pseudorandom excitation (Fig.28) within 0,1 dB

at discrete frequencies is often called the picket fence effect. (see for inst. Ref. [1] and Ref. [2]). However it is very important to keep in mind that this only works for linear systems. As the excitation signal is periodic, non-linearities in the system will also be excited periodically and their effects will therefore not be averaged out as was the case for random excitation. An optimum linear estimate of a non-linear system is therefore not found. Considerable differences in the Frequency Response Function estimates can thus be found between a random test and a pseudo-random test on a non-linear system.

47

5.3. Periodic Impulse Excitation A short impulse repeated every T, where T is the record length in the Analyzer, is called a periodic impulse signal, and is shown in Fig.30. It has a line spectrum where the components coincide with the lines in the Analyzer and can therefore be considered as a special case of a pseudo-random signal. The impulse is so short (one sample time ∆t) that the spectrum is nearly flat in the baseband frequency range, see Fig.30. Rectangular Weighting can be used and leakage in the analysis is avoided as for the pseudo-random signal. Another advantage of the periodic impulse technique is that gating of direct signal and the reflected signals can be performed in some applications, for example, in acoustics. The gating is done by use of a transient time weighting function whose delay and length in the data record is defined by the . For broadband analysis the delay time between direct sound and reflected sound however, should be relatively short as the record length T becomes short. In the Brüel & Kjær Analyzer Type 2032 or 2034 the record length T is ≈ 31 ms for a 25,6 kHz frequency span, so that the delay between direct sound and reflections should be less than ≈ 20 ms. The disadvantages however, are as follows: No linear approximation of a non-linear system is obtained for the same reasons as found in the discussion of pseudo-random noise. The signal has a very high crest factor (ratio of peak to RMS level). The high peak level might give problems in exciting non-linearities in some kind of systems. Also weighting functions (transient in channel A and transient or exponential in channel B) might have to be applied in order to increase the signal to

Fig. 30. Periodic Impulse excitation signal and its Autospectrum 48

noise ratio in the analysis. The periodic impulse signal should only be used in baseband analysis as the main energy in the signal always is concentrated at the low frequencies i.e the main lobe of the sin x x shaped envelope curve shown as a dotted line in Fig.30.

5.4. Periodic Random Excitation With the periodic random signal the advantages of the random test and the pseudo-random test are combined. As illustrated in Fig.31 it consists of a pseudo-random sequence (A) of length T which is repeated a couple of times followed by another sequence (B) of length T, which is independent of A repeated the same number of times. A third sequence (C) independent of the previous ones is now repeated and so on. As indicated in Fig.31 the first couple of periods of each sequence are used for the transient response of the system after the change of sequence i.e. change of phase and amplitude between all the sinewave components. The last sequence is then used for the analysis where the system is in a quasi stationary condition. Rectangular weighting should be used and there will be no leakage. In addition it will give the best linear approximation of the system, as the different pseudo-random sequences are independent of each other (changed randomly in phase and amplitude). Compared to pseudo-random testing the periodic-random testing will usually be slower since not all the blocks in the signals are analysed.

Fig. 31. Periodic random excitation signal 49

5.5. Sine Testing Excitation with a sinewave has been the traditional input signal in system analysis for many years and is still widely used. For Frequency Response Function testing, the sinewave is either stepped or swept through the frequency range of interest, while the output signal is being recorded and the input signal is controlled via a compressor loop to maintain a constant level. Some of the advantages of this type of testing is that the input signal level is well controlled, that the signal to noise ratio is good and that the crest factor of the signal is low. Also, a sine wave is the best input signal when non linearities have to be studied. For mechanical systems having non-linear spring elements swept sine testing can give very useful information about the behaviour of these "stiffness" elements. Examples of this kind are found in Ref. [9]. Study of harmonic distortion will also require a sinusoidal excitation signal. The main disadvantage of the traditional swept sine test is that it is slow compared to the other methods where all the frequencies are excited and analysed simultaneously. A special kind of swept sine testing is the so-called TDS. This technique is outside the scope of this article, but a discussion of TDS theory and its application is found in Technical Review Nos.1 & 2, 1983.

5.6. Impact Excitation A very popular and convenient excitation technique for mechanical structures is impact excitation using an impact hammer. Fig.32 shows the Brüel & Kjær Impact Hammer Type 8202, with a force transducer on which an impact tip is mounted. When the structure is excited by the hammer, energy is transferred to the structure in a very short period of time, giving a typical Input force signal as shown in Fig.32. The shape of this force signal depends upon the type of the hammer tip, mass of the hammer and the dynamic characteristics of the structure under investigation. As the frequency bandwidth of the force spectrum is determined by the length of the signal, these characteristics will determine the cut off frequency of the excitation signal. The more stiff the hammer tip and the structure is, the shorter will the signal be, and the wider will the frequency span be. Extra mass on the hammer, shown as dotted lines in Fig.32, will widen the force signal and therefore lower the cut off frequency. This is illustrated in Fig.33, where curves for steel-, plastic- and rubber-tip are shown with and without the extra mass on the hammer. The structure used was a stiff plate. 50

Fig. 32. Impact Hammer (Brüel & Kjær Type 8202) and a typical input force pulse

It should be kept in mind that the force measured by the force transducer is not the true force which is input to the structure. Denoting the effective mass of the tip as m and the effective mass of the hammer + m where F (minus tip) as M the true force F will be given by F = F M m is the measured force. The effective masses take into the material characteristics of the tip and the hammer. Therefore the impact hammer should always be calibrated using a proper calibration mass before the measurements are performed. This is described in Ref. [8]. The advantages of using impact testing are: a) It is very fast, only a few averages are needed; b) No elaborate fixturing as for shaker excitation is required; c) Easy to use in the field. However, the force signal has a high crest factor which can make this technique non suitable for nonlinear systems. The limited control of excitation bandwidth is also a disadvantage.

Fig. 33. Typical Autospectra for Impact Hammer force pulses when steel, plastic-, and rubbertip are used, with and without extra mass 51

Fig. 34. Analysis of a mechanical structure using Impact excitation a) Time ch.A and Time ch.B b) Autospectrum ch.A and Autospectrum ch.B c) Frequency Response Function H1(f) (magnitude) and Coherence Function

A practical application of the impact hammer technique on a mechanical structure is shown in Fig.34. A transient weighting function is applied to the force signal in ch.A, which will increase the signal to noise ratio and exclude other force signals which are not fed into the structure. Exponential weighting is applied to the response signal (acceleration) in channel B in order to avoid leakage effects due to truncation of the signal at the end of the record and in order to increase the signal to noise ratio in the analysis. This exponential weighting will introduce leakage but as it will have the same effect as additional damping in the system, correction can be made subsequently. If the time constant for the exponential weighting is τw the corresponding decay rate is σw = 1/τw. The measured decay rate σ thus has to be corrected in order to determine the true decay rate of the system σ. σ = σ + σw ⇒ σ = σ - σw. If the damping ratio (ratio of critical damping) 52

σ 1 , is used, the relation is ζ = ζ - ζw = ζ 2πf0 2πf0τw where f0 is the resonance frequency in question. This works very well as long as the damping of the system is not too small compared to the artificial damping caused by the weighting function. ζ =

In Fig.34.(b) the input force spectrum and the response spectrum are shown. A steel tip was used on the hammer and it can be seen that the force energy spectral density has dropped only 8 dB at 3,1 kHz compared to the levels at low frequencies. The magnitude of the estimated Frequency Response Function  H1(f)  and the Coherence γ2(f) are shown in Fig.34.c). Notice that the Coherence is one at all the resonance frequencies. The exponential weighting function used on the response signal will introduce leakage in the analysis, but since the signals repeat themselves the leakage Bileakage(f) will be the same for each analysis (see Section 3, Fig.11). The leakage is therefore not detected by the Coherence Function. The most important advantages and disadvantages of the different excitation techniques discussed in this section are summarized in Table 1.

Leakage in analysis

Best linear fit of non-linear system

Crest factor

Random

Yes

Yes

Pseudorandom

No

No

Periodic Impulse

Depends on weighting functions

Periodic Random Sine Impact

Signal to noise ratio

Control of excitation bandwidth

Speed

Medium

Fair

Good

Fast

Medium

Good

Good

Very fast

No

High

Poor

Limited (no zoom)

Very Fast

No

Yes

Medium

Good

Good

Slower than random

Can be avoided

No

Low

Good

Good

Slow

Depends on weighting functions

No

High

Poor

Limited (no zoom)

Very Fast

Table 1.

53

References [1] THRANE, N.

"The Discrete Fourier Transform and FFT Analyzers", B & K Technical Review, No.1 - 1979.

[2] BRIGHAM, E. O.

"The Fast Fourier Transform", Prentice Hall Inc., New Jersey, 1974.

[3] RANDALL, R. B.

"Application of B & K Equipment to Frequency Analysis", B & K, 1977.

[4] BENDAT, J. S. & PIERSOL, A. G.

"Engineering applications of Correlation and Spectral Analysis", John Wiley & Sons, N.Y. USA 1980.

[5] MITCHELL, L.D.

"Improved Methods for the FFT calculation of the Frequency Response Function", Journal of Mechanical Design, April 1982, Vol. 104.

[6] CAWLEY, P.

"The Reduction Bias Error in Transfer Function Estimates using FFT Based Analyzers", MSA - Session, ASME Conf., Dearborn Ml. Sept. 1983.

[7] BROWN, D., CARBON, G. RAMSEY, K.

&

"Survey of Excitation Techniques Applicable to the Testing of Automotive Structures", Society of Automotive Engineers, Detroit, Feb.-March, 1977.

[8] HALVORSEN, W.Q. & BROWN, D.L.

"Impulse Technique for Structural Frequency Response Testing." Sound and Vibration, Nov. 1977.

[9] BROCH, J. T., et.al.

"Mechanical Vibration and Shock Measurements", B & K 1980.

[10] MARKEL, J. D. & GRAY, A. H.,Jr.

"Linear Prediction of Speech", Springer - Verlag, New York, 1976.

[11] SCHAFER, R. W. & RABINER, L. R.

"Digital Representations of Speech Signals", Procc. IEEE Vol. 63, pp. 662 - 677, April 1975.

54

APPENDIX A Optimum Estimate of H(f) for the situation where there is noise in the input measurement (Section 4.3)

In ref.[4] it is shown that H1(f) gives the optimum estimate of H(f) in the situation with output noise (Section 4.2, Fig.15). The optimum estimate is defined as that estimate which minimizes the extraneous noise GNN(f). From this it can easily be proved that H2(f) is the optimum estimate of the Frequency Response Function in the situation where the measured input signal is contaminated with extraneous uncorrelated noise (section 4.3, Fig.16). Viewing the system in the opposite i.e. backward direction and defining an inverse Frequency Response Function Hinv.(f) given by 1 we have the situation shown in Fig.A.1. Hinv. (f) = H(f)

Fig. A.1. Situation with noise in the input measurement (Section 4.3, Fig.16). The system is viewed in the opposite i.e. backward direction and the inverse Frequency Response Function Hinv.(f) is defined by Hinv.(f) = 1 /H(f)

The relations are:

and

U(f) = Hinv.(f) B(f)

(A.1)

A(f) = Hinv.(f) · B(f) + M(f)

(A.2) 55

a(t) is considered as the measured output signal and b(t) is considered as the measured input signal. This situation corresponds to the situation where the output signal is contaminated with extraneous uncorrelated noise (Section 4.2, Fig.15). The estimate of Hinv.(f) which minimizes the Autospectrum GMM(f) of the noise will thus be given by H1(f) for the inverse situation in Fig.A.1 i.e. Hinv.(f) =

GBA(f) GBB(f)

(A.3)

The optimum estimate of the Frequency Response Function H(f) in the situation where the measured input signal is contaminated with extraneous noise (Fig.16, section 4.3) is therefore given by:

H(f) =

56

1 G (f) ≡ H2(f) = BB GAB(f) Hinv.(f)

(A.4)

PREVIOUSLY ISSUED NUMBERS OF BRÜEL & KJÆR TECHNICAL REVIEW

(Continued from cover page 2) 2-1977 Automated Measurements of Reverberation Time using the Digital Frequency Analyzer Type 2131. Measurement of Elastic Modulus and Loss Factor of PVC at High Frequencies. 1-1977 Digital Filters in Acoustic Analysis Systems. An Objective Comparison of Analog and Digital Methods of Real Time Frequency Analysis. 4-1976 An Easy and Accurate Method of Sound Power Measurements. Measurement of Sound Absorption of rooms using a Reference Sound Source. 3-1976 Registration of Voice Quality. Acoustic Response Measurements and Standards for MotionPicture Theatres. 2-1976 Free-Field Response of Sound Level Meters. High Frequency Testing of Gramophone Cartridges using an Accelerometer. 1-1976 Do We Measure Damaging Noise Correctly? 4-1975 On the Measurement of Frequency Response Functions. 3-1975 On the Averaging Time of RMS Measurements (continuation). SPECIAL TECHNICAL LITERATURE As shown on the back cover page, Brüel & Kjær publish a variety of technical literature which can be obtained from your local B & K representative. The following literature is presently available: Mechanical Vibration and Shock Measurements (English), 2nd edition Acoustic Noise Measurements (English), 3rd edition Acoustic Noise Measurements (Russian), 2nd edition Architectural Acoustics (English) Strain Measurements (English, German, Russian) Frequency Analysis (English) Electroacoustic Measurements (English, German, French, Spanish) Catalogs (several languages) Product Data Sheets (English, German, French, Russian) Furthermore, back copies of the Technical Review can be supplied as shown in the list above. Older issues may be obtained provided they are still in stock.

Printed in Denmark by Nærum Offset