Main Content

Automated Design of Audio Filters for Room Equalization

This example combines Optimization Toolbox™ and Audio Toolbox™ to develop an algorithm that automatically tunes a set of filter parameters.

There are many audio applications where it is desirable to compute parametric equalizer parameters to fit an arbitrary frequency response. For instance, one could fit a filter response to a measured impulse response (IR) to obtain a lower-order implementation of the same filter. Alternatively, one can apply correction to a measured loudspeaker response (anechoic or in-room) to smooth out any imperfections and create a perceptually flat frequency response. The latter is demonstrated here by designing an algorithm that automatically tunes the parameters of N parametric equalizers such that when the resulting EQ is applied to the speaker, the frequency response is perceived as flat in the room.

This example goes over the following steps:

  • Measure the in-room response of a loudspeaker using the Impulse Response Measurer

  • Compute a fractional-octave smoothed response

  • Take into account the microphone calibration data

  • Compute a target response that is perceptually optimal for the given loudspeaker system and room configuration

  • Optimize a set of filter parameters that modifies the response to better fit the target

  • Produce an audio filter or audio files to evaluate the results using headphones or listening room

In-Room Measurement

The first step is to obtain measurements for the system that needs to be improved.

Set up a full duplex sound interface so that it can both play on the loudspeaker, and record with a calibration microphone (such as a Behringer ECM8000 connected to a sound interface capable of supplying phantom power). Place the microphone on a stand so that you can move it into the listening position(s). You can start with the microphone centered in the listener's position, but measuring at several positions will help reduce the chances of overcorrecting for issues like a high frequency dip in the response that could only be present in a region smaller than a listener's head.

Launch the Impulse Response Measurer application.

impulseResponseMeasurer

Verify that the correct audio device is selected. Change the sample rate if desired (this example uses 96 kHz). Set the player and recorder channels to the loudspeaker and microphone, respectively.

Select the Swept Sine method. Set the number of runs to average several measurement together. This example uses five. Set the duration per run to have time for a long enough swept sine and a period of silence that is long enough for the reverberation to completely die off, which can be several seconds in a typical room.

In the advanced settings, set a pause between runs that allows time to move the microphone around the initial "center" position. You must keep silent during the "silence" part of the measurement, but you can move the microphone and make noise during this pause. The start frequency should be set below the range of the loudspeaker (10 Hz might be a good starting point). The stop frequency can be set to half the sampling rate, unless measuring a subwoofer with limited high-frequency range. Set the sweep duration to a few seconds, and make sure there is enough duration left for silence at the end.

You may test levels with the level meter or try a first capture with 1 short run. Set playback level loud enough to hide any background noise in the room and set the microphone level so that it is high but does not overload/clip.

Now you can capture and save export the data to a MAT file. The rest of this example uses a file provided here.

Import the Measurement

Import the last measurement that was exported by the application (by addressing with end). The data used here is in a similar (but compatible) format.

load('measured_ir_data.mat','measurementData');
Fs = measurementData.SampleRate(end);
ir = measurementData.ImpulseResponse(end).Amplitude;
frequency = measurementData.MagnitudeResponse(end).Frequency;
if isfield(measurementData.MagnitudeResponse,'PowerDb')
    magnitudeDB = measurementData.MagnitudeResponse(end).PowerDb;
else % Version R2022b of IRM has renamed this field to Magnitude (dB)
    magnitudeDB = measurementData.MagnitudeResponse(end).MagnitudeDB;
end

Using a helper function provided here, smooth the response by 1/24-octave sections. Since powerDB is a measurement, add extra smoothing (last argument set to true).

fullRange = [10,Fs/2]; % Full audio range (for the plots)
[powerFR,cfFR] = octaveAverage(frequency,db2mag(magnitudeDB),24,fullRange,true);
pdbFR = 20*log10(powerFR);

Plot the measurement and the smoothed response.

semilogx(frequency,magnitudeDB,'g:')
hold on
plot(cfFR,pdbFR,'b',LineWidth=2)
title('Measured Speaker Response')
legend('Raw Data','Octave Smoothed',Location='southwest')
xlabel('Frequency (Hz) \rightarrow')
ylabel('Magnitude (dB) \rightarrow')
yrange = [-50 0];
axis([30 22e3 yrange])
hold off
grid on

Figure contains an axes object. The axes object with title Measured Speaker Response contains 2 objects of type line. These objects represent Raw Data, Octave Smoothed.

Microphone Calibration

If there is calibration data available for the microphone, subtract it from the measurement.

In this case, generic calibration data for the Behringer ECM8000 microphone is used. Calling the helper function without capturing the output produces a convenience plot.

getMicCalibration()

Figure contains an axes object. The axes object with title Microphone Frequency Response contains 2 objects of type line. These objects represent Data Interpolation, Raw Data Points.

Subtract the microphone calibration data from the measured response, using the same frequency values.

micGainDB = getMicCalibration(cfFR);
pdbFRmic = pdbFR - micGainDB;

Determine a "range of interest" (ROI) for the optimization based on the measurement above and the manufacturer specifications for our loudspeaker. The bookshelf speaker measured above has a range of 60 Hz to 22 kHz according to the manufacturer. Set the ROI to a range of 40 Hz to 20 kHz to slightly increase the low end and to take into account the steep decline above 20 kHz.

ROI = [40 20e3]; % specify the region of interest for the optimization

Plot the corrected response and the ROI region.

semilogx(frequency,magnitudeDB,'g:')
hold on
plot(cfFR, pdbFRmic, 'b', LineWidth=2)
plot([ROI(1) ROI(1)], yrange, ':m', LineWidth=1.5);
plot([ROI(2) ROI(2)], yrange, ':m', LineWidth=1.5);
title('Measured Speaker Response')
legend('Raw Data','Octave Smoothed',...
       'ROI Low End','ROI High End',Location='southwest')
xlabel('Frequency (Hz) \rightarrow')
ylabel('Magnitude (dB) \rightarrow')
axis([30 24e3 -50 0])
hold off
grid on

Figure contains an axes object. The axes object with title Measured Speaker Response contains 4 objects of type line. These objects represent Raw Data, Octave Smoothed, ROI Low End, ROI High End.

Compute a Target Response

The next step is to determine a suitable target response for the system in the ROI. The main goal is to provide a response that is perceived as "flat", and potentially extend the low frequencies (within reasonable limits).

It is important to consider whether the response being optimized is a loudspeaker measurement in an anechoic chamber, or a loudspeaker in a room measured from a listening position. In the former case, the target could be a constant level (with a roll off outside of the loudspeaker's range). In the later case, which is the case for this example, the response at the listening position also depends on the room. The loudspeaker is perceived as "flat" if it exhibits a constant slope down. The degree of that slope depends on factors such as the distance of the listener. The perceived quality also depends on the low frequency cutoff point (approximately -10 dB), so the target curve can be used to extend the low frequencies if the boost remains moderate. Keep in mind that loudspeaker distortion increases in the lower frequencies, and the overall limits of the amplifier and loudspeaker should not be exceeded.

To compute a target response, fit a straight line (in the log-frequency domain) onto the frequency response (in dB) for a subset of the ROI. Add a roll off to the lower frequencies.

% Compute the octave average (only for the region of interest)
[power,cf] = octaveAverage(frequency,db2mag(magnitudeDB),24,ROI,true);
pdb = 20*log10(power) - getMicCalibration(cf);

% Set the range of the linear fit
lfCutOff = 2.5*ROI(1); % lowest frequency to linearize
hfMaxFit = 0.6*ROI(2); % max frequency to fit for

% Fit a straight line in the log-frequency domain
linIdx = cf>=lfCutOff & cf<=hfMaxFit;
if license('test','statistics_toolbox')
    pfit = robustfit(log(cf(linIdx)),pdb(linIdx));
else % if Statistics Toolbox is not available, use a simple linear regression
    pfit = [ones(nnz(linIdx),1) log(cf(linIdx))] \ pdb(linIdx);
end
targetResp = pfit(1)+pfit(2)*log(cf);

% Roll off the low frequencies, starting slightly above the linear range above
lfcutoff = 1.05*lfCutOff;
idx = cf<lfcutoff;
targetResp(idx) = targetResp(idx) - min(30,ROI(1)/2)*((lfcutoff-cf(idx))/lfcutoff).^2;

% Plot the target response
semilogx(frequency, magnitudeDB, 'g:');
hold on
axis([20 24e3 yrange]);
plot(cfFR, pdbFRmic, 'b', LineWidth=2)
plot(cf, targetResp, 'k', LineWidth=2.5);
plot([ROI(1) ROI(1)], yrange, ':m', LineWidth=1.5);
plot([ROI(2) ROI(2)], yrange, ':m', LineWidth=1.5);
title('Target Response for Loudspeaker');
legend('Measured Response (raw data)','Measured Response (octave smoothed)',...
       'Target Response','ROI Low','ROI High',Location='southwest')
xlabel('Frequency (Hz) \rightarrow')
ylabel('Magnitude (dB) \rightarrow')
grid on
hold off

Figure contains an axes object. The axes object with title Target Response for Loudspeaker contains 5 objects of type line. These objects represent Measured Response (raw data), Measured Response (octave smoothed), Target Response, ROI Low, ROI High.

The target response (in black) has a slight downward tilt and a roll off for the low frequencies that allows for some boost (10 to 12 dB).

Parametric Filter Overview

The variables being tuned by the optimization algorithm are typical audio parametric EQ parameters: Center Frequency, Filter Bandwidth, and Peak Gain.

Use the response to produce settings for a 12-band parametric filter (10 peak filters and 2 shelf filters).

Use the designParamEQ function from Audio Toolbox to design the filter. Use the lsqnonlin (Optimization Toolbox) function to perform the fit by tuning the parameters of the EQ bands until the speaker response is as flat as possible.

Before configuring the optimization algorithm, you can look at what a manual filter design looks like for 2 filters. Use the following controls to manually tune the filter parameters and observe the output response using fvtool. This allow us to visualize the parameters that the optimization algorithm automatically tunes.

gain = [4, ...
      -6.6]; 
  
centerFreq = [3520, ...
              7230]; 
          
bandwidth = [1920, ...
              4110];

Next, generate the filter coefficients using the specified parameters by calling designParamEQ:

[B,A] = designParamEQ(2, gain, (centerFreq/(Fs/2)), (bandwidth/(Fs/2)), Orientation='row');

Visualize the filter design.

fvtool([B,A],'Fs',Fs);

{"String":"Figure Figure 1: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB) contains an object of type line.","Tex":"Magnitude Response (dB)","LaTex":[]}

Parametric Filter Optimization

To produce the optimized parametric filters, call a helper function that sets starting values and limits for every filter parameter, then calls lsqnonlin. The optimizer uses the eqObjectiveFct function that computes the response of the given EQ and compares it to the desired response. The lsqnonlin optimizer attempts to minimize that error on every iteration.

% Run the optimization with a selected number of filters
numFilters = 10+2; % 10 peak filters, one low shelf, and one high shelf
[EQ,outputResp] = eqOptimizer(numFilters,frequency,pdb,targetResp,ROI,Fs);
Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to the parallel pool (number of workers: 4).

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0         37        0.162731                         0.425
     1         74       0.0545925             10          0.524      
     2        111       0.0290065             20          0.107      
     3        148       0.0284173        27.6894          0.218      
     4        185       0.0192782        6.92236          0.186      
     5        222        0.015537       0.352966           0.17      
     6        259       0.0116566        13.8447          0.408      
     7        296      0.00883451        13.8447        0.00624      
     8        333      0.00883451        27.6894        0.00624      
     9        370      0.00853245        6.92236          0.102      
    10        407      0.00708075        6.92236         0.0825      
    11        444      0.00627381        13.8447        0.00633      
    12        481      0.00578382        27.6894         0.0409      
    13        518      0.00578382        27.6894         0.0409      
    14        555       0.0051447        6.92236           0.05      
    15        592      0.00494765        13.8447          0.053      
    16        629      0.00470029        7.82709         0.0469      
    17        666      0.00443052        1.95677         0.0821      
    18        703      0.00412292        3.91355         0.0278      
    19        740      0.00376075        7.82709         0.0222      
    20        777      0.00353988        1.91967         0.0206      
    21        814      0.00353988        13.1488         0.0206      
    22        851      0.00353514        3.28721         0.0111      
    23        888      0.00344613       0.821802         0.0153      
    24        925      0.00343122         1.6436        0.00837      
    25        962      0.00339164        1.58048        0.00135      
    26        999      0.00339164        3.28721        0.00135      
    27       1036      0.00339095       0.821802        0.00111      
    28       1073      0.00339053       0.821802        0.00654      
    29       1110      0.00338997       0.205451        0.00349      
    30       1147      0.00338909       0.410901       0.000928      
    31       1184      0.00338825       0.704818        0.00228      
    32       1221      0.00338771       0.176205        0.00082      
    33       1258      0.00338757       0.352409       0.000445      
    34       1295      0.00338743       0.315369       0.000232      
    35       1332      0.00338743       0.352409       0.000232      
    36       1369      0.00338741      0.0881023       0.000206      
    37       1406       0.0033874      0.0881023        0.00027      
    38       1443      0.00338739      0.0881023        0.00013      
    39       1480      0.00338738       0.176205       0.000224      
    40       1517      0.00338736      0.0440511       0.000127      
    41       1554      0.00338735      0.0881023       0.000183      
    42       1591      0.00338734      0.0881023        6.4e-06      
    43       1628      0.00338734      0.0881023       6.56e-05      
    44       1665      0.00338734      0.0881023       2.51e-05      
    45       1702      0.00338733      0.0220256       0.000102      

Optimization stopped because the relative sum of squares (r) is changing
by less than options.FunctionTolerance = 1.000000e-08.

Results

Examine the results. Start by computing the filter frequency response over a larger frequency range.

freqzResp = eqFreqz(EQ,frequency,Fs);
pFR = octaveAverage(frequency,abs(freqzResp),24,fullRange,false);
filtRespFR = 20*log10(pFR);
outputRespFR = pdbFRmic + filtRespFR;

Plot the system response.

% Plot the system response
semilogx(frequency, magnitudeDB, 'g:');
hold on
axis([20 24e3 yrange]);
plot(cfFR, pdbFRmic, 'b', LineWidth=2)
plot(cf, targetResp, 'k', LineWidth=2.5);
plot([ROI(1) ROI(1)], yrange, ':m', LineWidth=1.5);
plot([ROI(2) ROI(2)], yrange, ':m', LineWidth=1.5);
plot(cfFR, outputRespFR, 'r', LineWidth=2);
title('Measured Speaker Response with and without EQ');
legend('Measured Response (raw data)',...
       'Measured Response (octave smoothed)',...
       'Target Response','ROI Low','ROI High',...
       'Corrected Response',Location='southwest')
xlabel('Frequency (Hz) \rightarrow')
ylabel('Magnitude (dB) \rightarrow')
grid on
hold off

Figure contains an axes object. The axes object with title Measured Speaker Response with and without EQ contains 6 objects of type line. These objects represent Measured Response (raw data), Measured Response (octave smoothed), Target Response, ROI Low, ROI High, Corrected Response.

Also plot error relative to target, but invert the error curves to make it easier to compare them with the optimized EQ.

% Plot error relative to target
errorOld = targetResp - pdb;
errorNew = targetResp - outputResp;
semilogx(cf([1,end]), [0 0], 'k', LineWidth=2);
hold on;
plot(cf, errorOld, 'b', LineWidth=2);
plot(cf, errorNew, 'g', LineWidth=2);
plot(cfFR, filtRespFR, 'r', LineWidth=1.5);
hold off; grid on;
axis([20 24e3 -15 15])
title('Error Relative to Target');
xlabel('Frequency (Hz) \rightarrow');
ylabel('Magnitude (dB) \rightarrow');
legend('Target','Original Error (inverted)',...
       'Corrected Error (inverted)',...
       'Optimized EQ',Location='southwest');

Figure contains an axes object. The axes object with title Error Relative to Target contains 4 objects of type line. These objects represent Target, Original Error (inverted), Corrected Error (inverted), Optimized EQ.

Sort peak filters by center frequency and convert to table format.

EQpk = EQ(1:end-2,:);
[~,idx] = sort(EQpk(:,3));
EQ(1:end-2,:) = EQpk(idx,:);

% Create a table with all the EQ settings
filterType =  [repmat("PK",numFilters-2,1);"LS";"HS"];
EQt = table(filterType,EQ(:,3),EQ(:,1),EQ(:,2),VariableNames=...
                   {'Type','Frequency (Hz)','Gain (dB)','Q/S'})
EQt=12×4 table
    Type    Frequency (Hz)    Gain (dB)      Q/S  
    ____    ______________    _________    _______

    "PK"        80.535          7.2637     0.65291
    "PK"        120.32         -8.7422      2.2075
    "PK"        207.13         -7.8682      1.1868
    "PK"        489.96         -9.1706      6.7641
    "PK"        491.03          4.8432      14.672
    "PK"        1098.6          6.3122      4.5365
    "PK"        2034.4         -6.4213      2.1366
    "PK"        6796.3          6.7913      1.1706
    "PK"        9547.3          -9.256      1.2683
    "PK"         12025         -1.5581       5.283
    "LS"          3312          3.0588           5
    "HS"         16571          -2.419       4.166

Compute an overall gain that avoids any tones increasing in amplitude. This is not a guarantee that signals cannot clip but is generally more than sufficient in practice.

gainDB = -max(filtRespFR)-.1;

Instantiate a multibandParametricEQ object with the EQ settings. You can use this object to visualize the filter, create an audio plugin, or load either the object or plugin into the Audio Test Bench.

% Create EQ object
mbpeq = multibandParametricEQ(HasLowShelfFilter=true, HasHighShelfFilter=true, ...
                              NumEQBands=numFilters-2, EQOrder=2, SampleRate=Fs, ...
                              Frequencies=EQ(1:end-2,3)', QualityFactors=EQ(1:end-2,2)', PeakGains=EQ(1:end-2,1)', ...
                              LowShelfCutoff=EQ(end-1,3), LowShelfSlope=EQ(end-1,2), LowShelfGain=EQ(end-1,1), ...
                              HighShelfCutoff=EQ(end,3), HighShelfSlope=EQ(end,2), HighShelfGain=EQ(end,1));

Produce output files to subjectively evaluate the results, either with headphones or in the actual listening room. The IR is included for headphone evaluation but should be omitted when testing in the actual room (which produces that response itself).

[rock,fileFs] = audioread('RockDrums-44p1-stereo-11secs.mp3');

% Resample the test file to match the sample rate of the IR measurement
rock = resample(rock,Fs,fileFs,100);

% Convolve the IR with the test audio. This will simulate the
% effect of the room when evaluating the EQ using headphones.
rockIR = conv(rock(:,1),ir);
rockIR(:,2) = conv(rock(:,2),ir);
rockIR = rockIR*.97/max(abs(rockIR),[],'all');
audiowrite('RockDrums-with-IR.wav',rockIR,Fs,Comment=...
       'Convolution of rock drums with impulse response (for simulated evaluation on headphones)');

Try the plugin in the Audio Test Bench application.

audioTestBench(mbpeq)

In the Audio Test Bench, click the "Audio File Reader" button with the gear icon and select 'RockDrums-with-IR.wav' if using headphones, or simply use 'RockDrums-44p1-stereo-11secs.mp3' if playing back over the same system that was used to measure the impulse response. With the Audio Test Bench, you can toggle the EQ on and off, and you can even further tune the EQ settings to your liking.

You can also process files with the EQ to listen outside of the Audio Test Bench.

% Apply the EQ to the original file (for use in the measured room).
% Use the gain that was computed to avoid clipping.
rockEQ = db2mag(gainDB)*mbpeq(rock);
reset(mbpeq); % reset the EQ before processing another file
audiowrite('RockDrums-with-Correction-only.wav',rockEQ,Fs,Comment=...
          ['Convolution of rock drums with correction (for listening '...
           'in the same environment the IR was measured in)']);

% Apply the EQ to the IR-processed file (for use with headphones).
% Use an output level that avoids clipping.
rockIREQ = mbpeq(rockIR);
reset(mbpeq); % reset the EQ so it is free to process another file
rockIREQ = rockIREQ/max(abs(rockIREQ),[],'all')*.97;
audiowrite('RockDrums-with-IR-and-Correction.wav',rockIREQ,Fs,Comment=...
          ['Convolution of rock drums with impulse response and '...
           'correction (IR simulation for evaluation over headphones)']);

Alternatively, to apply the EQ to any playback on a selected device, Windows users can export the EQ settings in a Room EQ Wizard (REW) compatible format and load it in Equalizer APO.

eqExport2APO("myroomeq.txt",EQ,gainDB);
type("myroomeq.txt")
Filter Settings file

MATLAB Filter Export
Dated: 31-Aug-2022 02:18:11

Notes: Generated using MATLAB example

Equalizer: Generic
Preamp: -8.7 dB
Filter: ON  PK  Fc    80.53 Hz  Gain   7.26 dB  Q 0.65291
Filter: ON  PK  Fc   120.32 Hz  Gain  -8.74 dB  Q 2.20753
Filter: ON  PK  Fc   207.13 Hz  Gain  -7.87 dB  Q 1.18679
Filter: ON  PK  Fc   489.96 Hz  Gain  -9.17 dB  Q 6.76410
Filter: ON  PK  Fc   491.03 Hz  Gain   4.84 dB  Q 14.67229
Filter: ON  PK  Fc  1098.56 Hz  Gain   6.31 dB  Q 4.53645
Filter: ON  PK  Fc  2034.43 Hz  Gain  -6.42 dB  Q 2.13662
Filter: ON  PK  Fc  6796.28 Hz  Gain   6.79 dB  Q 1.17063
Filter: ON  PK  Fc  9547.32 Hz  Gain  -9.26 dB  Q 1.26832
Filter: ON  PK  Fc 12024.72 Hz  Gain  -1.56 dB  Q 5.28302
Filter: ON  LS  Fc  3312.01 Hz  Gain   3.06 dB  Q 1.63270
Filter: ON  HS  Fc 16570.69 Hz  Gain  -2.42 dB  Q 1.46596