# Calculating cross-correlation using Walsh-Hadamard transform

I am trying to implement MLS method of measuring impulse responses. There is an article describing the method: http://www.commsp.ee.ic.ac.uk/~mrt102/projects/mls.... As I understand, to get an impulse response of a system I need to calculate cross-calculation of input signal - MLS and output signal of a system.

I tried to do measurements on simple lowpass RC-filter. The result I get is strange, it is a noisy signal that does not look like lowpass filter signal response. I tried to use convolution and circular correlation but it didn't help.

I suspect the problem is related with calculation methods. I calculate correlation using the convolution theorem that allows to calculate convolution using Fourier transforms of signals. But, the article describes Walsh-Hadamard transform to use in the method.

Also I checked the form of MLS, the generated signal is correct, but I am not sure should it be periodic or not.

My questions are:

Can I get wrong results because I don't use Walsh-Hadamard transform during calculations?

Could you, please, give a sample of code in Matlab or Python that calculate correlation using Walsh-Hadamard transform?

Can there be other sources of problems, for example, distortion related with periodicity of MLS-signal?

P.S. I can show screnshots of wrong impulse responses I got if it is required.

• Why not Google search this code example tinyurl.com/c4z8xxv8 Apr 10 at 6:40
• What is MLS? What is the system? Please include some details in your question instead of expecting readers to check an external link for relevant data. Apr 10 at 15:13

The link provided as a reference is broken. But, regardless of what the OP's source says, the Walsh-Hadamard transform method (multiply the data vector with a $$2^n\times 2^n$$ Hadamard matrix $$H_n$$ in Sylvester form $$H_n = \left[\begin{matrix} H_{n-1} & H_{n-1}\\ H_{n-1} & -H_{n-1}\end{matrix}\right]$$ does not give the cross-correlation of the data vector with a binary maximum-length linear feedback shift register sequence (presumably the MLS referred to by the OP) of period $$2^n-1$$. It is possible to permute the data vector entries to get a close approximation of the result, indeed even the exact result, of the desired cross-correlation, but just applying the Walsh-Hadamard transform to the unpermuted data vector results in garbage, as the OP has discovered.

A description of how the MLS method works can be found in this answer of mine. The necessary permutations to make the Walsh-Hadamard transform method work are probably there somewhere on this forum.

• Thank you very much. Also I'd like to ask a few questions more. Is calculating using Walsh-Hadamard transform mandatory? Can't I just calculate cross-correlation between output and MLS and get periodic impulse response? Apr 10 at 19:19
• "Can't I just calculate cross-correlation between output and MLS and get periodic impulse response?" Sure you can, but be absolutely sure that you are using the raw unadulterated definition of cross-correlation , and none of this nonsense of doing circular convolution via FFTs and so on. Ignore everything you read after the definition of $R_{x,y}(t)$ or $R_{x,y}[k]$ as the integral or sum of a bunch of terms that look like $x(\tau)y(t+\tau) d\tau$ or $x[n]y[n+k]$, especially stuff saying "Hey, lookit! We can calculate these by taking Fourier transforms, multiplying, and inverting!" Apr 10 at 19:32
• Do I understand correctly, I can't use the convolution theorem (that uses exactly FFT) to calculate cross-correlation? Apr 10 at 19:49
• Well, it is possible to use FFTs to do the needed calculations but its use should be restricted to only those who thoroughly understand what they are doing; else the results will be pretty meaningless. I urge you in the strongest possible terms to avoid using FFTs in the calculations. Apr 10 at 20:39

In my opinion, using fast Fourier transform (FFT) to circularly deconvolve a maximum length sequence (MLS) makes sense from a development point of view because FFT is typically readily available and because using FFT makes switching to another test signal easy. If only power-of-two-size FFTs are available, then for a length $$2^m-1$$ MLS, a length $$2^{m+1}$$ FFT can be used with repetition twice of one of the sequences to be convolved and zero-padding of the other.

Using the Walsh–Hadamard transform can be advisable, for power saving reasons if the computation will be done often, or when computing resources are limited.

Here is a Python script for resolving the impulse response by using as the test signal an adjusted MLS that is still two-level, but can be deconvolved to an impulse by circularly convolving it by a time-shifted reverse of the adjusted MLS. The convolution is done using an FFT of length $$2^m-1$$, the same length as the MLS:

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

nbits = 5 # Adjustable parameter: Number of bits in MLS generator linear feedback shift register.

mls_len = 2**nbits - 1 # MLS length
print("MLS length:")
print(mls_len)

mls = signal.max_len_seq(nbits)[0].astype("float64") # Maximum length sequence, each value is either 0 or 1
mls = (1/(2**(nbits/2)+1) + 2*mls - 1) # Fix the DC term problem
mls /= np.sqrt(np.mean(mls**2)) # Normalize
plt.plot(mls, "_")
plt.show()

system_mls_response = mls; # System response to MLS input. For a pass-through system that we use for testing, we have MLS as the response.
# system_mls_response -= np.mean(system_mls_response) # Remove DC

print("Impulse response:")
system_ir = np.fft.ifft(np.fft.fft(system_mls_response)*np.fft.fft(np.roll(np.flip(mls), 1)))/mls_len # FFT-accelerated deconvolution
plt.semilogy(np.abs(np.real(system_ir)), "_")
#plt.ylim([1e-2, 1e0])
plt.show()

print("Impulse amplitude in system impulse response:")
print(np.real(system_ir[0]))
print("Root mean square noise in system impulse response:")
print(np.sqrt(np.mean(np.real(system_ir[1:mls_len])**2)))


Reformatted output:

Figure 1. An adjusted maximum length sequence of length 31. This sequence fixes the DC bin problem (see Fig. 3) and is normalized so that in deconvolution implemented by convolution one only needs to normalize the output by dividing by the MLS length.

Figure 2. Absolute value of the resolved impulse response of a pass-through (identity) system, with a logarithmic vertical axis. The numerical peak amplitude is 1.0 and the root mean square of the noise floor is 3.9983021680982406e-17. These numbers are good enough to conclude that any error is ultimately due to floating point rounding errors propagating to the final result.

Let's also see what we get with just a MLS instead of an adjusted MLS:

Figure 3. The result of circular convolution of an MLS with its time-shifted reverse, with mean removed. The shift in the level illustrates that either an adjusted MLS sequence or other more careful management of the DC term (0 Hz, zero-frequency term) are needed in order to resolve the impulse response. It would work to do deconvolution by frequency-domain division. Also using a longer MLS would reduce the error, but would not allow to resolve the DC bin. To create this plot, the line with comment Fix the DC term problem was commented out, and the line with comment Remove DC was uncommented. The plot vertical range was adjusted by uncommenting plt.ylim([1e-2, 1e0]).