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.
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,':',Color="#77AC30") hold on plot(cfFR,pdbFR,'b',LineWidth=2,Color="#0072BD") 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.
getMicCalibration()
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 optimizationPlot the corrected response and the ROI region.
semilogx(frequency,magnitudeDB,':',Color="#77AC30") hold on plot(cfFR,pdbFRmic,Color="#0072BD",LineWidth=2) 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); 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,':',Color="#77AC30"); hold on axis([20 24e3 yrange]); plot(cfFR,pdbFRmic,Color="#0072BD",LineWidth=2) plot(cf,targetResp,Color="#000000",LineWidth=2.5); 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, ...
-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)
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.
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,':',Color="#77AC30"); hold on axis([20 24e3 yrange]); plot(cfFR,pdbFRmic,Color="#0072BD",LineWidth=2) plot(cf,targetResp,Color="#000000",LineWidth=2.5); plot([ROI(1) ROI(1)],yrange,':',Color="#A2142F",LineWidth=1.5); plot([ROI(2) ROI(2)],yrange,':',Color="#A2142F",LineWidth=1.5); plot(cfFR,outputRespFR,Color="#D95319",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
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],Color="#000000",LineWidth=2); hold on; plot(cf,errorOld,Color="#0072BD",LineWidth=2); plot(cf,errorNew,Color="#77AC30",LineWidth=2); plot(cfFR,filtRespFR,Color="#D95319",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');
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'); audiowrite('RockDrums-with-IR.wav',rockIR,Fs,Comment=... '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 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: 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
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)