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])
.