Calculating cross-correlation using Walsh-Hadamard transform

Question

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?

Thank you in advance.

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

Answer 1

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.

Answer 2

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)

print("Adjusted MLS:")
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]).