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.

Start the Impulse Response Measurer app by entering impulseResponseMeasurer at the command prompt. You can also click the app icon on the Apps tab of the MATLAB® Toolstrip.

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.

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;

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.

hold on
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

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 with no output arguments plots the microphone frequency response.


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.

hold on
plot([ROI(1) ROI(1)],yrange,':',Color="#A2142F",LineWidth=1.5);
plot([ROI(2) ROI(2)],yrange,':',Color="#A2142F",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

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 this toolbox is not available, use a simple linear regression
    pfit = [ones(nnz(linIdx),1) log(cf(linIdx))] \ pdb(linIdx);
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
hold on
axis([20 24e3 yrange]);
plot([ROI(1) ROI(1)],yrange,':',Color="#A2142F",LineWidth=1.5);
plot([ROI(2) ROI(2)],yrange,':',Color="#A2142F",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

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, ...
centerFreq = [3520, ...
bandwidth = [1920, ...

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.


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);
                                            Norm of      First-order 
 Iteration  Func-count      Resnorm            step       optimality
     0         37          0.162731                            0.425
     1         74         0.0545925              10            0.524      
     2        111         0.0289868              20            0.108      
     3        148         0.0272225         27.8857            0.215      
     4        185         0.0187429         6.97142            0.188      
     5        222         0.0154992        0.345245            0.174      
     6        259          0.014701         13.9428            0.302      
     7        296         0.0105981         3.48571            0.195      
     8        333        0.00902063         6.97142           0.0335      
     9        370        0.00788068         6.97142          0.00907      
    10        407        0.00635659         13.9428            0.359      
    11        444        0.00635659         30.7718            0.359      
    12        481         0.0055785         6.97142           0.0129      
    13        518        0.00492435         13.9428          0.00426      
    14        555        0.00452565         27.8857           0.0686      
    15        592        0.00425527         12.6583           0.0446      
    16        629        0.00425527         16.7565           0.0446      
    17        666         0.0040533         4.18912           0.0339      
    18        703        0.00384572         8.37823          0.00953      
    19        740        0.00384572         4.72211          0.00953      
    20        777        0.00374876         1.18053           0.0102      
    21        814        0.00357736         2.36106          0.00508      
    22        851        0.00357736         4.41566          0.00508      
    23        888        0.00351695         1.10391           0.0078      
    24        925        0.00334974         2.20783            0.012      
    25        962        0.00334974         4.41566            0.012      
    26        999        0.00323091         1.10391          0.00812      
    27       1036        0.00309264         2.20783          0.00466      
    28       1073        0.00285452         2.20783           0.0865      
    29       1110        0.00262239         4.01131           0.0678      
    30       1147        0.00262239         4.33015           0.0678      
    31       1184        0.00255327         1.08254           0.0425      
    32       1221        0.00245483         2.40821           0.0338      
    33       1258        0.00243857         4.33015          0.00585      
    34       1295        0.00243857         3.19387          0.00585      
    35       1332         0.0024362          0.2198          0.00533      
    36       1369        0.00243225        0.798467          0.00404      
    37       1406        0.00243037        0.798467          0.00285      
    38       1443        0.00242846         1.59693          0.00203      
    39       1480        0.00242656         1.59693          0.00116      
    40       1517        0.00242605         1.59693         0.000984      
    41       1554        0.00242549        0.864994         0.000378      
    42       1591        0.00242537         1.26553          0.00064      
    43       1628        0.00242512        0.795897         0.000308      
    44       1665        0.00242512          1.5253         0.000308      
    45       1702        0.00242505        0.381326         0.000439      
    46       1739        0.00242496        0.762652         0.000113      
    47       1776        0.00242496          1.5253         0.000113      
    48       1813        0.00242492        0.381326         0.000196      
    49       1850        0.00242484        0.762652         9.39e-05      
    50       1887        0.00242484          1.7062         9.39e-05      
    51       1924         0.0024248        0.381326         0.000186      
    52       1961        0.00242474        0.762652         0.000158      
    53       1998        0.00242474          1.2928         0.000158      
    54       2035        0.00242471        0.323201         0.000187      
    55       2072        0.00242466        0.646402         0.000105      
    56       2109        0.00242466          1.1841         0.000105      
    57       2146        0.00242464        0.296026         0.000128      
    58       2183         0.0024246        0.592051         8.08e-05      
    59       2220        0.00242459        0.973815         0.000327      
    60       2257        0.00242454        0.243454         0.000198      
    61       2294        0.00242453        0.486907          8.3e-05      
    62       2331        0.00242452        0.372546         4.24e-05      

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


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
hold on
axis([20 24e3 yrange]);
plot([ROI(1) ROI(1)],yrange,':',Color="#A2142F",LineWidth=1.5);
plot([ROI(2) ROI(2)],yrange,':',Color="#A2142F",LineWidth=1.5);
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

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;
hold on;
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');

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"         72.52          3.9074      2.1516
    "PK"        122.84         -5.8932      3.0223
    "PK"        211.61          -9.727     0.92132
    "PK"        492.18           -6.61      3.7217
    "PK"           779         -4.2525      7.8339
    "PK"        1102.3          5.1899       4.358
    "PK"        2127.6         -15.942      1.1218
    "PK"        4842.4          8.0688     0.48624
    "PK"        9175.5         -6.1896      2.0737
    "PK"         11710         -6.4538      2.1912
    "LS"        3194.1          6.7385           5
    "HS"         15580         -3.5357      3.4374

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');
       'Convolution of rock drums with impulse response (for simulated evaluation on headphones)');

Type audioTestBench(mbpeq) at the command prompt to try the plugin in the Audio Test Bench app.

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
          ['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;
          ['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.

Filter Settings file

MATLAB Filter Export
Dated: 21-Dec-2022 18:26:34

Notes: Generated using MATLAB example

Equalizer: Generic
Preamp: -8.5 dB
Filter: ON  PK  Fc    72.52 Hz  Gain   3.91 dB  Q 2.15162
Filter: ON  PK  Fc   122.84 Hz  Gain  -5.89 dB  Q 3.02225
Filter: ON  PK  Fc   211.61 Hz  Gain  -9.73 dB  Q 0.92132
Filter: ON  PK  Fc   492.18 Hz  Gain  -6.61 dB  Q 3.72172
Filter: ON  PK  Fc   779.00 Hz  Gain  -4.25 dB  Q 7.83389
Filter: ON  PK  Fc  1102.35 Hz  Gain   5.19 dB  Q 4.35797
Filter: ON  PK  Fc  2127.60 Hz  Gain -15.94 dB  Q 1.12182
Filter: ON  PK  Fc  4842.40 Hz  Gain   8.07 dB  Q 0.48624
Filter: ON  PK  Fc  9175.50 Hz  Gain  -6.19 dB  Q 2.07367
Filter: ON  PK  Fc 11710.49 Hz  Gain  -6.45 dB  Q 2.19118
Filter: ON  LS  Fc  3194.10 Hz  Gain   6.74 dB  Q 1.89623
Filter: ON  HS  Fc 15580.23 Hz  Gain  -3.54 dB  Q 1.34551