# Compensating Loudspeaker frequency response in an audio signal

I have been working on a project in which I was required to work on the audio signals recorded from the loudspeaker kept in front of a filter. So, to simply explain it:

$$\boxed{\rm LoudSpeaker} \longrightarrow \boxed{\rm Filter} \longrightarrow\boxed{\rm Microphone}$$

Now the project depends on finding how the filter reacts to the signals from the loudSpeaker. The loudSpeaker outputs a frequency sweep.

When I saw the datasheet of the loudspeaker I found that it had a certain frequency response which was necessary to compensate within the sweep. But now that I have already taken the readings into the microphone I have to subtract the loudSpeaker frequency response from the final spectrum.

A possible way to that is Deconvolution. But I can only explain that in theory.

Can someone help me how to implement these functions in MATLAB? Or a different way to solve this problem???

• Measure without the filter, then measure again with the filter. The filter's linear effect is then the quotient of the two measured frequency responses. Commented Jun 6, 2016 at 9:03
• What is the actual problem though? Your speaker's datasheet is already telling you everything you need to know. Just rectify->integrate your received chirp to get its envelope which is basically your frequency response and then boost / attenuate this curve by the appropriate amount indicated by your speaker's frequency response. In short, all that you have to do is apply the inverse filter which you already know (almost) from the datasheet.
– A_A
Commented Jun 6, 2016 at 15:23

Yes, you can do this with a Wiener Filter which uses the Wiener-Hopf equation to determine the least squared solution to the filter that would compensate for your channel, using the known transmit and receive sequences. The channel is the unknown being solved, and the tx and rx sequences are known. The Wiener Filter is a Minimum Mean Square Error (MMSE) equalizer as it will determine the filter coefficients to minimize the mean square error for a given filter length. (Similarly "LMS Equalizers", which this is not, use an adaptive and iterative algorithm and converge to the same criteria of minimum least squares error).

BOTTOM LINE:

Here is the Matlab function with error checking removed:

function coeff = equalize(tx,rx,depth,ntaps)
%Determines equalizer coefficients using the Wiener-Hopf equations
%TX = Transmitted (Desired) waveform, row vector, length must be > depth+2*ntaps
%RX = Received (Distorted) waveform, row vector, length must be >=depth
%DEPTH = Depth of solution matrix (recommend 10x ntaps but based on duration of stationarity)
%NTAPS = Number of taps for equalizer filter

%force row vectors
tx= tx(:)';
rx= rx(:)';

delay=floor(ntaps/2);
A=convmtx(rx(1:depth).',ntaps);
R=A'*A;
X=[zeros(1,delay) tx(1:depth) zeros(1,ceil(ntaps/2)-1)].';
ro=A'*X;
coeff=(inv(R)*ro);


USE:

Once the coeff for the FIR filter are determined using the function above, then the Matlab filter function can process the receive sequence:

tx_recovered = filter(coeff, 1, rx)

If you want to see the channel response of the filter use:

freqz(coeff)

If you want the solution to be the estimate of the channel instead of the compensation filter that undoes the channel response, simply swap tx and rx:

coeff = equalize(rx,tx,depth,ntaps)

DETAILS FOR THE VERY INTERESTED:

See my slides below giving a high level overview / derivation of the process, This in general form is the Normal Equation (http://mathworld.wolfram.com/NormalEquation.html) used for least squared curve fitting and other applications. I believe I was first introduced to this viewpoint in demonstrating how the Normal Equation is performing deconvolution from the book "Theory and Practice of Modem Design" by John A.C. Bingham.

In practice, I typically do a cross correlation first to determine the channel response time (delay spread) and initial time alignment, and then use an initial equalizer FIR length (# of taps) that exceeds the delay spread (not knowing if leading or trailing echos dominate I will typically start with 2x the delay spread for the FIR length). Once I see the result, the filter size can be reduced if desired based on insignificant magnitudes of the coefficients at the edges of the filter. If the sequences are not exactly aligned, but still within the span of the filter, then the dominant tap will be offset accordingly- so not critical to align beforehand and this gives you insight into what happens if they are grossly misaligned.

Here is an interesting example of the equalizer function I used recently on a sound file from Dalen to equalize the waveforms received by the left and right channels as received by two microphones (treating left as transmit and right as receive and ignoring the actual third party transmitter for the two). The two channels are not recognizable prior to equalization, and completely aligned in amplitude, phase and characteristic after.

Here is a plot of the left and right channels prior to equalization:

Here is the same plot after equalization, right was filtered with the equalizer, and left was filtered with a simple filter just as long as the equalizer with a single unity gain tap in the center and zero elsewhere (to match the delay as the equalizer assumes nominal delay is in the center of the equalizer filter):

This is a zoom in plot of the waveforms after equalization showing how identical the two sequences have become:

How do I know this? I teach courses on DSP and Python related to wireless comm through dsprelated.com and the ieee with new courses running soon! These slides shown are part of my "DSP for Software Radio" course.

• How do I choose or determine (what value to use in) the ntaps value? Could you please add that detail here? Commented Oct 7, 2020 at 9:20
• That would be based on how stationary your signal is. To really understand this the Allan Variance (2 sample variance) gives a lot of insight- I detailed this in another post that you can probably find by searching here under ADEV (I will look for it later otherwise). But to simplify consider a stationary random white process: this we could average for all time and the resulting noise out of the average would keep reducing (by square root of N in standard deviation).... but if the signal wasn’t stationary the standard deviation would start getting worst the longer we average it. Commented Oct 7, 2020 at 13:16
• @skrowten_hermit for cases of 1/f, random walk, drift noise there will be an optimum time that minimizes the average error —- same thing with the related dynamics of a comm channel - hope that helps Commented Oct 7, 2020 at 13:19
• Okay. I see. So, if I have 2 input wav files (mono channel), with N and M samples in them respectively, as mentioned here, aren't they a stationary signal? What value would I pass as ntaps in that case? I just read your answer here (is this the ADEV post you were referring to?), which is concerned with accuracy of phase result in DFT bins. Is that relevant to my problem? Commented Oct 7, 2020 at 14:15
• @mohammadsdtmnd I do not recommend equalizing with IIR: If the channel is not minimum phase (which is often the case), the IIR solution is not stable. Commented Aug 22, 2023 at 10:38