Impulse response measurements using MLS

by

Jens Hee http://jenshee.dk

August 2003

Change log
20. may 2003 1. Document started. 11. july 2003 1. First version released. 30. july 2003 1. The section about primitive polynomials is replaced by a section about nite elds. 2. A program example is added. 1. august 2003 1. Extension of reordering example in section Deconvolution step by step. 2. Bibliography added. 6. august 2003 1. Formula for deconvolution is corrected for DC oset

Contents
1 Introduction 1.1 What is MLS? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Measurement techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 The 2.1 2.2 2.3 Maximum Length Sequence Finite elds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . MLS by recursion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Autocorrelation function and crest factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 3 3 4 4 5 5 6 6 7 7 9

3 Computing the impulse response 3.1 Deconvolution . . . . . . . . . . . 3.2 The M -sequence matrix . . . . . 3.3 Decomposition of the M -sequence 3.4 Hadamard matrix equivalence . . 3.5 Deconvolution step by step . . . . 3.6 Program example . . . . . . . . .

ii

Chapter 1 Introduction
1.1 What is MLS?

MLS (Maximum length sequence) for signal processing purposes is the name of a pseudo random signal consisting of the values 1 and -1. It is periodic with the period P = 2N 1 and has a at frequency response (see Chapter 2). MLS is also the name of a computational cost eective technique for measuring the impulse response of a linear system using the MLS signal as an input to the system (see Chapter 3).

1.2

Measurement techniques

The impulse response and the frequency response of a system are the Fourier transforms of each other. It is therefore only necessary to measure one of them, then the other one is known. In the following some of the more common techniques are briey reviewed.

Impulse excitation
A strait forward method is to use an impulse as the excitation signal. The output from the system is then by denition the impulse response. The draw back of this method is that the impulse must be short enough to have a at frequency response over the frequency range of interest and at the same time produce enough energy to give a reasonable signal to noise ratio without overloading the system. The signal to noise ratio can be improved by averaging, but the measurement time will increase correspondingly.

Linear sweep excitation

If a linear sweep is used for excitation, the frequency response of the system under consideration can be obtained by mixing the output of the system with the applied sweep signal followed by a low-pass ltering. In practice the method is more complicated than indicated, but the details are not given here. Since a sweep is used for excitation a high signal to noise ratio is obtained. Lowering the sweep rate can increase the S/N ratio. The method can also be used for measuring the distortion of non-linear systems. 1

Stepped sine excitation

Stepped sine analysis is carried out similar to the linear sweep method except that the excitation signal is stepped through a sequence of frequencies. Normally this method is used to obtain a logarithmic frequency axis and is not suitable for measuring the impulse response.

MLS excitation
The Impulse response of a system can be measured by applying an MLS signal and process the output by a modied Hadamard transform. The details are given below. The method has several advantages. The processing is simple (although theoretically complex), which can be important for small computer systems. The method gives a good S/N ratio due to the low crest factor of the MLS signal.

Two-channel FFT
The previous methods all use a particular excitation signal. It is also assumed that the measurement channel has a linear amplitude and phase characteristic. If both input and output of the system are measured and processed using FFT, the frequency response of the system can be obtained using almost any input signal and only amplitude and phase match between channels is required. The method is computational more complex than the MLS method, but gives a high degree of freedom for selecting the excitation signal.

Chapter 2 The Maximum Length Sequence

As mentioned above an MLS signal is a pseudo random periodic sequence. It can be generated by a one-bit linear feedback shift-register whose characteristic polynomial is primitive (see below), but for signal processing, the binary states of 0 and 1 are mapped into 1. The theory about primitive polynomials is part of the rather extensive theory about nite elds also called Gallois elds. The next section gives some facts from this theory.

2.1

Finite elds

A nite eld is a nite set where addition, subtraction, multiplication and division are dened. Finite elds exist only when the number of elements are pN where p is a prime and N any positive integer. Only elds with p = 2 are of interest for MLS signals. It can be shown that the elements in such a eld can be represented by the 2N polynomials of degree N 1 with coecients 0 and 1 and the multiplication and division is dened modulo an irreducible polynomial of degree N . When adding and subtracting coecients the operation is modulo-2. For N = 3 an irreducible polynomial is x3 + x + 1 because neither 0 nor 1 are roots. The elements in the eld are: 0 000 1 001 x 010 x2 100 x + 1 011 x2 + x 110 2 x + x + 1 111 x2 + 1 101 It is seen that apart from the rst element each element is equal to the previous one multiplied by x. When an irreducible polynomial has this property it is called a primitive polynomial and x is called a primitive element or generating element. For N = 3 both the irreducible polynomials are primitive but for N = 4 only two of the three irreducible polynomials are primitive. It is not simple to nd the primitive polynomials, but several textbooks give lists of primitive polynomials up to a very high degree. Some of them can also be recovered from the program example (see below). 3

2.2

MLS by recursion
p(x) = xn + an1 xn1 + ... + a1 x + a0

If the primitive polynomial is given by:

a maximum length sequence can be calculated from the recursion formula: sk+n = an1 sk+n1 + ... + a0 sk As an example the polynomial x3 + x + 1 is primitive and the recursion sk+3 = sk+1 + sk can be used to generate the sequence: 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, ... Note that the sequence repeats itself after 23 1 = 7 samples. If an irreducible but not primitive polynomial had been used the resulting repetition would have been shorter. It can be shown that the repetition can be no longer. That is why sequences generated by primitive polynomials are called maximum length sequences. By replacing 0 by 1 and 1 by -1 one gets an MLS signal with period 7: 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...

2.3

Autocorrelation function and crest factor

1 for k = 0 1 P for 0 < k < P 1

The autocorrelation function of the MLS signal is given by: (k) =

It is seen that the spectrum of the sequence is at except for the DC component. The crest factor cf = 1 which is 3 dB lower than for a sine wave.

Chapter 3 Computing the impulse response

3.1 Deconvolution
y =hx If y is convolved with a sequence z: yz =hxz =h If: 1. x z is a sequence of delta functions with the periodicity of x 2. The length of h is shorter than the periodicity of x. then h is a periodic version of the impulse response h. If zn = xn then xz is the autocorrelation function of the sequence x and it is seen from the above that when x is an MLS signal the rst clause is fullled except for a DC oset. The periodicity of the MLS signal must be adjusted in order to fulll the second clause. For an MSL signal h is given by: h (n) =
P P 1 ( y(k)x(k n) y(k)) P + 1 k=1 k=1

If the impulse response of a system is h then the input output relation is given by the convolution:

where P is the periodicity of the sequence. In matrix form the deconvolution term for a 7 sample sequence can be written:

h1 h2 h3 h4 h5 h6 h7

x1 x2 x3 x4 x5 x6 x7

x2 x3 x4 x5 x6 x7 x1

x3 x4 x5 x6 x7 x1 x2

x4 x5 x6 x7 x1 x2 x3 5

x5 x6 x7 x1 x2 x3 x4

x6 x7 x1 x2 x3 x4 x5

x7 x1 x2 x3 x4 x5 x6

y1 y2 y3 y4 y5 y6 y7

3.2

The M -sequence matrix

As mentioned above the MLS signal only consists of +1 and 1 and the matrix can be written: h1 h2 h3 h4 h5 h6 h7
=

1 1 1 1 1 1 1

1 1 1 1 1 1 1

1 1 1 1 1 1 1

1 1 1 1 1 1 1

1 1 1 1 1 1 1

1 1 1 1 1 1 1

1 1 1 1 1 1 1

y1 y2 y3 y4 y5 y6 y7

It is clear that no multiplications are required and the number of additions is P (P 1). The number of additions can be reduced to approximately P log2 (P ) by reordering the matrix and using the fast Hadamard transform as shown in the next section.

3.3

Decomposition of the M -sequence matrix

Returning to the binary formulation in Chapter 2 it is easy to see that the matrix M is symmetric and thus not only the rows but also the columns satisfy the recursion formula. Therefore every row (column) of M can be expressed as a linear modulo-2 combinations of the rst N rows (columns) and we can write: M = LS = L S where L is a matrix of order P N , S is the N P matrix formed by the rst N rows of M and the primes denote matrix transposition. M can then be written:
M = LS =

1 1 1 0 0 1 0

1 1 0 0 1 0 1

1 0 0 1 0 1 1

0 0 1 0 1 1 1

0 1 0 1 1 1 0

1 0 1 1 1 0 0

0 1 1 1 0 0 1

1 0 0 1 0 1 1

0 1 0 1 1 1 0

0 0 1 0 1 1 1

1 1 1 1 1 0 1 0 0

0 0 1 0 0 1 0 1 1 0 1 1

Since all rows of M are distinct all rows of L are distinct. Moreover since S consists of the rst N rows of M the rst N rows of L form an identity matrix. If the the rst N columns of S are considered one can write: L = S and since every nonzero binary N -vector appears in both L and S we get: L = S 1 Consequently L consists of columns in M which forms an identity matrix by the rst N rows.

3.4

Hadamard matrix equivalence

H = BB =

It can be shown that the Hadamard matrix can be written in a similar way to the M matrix: 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 0 1 1 0 0 1 1 0 1 1 0 0 1 1 0 0 0 0 0 1 1 1 1 0 1 0 1 1 0 1 0 0 0 1 1 1 1 0 0 0 1 1 0 1 0 0 1
=

0 0 0 0 1 1 1 1

0 0 1 1 0 0 1 1

0 1 0 1 0 1 0 1

0 0 0 0 0 0 1 1 0 1 0 1

1 1 1 1 0 0 1 1 0 1 0 1

Note that H has P + 1 rows and P + 1 columns. If a zero row and a zero column are bordered to the left and top of M , M can be transformed into H by reordering the columns and rows according to L and S. The fast Hadamard transform can then be used to compute the matrix product involving M and the number of additions can be reduced to P log2 (P ).

3.5

Deconvolution step by step

Although theoretically complex the steps necessary to perform the deconvolution are fairly simple and are as follows: 1. Generate the MLS signal. 2. Reorder the output sequence from the system under consideration according to the columns of S and add a zero element in front. If the system under consideration is DC coupled yn should be used as the rst element. 3. Apply the fast Hadamard transform to the reordered output sequence. 4. Omit the rst element of the output from the Hadamard transform, reorder according to the rows of L and divide by P + 1. The reordering of y in the example above gives the following:

1 1 1 7

1 1 0 6

1 0 0 4

0 0 1 1

0 1 0 2

1 0 1 5

0 1 1 3

y1 y2 y3 y4 y5 y6 y7

y4 y5 y7 y3 y6 y2 y1

The reordering of the output of the Hadamard transform gives the impulse response h:

1 0 0 1 0 1 1

0 1 0 1 1 1 0

0 0 1 0 1 1 1

4 2 1 6 3 7 5

h4 h2 h1 h6 h3 h7 h5

h1 h2 h3 h4 h5 h6 h7

3.6

Program example

#include "stdio.h" void GenerateSignal(bool *mls, double *signal, long P) { long i; double *input = new double[P]; for (i = 0; i < P; i++) { input[i] = -2 * mls[i] + 1; } for (i = 0; i < P; i++) { signal[i] = 2.0 * input[(P + i + 0.4 * input[(P + i + 0.2 * input[(P + i - 0.1 * input[(P + i - 0.8 * input[(P + i } delete []input; } void GenerateMls(bool *mls, long P, long N) { const long maxNoTaps = 18; const bool tapsTab[16][18] = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, { 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, // Change 0 to 1 and 1 to -1

// Simulate a system with h = {2, 0.4, 0.2, -0.1, -0.8}, just an example

0) 1) 2) 3) 4)

% % % % %

P] P] P] P] P];

1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

bool long bool long

taps[maxNoTaps]; i, j; *delayLine = new bool[maxNoTaps]; sum;

for (i = 0; i < N; i++) // copy the Nth taps table { taps[i] = tapsTab[maxNoTaps - N][i]; delayLine[i] = 1; } for (i = 0; i < P; i++) // Generate an MLS by summing the taps mod 2 { sum = 0; for (j = 0; j < N; j++) { sum += taps[j] * delayLine[j]; }

sum &= 1;

// mod 2

mls[i] = delayLine[N - 1]; for (j = N - 2; j >= 0; j--) { delayLine[j + 1] = delayLine[j]; } delayLine[0] = *(bool*)&sum; } delete []delayLine; } void FastHadamard(double *x, long P1, long N) { long i, i1, j, k, k1, k2; double temp; k1 = P1; for (k = 0; k < N; k++) { k2 = k1 >> 1; for (j = 0; j < k2; j++) { for (i = j; i < P1; i = i + k1) { i1 = i + k2; temp = x[i] + x[i1]; x[i1] = x[i] - x[i1]; x[i] = temp; } } k1 = k1 >> 1; } } void PermuteSignal(double *sig, double *perm, long *tagS, long P) { long i; double dc = 0; for (i = 0; i < P; i++) dc += sig[i]; perm[0] = -dc; for (i = 0; i < P; i++) perm[tagS[i]] = sig[i]; } void PermuteResponse(double *perm, double *resp, long *tagL, long P) { long i; const double fact = 1 / double(P + 1);

// Just a permutation of the measured signal

for (i = 0; i < P; i++) // Just a permutation of the impulse response { resp[i] = perm[tagL[i]] * fact; } resp[P] = 0; } void GeneratetagL(bool *mls, long *tagL, long P, long N) { long i, j;

10

long *colSum = new long[P]; long *index = new long[N]; for (i = 0; i < P; i++) // Run through all the columns in the autocorr matrix { colSum[i] = 0; for (j = 0; j < N; j++) // Find colSum as the value of the first N elements regarded as a binary number { colSum[i] += mls[(P + i - j) % P] << (N - 1 - j); } for ( j = 0; j < N; j++) // Figure out if colSum is a 2^j number and store the column as the jth index { if (colSum[i] == (1 << j)) index[j] = i; } } for (i = 0; i < P; i++) // For each row in the L matrix { tagL[i] = 0; for ( j = 0; j < N; j++) // Find the tagL as the value of the rows in the L matrix regarded as a binary number { tagL[i] += mls[(P + index[j] - i) % P] * (1 << j); } } delete []colSum; delete []index; } void GeneratetagS(bool *mls, long *tagS, long P, long N) { long i, j; for (i = 0; i < P; i++) // For each column in the S matrix { tagS[i] = 0; for (j = 0; j < N; j++) // Find the tagS as the value of the columns in the S matrix regarded as a binary number { tagS[i] += mls[(P + i - j) % P] * (1 << (N - 1 - j)); } } } void main() { const long N = 18; const long P = (1 << N) - 1; long i; bool *mls = new bool[P]; long *tagL = new long[P]; long *tagS = new long[P]; double *signal = new double[P]; double *perm = new double[P + 1]; double *resp = new double[P + 1]; GenerateMls(mls, P, N); GeneratetagL(mls, tagL, P, N); GeneratetagS(mls, tagS, P, N); GenerateSignal(mls, signal, P); // Generate the Maximum length sequence // Generate tagL for the L matrix // Generate tagS for the S matrix // Do a simulated measurement and get the signal

11

PermuteSignal(signal, perm, tagS, P); FastHadamard(perm, P + 1, N); PermuteResponse(perm, resp, tagL, P); printf("Impulse response:\n"); for (i = 0; i < 10; i++) printf("%10.5f\n", resp[i]); delete delete delete delete delete delete } []mls; []tagL; []tagS; []signal; []perm; []resp;

// Permute the signal according to tagS // Do a Hadamard transform in place // Permute the impulseresponse according to tagL

12

Bibliography
[1] W. T. Chu, Impulse Response and Reverberation-Decay Measurements Made by Using a Periodic Pseudorandom Sequence, Applied Acoustics, vol. 29, pp. 193-205, 1990. [2] M.Cohn and A. Lempel, On Fast M-Sequence Transform, IEEE Trans. Inf.Theory, IT-23, pp 135-137 1977. [3] F. J. MacWilliams and N. J. A. Sloane, Pseudo-random Sequences and Arrays, Proc. Inst Elec. Eng., vol. 64, no 12, pp. 1715-1729, 1976. [4] M. R. Schroeder, Number Theory in Science and Communication, Third Edition 1997, Springer. [5] R. J. McEliece, Finite Fields for Computer Scientists and Engineers, 1987, Kluwer Academic Publishers. [6] R. Lidl & H. Niederreiter, Finite Fields, Enclyclopedia of Mathematics and its Applications 20, 1997, Cambridge University Press.

13