Sparse Graphic Equalizer Design

—A typical graphic equalizer frequency resolution is one-third octave comprising 31 bands. A previous design based on a least-squares optimization of the band-ﬁlter gains with a single second-order section per band has an accuracy of 1 dB. However, the design always uses all the band ﬁlters even when a small number of gains is adjusted. This letter proposes a sparse design of a one-third-octave graphic equalizer, where the number ofactivebandsisminimizedusingorthogonalmatchingpursuitandlinearprogrammingbeforetheleast-squaresgainoptimization.Inaddition,thebandwidthsofthehigh-frequencybandﬁltersaregaindependentinordertooptimizetheirshaperelativetoananalogprototype.Theoptimalbandwidthisobtainedwithlinearinterpolationduringtheﬁlterdesign.Theproposeddesignachievesapproximatelythesameorbetteraccuracyincomparisontothe state-of-the-art non-sparse design and can be used to automatically reducethecomputationalloadofequalizationwhenonlysomebandﬁltersareactive.


I. INTRODUCTION
G RAPHIC equalizers (GEQs) that are commonly used in audio technology consist of parametric filters boosting or attenuating the signal gain in predefined frequency bands [1]- [3]. State-of-the-art GEQs [4], [5] require one second-order filter per band to achieve an accuracy of ±1 dB from the target gain values, which is sufficient for human hearing [6], [7]. This letter presents a sparse GEQ design method that reduces the number of band filters when some frequency band gains are unused, requiring neither boost nor cut.
Early GEQs were built using analog band filters connected in parallel [2]. Digital GEQ designs are based on either the parallel [8]- [10] or cascade infinite impulse response (IIR) filter design [11]- [13]. In 2017, Välimäki and Liski [4] proposed an accurate method for designing a cascade octave-band GEQ with 1-dB accuracy using the least squares (LS) method. Their design has subsequently been extended to one-thirdoctave design [5] and to a parallel structure [14]. Recently may be removed and the bandwidths of those remaining are increased [15]. This paper focuses on a common special case in which one or more band gains are 0 dB. This means that the corresponding filter may not be needed, unless the band filter must compensate for gain leakage from neighboring bands. An inactive band filter can be replaced with unity gain, which leads to computational savings during real-time filtering. We propose using a combination of matching pursuit (MP) and linear programming algorithms to optimize the filter gains, so that the sparsest possible solution fulfilling the accuracy criterion is obtained. The proposed design is shown to be at least as accurate as the non-sparse baseline method.
This letter is structured as follows. Section II recapitulates the baseline design, Section III introduces the new sparse approximation method, Section IV compares the proposed and baseline methods, and Section V concludes the letter.

II. BACKGROUND
An accurate cascade GEQ design with one-third-octave resolution [5] achieves a high accuracy with a single second-order IIR filter per band by allowing the filter gains to differ from the target gains, thus accounting for the leakage of the filters to their neighboring bands. This is accomplished with the LS design based on an interaction matrix B: where g are the optimized filter gains, vector t 1 contains the user-set gains in dB in odd rows and their linearly interpolated values in even rows, and T is the transpose operation [5]. The elements of the M × N interaction matrix B, where N is the number of bands and M = 2N − 1 is the number of design points, are obtained by designing each of the N band filters with a prototype gain g p and evaluating their dB responses at M points comprising the filter center frequencies and the geometric means between neighboring center frequencies. Each dB response is normalized by its prototype gain to obtain the elements of B, which are real-valued numbers between 0 and 1.
The design method [5] utilizes symmetric IIR band filters proposed by Orfanidis [16]. The symmetry is achieved by allowing the Nyquist gain (the gain at the Nyquist limit) to be arbitrary instead of 0 dB. The following equations include the Nyquist gain G Nq as a free parameter, instead of it being determined during the filter design based on other input values [13]. The transfer function of each band filter is where the gain factor k 0 = V (G Nq + W 2 + A 2 ) and the filter coefficients are and where G is the filter linear peak gain, G B is the gain at the band edges, B = 2πf B /f s is the bandwidth, ω c = 2πf c /f s is the center frequency, and f s = 44.1 kHz is the sampling rate.
In the one-third-octave design, the number of bands is N = 31. The recommended center frequencies are obtained by starting from 1 kHz and moving equidistantly on a logarithmic scale in both directions. In order to have accurate control of the filters at the band center frequencies, the bandwidth boundaries are set equal to the neighboring center frequencies, which results in the nominal values B = ( For improved accuracy, gain-dependent bandwidths are an option, i.e., each band filter is designed with varying peak gains (e.g., 1-33 dB; only the positive gains are needed due to symmetry), and the shapes of the filters are compared to greatly oversampled versions (f s,O = 10 MHz). By minimizing the maximum absolute difference between the two filters, we obtain a factor K B that determines the bandwidth: B = K B ω c . The factors K B are stored in a table, and the bandwidths are obtained by linearly interpolating the tabulated factor values.
The design equations in (4) enable control over G B . It is beneficial to link the gain at the band edges to the filter peak gain g B = cg, where g B and g are the respective values of G B and G in dB and c is a free design parameter. Rämö et al. [5] show that c = 0.38 results in an accurate design.
The Nyquist gain G Nq , also being a design parameter, and because its value is gain-dependent and differs for each band filter, can be determined using polynomial approximation [5]. However, for improved accuracy the gain-dependent Nyquist gain for each band can be obtained from the highly oversampled filters by evaluating their response at 22.05 kHz, i.e., at the Nyquist limit for the sample rate of 44.1 kHz. These values are then stored in a table, and during filter design, the correct Nyquist gains are approximated with linear interpolation.
The two final design parameters are the prototype gain g p and the number of iterative steps. For fast GEQ design, the number of iterative steps can be set to zero. The resulting inaccuracy can be decreased with the band-dependent g p . In this work, the iterative method is used, and thus the value of the prototype gain is not as critical; g p = 11 dB is applied [5]. The initial interaction matrix B is used to find a sparse approximation of the desired magnitude response. Finally, step (1) is performed using the second interaction matrix B 2 , obtained with the first set of optimized gains g opt , which now act as band-dependent prototype gains. A column of B 2 can be zero if the optimal gain for the corresponding band is zero.

III. SPARSE APPROXIMATION OF THE TARGET RESPONSE
In this section, two design procedures are presented: the first based on MP and the second based on a linear programming model. Thereafter, a novel design procedure combining the benefits of the two approaches is presented.

A. Sparse Approximation With Matching Pursuit
In [5], the gains of the filter cascade are optimized using an interaction matrix B ∈ R M ×N , where M ≥ N . As a result, the final magnitude response is represented as a linear combination of the columns of the interaction matrix, where the gains are the coefficients of the summation. This section presents our proposed solution to the problem of finding a vector of coefficients that is the sparsest approximation of the target response.
The term sparse refers to a measurable property of a vector. In particular, the sparsity of a vector is measured using the zero-norm, which is a nonlinear function counting the number of non-zero entries in the vector. The zero-norm of a vector x is defined as x 0 = card{x i : x i = 0}. Finding a sparse approximation of the target dB gains can be formulated as a constrained nonlinear optimization problem: where t 1 ∈ R M is a vector containing the user-set gains in dB in odd rows and their linearly interpolated values in even rows [5], and ξ is the acceptable error for the approximation. Problem (5) is a combinatorial optimization problem, it is defined NP-Hard [17], and is a general problem that arises in various fields such as machine learning [18] and signal processing [19], [20]. Although MP algorithms [21] typically solve problems in the form min x x 0 , s.t. Ax = b, where the matrix A is a full rank over-complete matrix (M ≤ N ), the pursuit algorithm can be modified to allow error tolerance and therefore address the cases where the dictionary matrix is square or over-determined (in this case some properties of MP, such as the uniqueness of the solution, do not apply). In the following, a special pursuit algorithm, named Orthogonal Matching Pursuit (OMP) [22], is presented (see Algorithm 1).

B. Sparse Approximation With Linear Programming
Problem (5) may have no solution; to find the sparsest approximation of the target response, (5) can be rewritten as where λ 1 is a penalty term used to favor accurate solutions. To make the problem solvable with a standard optimization algorithm, the zero-norm, which is a non-convex discontinuous function, can be replaced by the norm 1 := The minimization of the 1 over convex set min x x 1 , s.t. x ∈ P , can be modeled as a linear programming problem min x,y i y i , s.t. x ∈ P, −y i ≤ x i ≤ y i . Thus, the following problem is obtained: where h i are slack variables required for modeling the 1 minimization problem as a linear programming problem. Linear optimization algorithms usually require that the problem is presented in the following standard form: min x c t x, s.t. Ax = b, x ≥ 0. Here, problem (7) is transformed to the standard form, and a practical optimization algorithm suitable for solving the problem is presented. Let where 1 N is a size N vector of ones and 0 N is a size N vector of zeros. The inequalities in (7) can then be written in matrix form:  (8) and (10) 2 Solve the problem (10) by using the Mehrotra's algorithm, let x * be the solution where 0 M ×N is a size M × N matrix of zeros and I N is the identity matrix of size N . To achieve the standard form, all the variables are constrained to be non-negative. For this purpose, the variables z are split into two sets, the positives and the negatives, i.e., z = z + − z − , where z + = max(z, 0) ≥ 0, z − = max(−z, 0) ≥ 0. Problem (7) is rewritten in the following intermediate form: The inequality constraints are converted to equalities by introducing a vector of slack variables s constrained to be nonnegative. As a result the inequalities in (7) are transformed to where p = 4N + 2, q = 2M + 2N + 1, Mehrotra's practical algorithm [23] has been implemented in MATLAB to solve the optimization problem (10) (see Algorithm 2).

C. Proposed Design Procedure
Algorithm 1 is a greedy algorithm that finds a sparse approximation by using an incremental approach until the required accuracy (ξ) is achieved or the maximum number of iterations have been performed. As a result, the solution may not be feasible. On the other hand, Algorithm 2 finds the sparsest solution that minimizes the approximation error of the desired user-set gains with an increased computational cost. In fact, Algorithm 2 requires the definition of additional slack variables to achieve the standard form of the linear programming problem, and as a consequence, the solution space of Algorithm 2 is higher than that of Algorithm 1.
Thus, in order to utilize the benefits of both algorithms (see [19] for a review of sparse approximation algorithms), a new procedure is proposed, which calls the first two algorithms as sub-algorithms: Algorithm 1 is first used to obtain the sparse solution. If the result is a feasible solution, g opt is obtained, and the algorithm stops. In the case of an unfeasible solution, Algorithm 2 is utilized. The algorithm stops on obtaining g opt . At the end of the proposed algorithm, the equalizer band gains are calculated with (1), using the interaction matrix B 2 computed with the band-dependent prototype gain g opt . As a result, some columns of the matrix can be null, and in such cases a sparse solution has been found, allowing the corresponding band filters to be omitted from the filter chain.

IV. COMPARISON
This section compares the novel design to a baseline stateof-the-art GEQ design by Rämö et al. [5]. The first clear advantage is the reduced computing time: since the sparse approximation can lead to unused bands, less filters are used during the equalization in the proposed design compared to the baseline GEQ, where all bands are always active. Fig. 1 shows the mean processing time of a 30-s pink-noise segment filtered using a GEQ with an increasing number of active gains (at 3 dB) compared to the baseline with all bands active. The mean computational time (over 100 runs) is reported as a function of the number of active bands; the number of active bands typically differs from the number of user-set gains, and not all values are necessarily obtained, as shown in the figure. The computing time is seen to increase approximately linearly with the number of active band filters. Thus, computational savings are obtained when some band filters are not used.
The accuracy of the GEQ is paramount. Here, an application example is presented to compare the proposed design to the state-of-the-art GEQ [5]. Suitable application examples are those in which a limited bandwidth is equalized, such as low-frequency equalization in noise [24] and headphone equalization [25]. The latter is considered here: we take a measured response from a pair of Sennheiser HD-650 headphones and set the target to the Harman target curve [26].
The measured headphone response is shown in Fig. 2(a) together with the target curve. Both responses are normalized so that they are equal to 0 dB at 1 kHz. As is shown, the largest differences between the two curves are at low and high frequencies. However, these frequency regions are left unprocessed and the correction is limited to 100 Hz-10 kHz for the following reasons. At frequencies below 100 Hz, the headphone response is highly attenuated due to the mechanical limitations of the headphones. Boosting the gain would introduce unwanted distortion. Above 10 kHz, the measurement uncertainty caused by variations in the headphone placement increases, and boosting high frequencies leads to coloration artifacts [27]. Furthermore, dummy heads used for measuring headphones simulate the behavior of a normal adult human ear in the frequency range 100 Hz-10 kHz [28], [29], and thus the measurement may be erroneous at low and high frequencies.
The target gains are obtained by sampling the onethird-octave smoothed difference of the measured headphone response and the target curve at the GEQ center frequencies between 100 Hz-10 kHz. Fig. 2(b) shows the target gains and the realized GEQ responses, where the proposed method uses an accuracy parameter value of ξ = 0.2 dB and λ = 1000. Both GEQ responses match the target very accurately, since their maximum errors at any band filter center frequency between 100 Hz-10 kHz are 0.20 dB, which occurs at 8 kHz for the baseline design, and 0.13 dB at 315 Hz for the proposed sparse design. The state-of-the-art GEQ uses all 31 band filters, whereas the proposed design only has 17 active bands.
Thus, the proposed design maintains the desired accuracy using 45% less band filters. Note that the inactive bands need not be at the edges of the equalization band, but they can be in the middle as well, as shown in this example. A similar behavior is observed also with other tested target gain settings, where the proposed method results in approximately the same or improved accuracy compared to the state-of-the-art design.

V. CONCLUSION
This letter proposes a new GEQ design utilizing sparse approximation. The one-third-octave GEQ comprises a single second-order section per band, and the band gains are optimized with an LS solution. The proposed method minimizes the number of active bands using a sparse approximation algorithm utilizing either OMP or linear programming. The band filters correspoding to the non-active gains can be left out of the cascade GEQ. The accuracy of the design is increased with gaindependent bandwidth and Nyquist gain values, and, depending on the application, the proposed design results in approximately the same or improved accuracy when compared to a previous state-of-the-art GEQ, achieving a ±1 dB accuracy while potentially using fewer bands. The proposed design is suitable in all applications where a limited bandwidth is processed or some filter gains are 0 dB. The relevant MATLAB code is available online [30].