Free Books
   Free Books    Physical Audio Signal Processing  

Frequency-Response Matching Using
Digital Filter Design Methods

Given force inputs and velocity outputs, the frequency response of an ideal mass was given in Eq.$ \,$(7.1.2) as

$\displaystyle \Gamma_m(j\omega) \eqsp \frac{1}{m j\omega}, $

and the frequency response for a spring was given by Eq.$ \,$(7.1.3) as

$\displaystyle \Gamma_k(j\omega) \eqsp \frac{j\omega}{k}. $

Thus, an ideal mass is an integrator and an ideal spring is a differentiator. The modeling problem for masses and springs can thus be posed as a problem in digital filter design given the above desired frequency responses. More generally, the admittance frequency response ``seen'' at the port of a general $ N$th-order LTI system is, from Eq.$ \,$(8.3),

$\displaystyle H(s) \isdefs \frac{B(s)}{A(s)} \isdefs \frac{b_M s^M + \cdots b_1 s + b_0}{a_N s^N + \cdots a_1 s + a_0} \protect$ (9.14)

where we assume $ M<N$. Similarly, point-to-point ``trans-admittances'' can be defined as the velocity Laplace transform at one point on the physical object divided by the driving-force Laplace transform at some other point. There is also of course no requirement to always use driving force and observed velocity as the physical variables; velocity-to-force, force-to-force, velocity-to-velocity, force-to-acceleration, etc., can all be used to define transfer functions from one point to another in the system. For simplicity, however, we will prefer admittance transfer functions here.

Ideal Differentiator (Spring Admittance)

Figure 8.1 shows a graph of the frequency response of the ideal differentiator (spring admittance). In principle, a digital differentiator is a filter whose frequency response $ H(e^{j\omega T})$ optimally approximates $ j\omega $ for $ \omega T$ between $ -\pi$ and $ \pi$. Similarly, a digital integrator must match $ 1/j\omega$ along the unit circle in the $ z$ plane. The reason an exact match is not possible is that the ideal frequency responses $ j\omega $ and $ 1/j\omega$, when wrapped along the unit circle in the $ z$ plane, are not ``smooth'' functions any more (see Fig.8.1). As a result, there is no filter with a rational transfer function (i.e., finite order) that can match the desired frequency response exactly.

Figure 8.1: Imaginary part of the frequency response $ H(e^{j\omega T})=j\omega $ of the ideal digital differentiator plotted over the unit circle in the $ z$ plane (the real part being zero).
\includegraphics[scale=0.9]{eps/f_ideal_diff_fr_cropped}

The discontinuity at $ z=-1$ alone is enough to ensure that no finite-order digital transfer function exists with the desired frequency response. As with bandlimited interpolation (§4.4), it is good practice to reserve a ``guard band'' between the highest needed frequency $ f_{\mbox{\tiny max}}$ (such as the limit of human hearing) and half the sampling rate $ f_s/2$. In the guard band $ [f_{\mbox{\tiny max}},f_s/2]$, digital filters are free to smoothly vary in whatever way gives the best performance across frequencies in the audible band $ [0,f_{\mbox{\tiny max}}]$ at the lowest cost. Figure 8.2 shows an example. Note that, as with filters used for bandlimited interpolation, a small increment in oversampling factor yields a much larger decrease in filter cost (when the sampling rate is near $ 2f_{\mbox{\tiny max}}$).

In the general case of Eq.$ \,$(8.14) with $ s=j\omega$, digital filters can be designed to implement arbitrarily accurate admittance transfer functions by finding an optimal rational approximation to the complex function of a single real variable $ \omega $

$\displaystyle H(e^{j\omega}) \eqsp \frac{B(j\omega)}{A(j\omega)} \eqsp \frac{b_... ...ega)^M + \cdots b_1 j\omega + b_0}{a_N (j\omega)^N + \cdots a_1 j\omega + a_0} $

over the interval $ -\omega_{\mbox{\tiny max}}\leq \omega \leq \omega_{\mbox{\tiny max}}$, where $ \omega_{\mbox{\tiny max}}T<\pi$ is the upper limit of human hearing. For small guard bands $ \delta\isdeftext \pi-\omega_{\mbox{\tiny max}}T$, the filter order required for a given error tolerance is approximately inversely proportional to $ \delta$.

Digital Filter Design Overview

This section (adapted from [428]), summarizes some of the more commonly used methods for digital filter design aimed at matching a nonparametric frequency response, such as typically obtained from input/output measurements. This problem should be distinguished from more classical problems with their own specialized methods, such as designing lowpass, highpass, and bandpass filters [343,362], or peak/shelf equalizers [559,449], and other utility filters designed from a priori mathematical specifications.

The problem of fitting a digital filter to a prescribed frequency response may be formulated as follows. To simplify, we set $ T=1$.

Given a continuous complex function $ H(e^{j\omega}),\,-\pi < \omega \le \pi$, corresponding to a causal desired frequency response,9.8 find a stable digital filter of the form

$\displaystyle {\hat H}(z) \isdefs \frac{{\hat B}(z)}{ {\hat A}(z)},$

where
$\displaystyle {\hat B}(z)$ $\displaystyle \isdef$ $\displaystyle {\hat b}_0 + {\hat b}_1 z^{-1} + \cdots + {\hat b}_{\hat{N}_b}z^{-{\hat{N}_b}}$ (9.15)
$\displaystyle {\hat A}(z)$ $\displaystyle \isdef$ $\displaystyle 1 + {\hat a}_1 z^{-1} + \cdots + {\hat a}_{\hat{N}_a}z^{-{\hat{N}_a}} ,$ (9.16)

with $ {\hat{N}_b},{\hat{N}_a}$ given, such that some norm of the error

$\displaystyle J(\hat{\theta}) \isdefs \left\Vert\,H(e^{j\omega}) - {\hat H}(e^{j\omega})\,\right\Vert \protect$ (9.17)

is minimum with respect to the filter coefficients

$\displaystyle \hat{\theta}^K\isdefs [{\hat b}_0,{\hat b}_1,\ldots\,,{\hat b}_{\hat{N}_b},{\hat a}_1,{\hat a}_2,\ldots\,,{\hat a}_{\hat{N}_a}]. $

The filter coefficients are constrained to lie in some subset $ \hat{\Theta}\subset\Re ^{{\hat N}}$, where $ {\hat N}\isdef {\hat{N}_a}+{\hat{N}_b}+1$. The filter coefficients may also be complex, in which case $ \hat{\Theta}\subset{\bf C}^{{\hat N}}$.

The approximate filter $ {\hat H}$ is typically constrained to be stable, and since $ {\hat B}(z)$ is causal (no positive powers of $ z$), stability implies causality. Consequently, the impulse response of the model $ {\hat h}(n)$ is zero for $ n<0$.

The filter-design problem is then to find a (strictly) stable $ {\hat{N}_a}$-pole, $ {\hat{N}_b}$-zero digital filter which minimizes some norm of the error in the frequency-response. This is fundamentally rational approximation of a complex function of a real (frequency) variable, with constraints on the poles.

While the filter-design problem has been formulated quite naturally, it is difficult to solve in practice. The strict stability assumption yields a compact space of filter coefficients $ \hat{\Theta}$, leading to the conclusion that a best approximation $ \hat{H}^\ast $ exists over this domain.9.9Unfortunately, the norm of the error $ J(\hat{\theta})$ typically is not a convex9.10function of the filter coefficients on $ \hat{\Theta}$. This means that algorithms based on gradient descent may fail to find an optimum filter due to their premature termination at a suboptimal local minimum of $ J(\hat{\theta})$.

Fortunately, there is at least one norm whose global minimization may be accomplished in a straightforward fashion without need for initial guesses or ad hoc modifications of the complex (phase-sensitive) IIR filter-design problem--the Hankel norm [155,428,177,36]. Hankel norm methods for digital filter design deliver a spontaneously stable filter of any desired order without imposing coefficient constraints in the algorithm.

An alternative to Hankel-norm approximation is to reformulate the problem, replacing Eq.$ \,$(8.17) with a modified error criterion so that the resulting problem can be solved by linear least-squares or convex optimization techniques. Examples include

  • Pseudo-norm minimization: (Pseudo-norms can be zero for nonzero functions.) For example, Padé approximation falls in this category. In Padé approximation, the first $ {\hat{N}_a}+{\hat{N}_b}+1$ samples of the impulse-response $ h(n)$ of $ H$ are matched exactly, and the error in the remaining impulse-response samples is ignored.

  • Ratio Error: Minimize $ \vert\vert\,H(e^{j\omega})/{\hat H}(e^{j\omega})\,\vert\vert $ subject to $ {\hat B}(z)=1$. Minimizing the $ L2$ norm of the ratio error yields the class of methods known as linear prediction techniques [20,296,297]. Since, by the definition of a norm, we have $ \vert\vert\,e^{j\theta}E(e^{j\omega})\,\vert\vert = \vert\vert\,E(e^{j\omega})\,\vert\vert $, it follows that $ \vert\vert\,H/{\hat H}\,\vert\vert = \vert\vert\,\vert H\vert/\vert{\hat H}\vert\,\vert\vert $; therefore, ratio error methods ignore the phase of the approximation. It is also evident that ratio error is minimized by making $ \vert{\hat H}(e^{j\omega})\vert$ larger than $ \vert H(e^{j\omega})\vert$.9.11 For this reason, ratio-error methods are considered most appropriate for modeling the spectral envelope of $ \vert H(e^{j\omega})\vert$. It is well known that these methods are fast and exceedingly robust in practice, and this explains in part why they are used almost exclusively in some data-intensive applications such as speech modeling and other spectral-envelope applications. In some applications, such as adaptive control or forecasting, the fact that linear prediction error is minimized can justify their choice.

  • Equation error: Minimize

    $\displaystyle \left\Vert\,{\hat A}(e^{j\omega})H(e^{j\omega})-{\hat B}(e^{j\ome... ...}(e^{j\omega})\left[ H(e^{j\omega})-{\hat H}(e^{j\omega})\right]\,\right\Vert. $

    When the $ L2$ norm of equation-error is minimized, the problem becomes solving a set of $ {\hat N}={\hat{N}_a}+{\hat{N}_b}+1$ linear equations.

    The above expression makes it clear that equation-error can be seen as a frequency-response error weighted by $ \vert{\hat A}(e^{j\omega})\vert$. Thus, relatively large errors can be expected where the poles of the optimum approximation (roots of $ {\hat A}(z)$) approach the unit circle $ \vert z\vert=1$. While this may make the frequency-domain formulation seem ill-posed, in the time-domain, linear prediction error is minimized in the $ L2$ sense, and in certain applications this is ideal. (Equation-error methods thus provide a natural extension of ratio-error methods to include zeros.) Using so-called Steiglitz-McBride iterations [287,449,288], the equation-error solution iteratively approaches the norm-minimizing solution of Eq.$ \,$(8.17) for the L2 norm.

    Examples of minimizing equation error using the matlab function invfreqz are given in §8.6.3 and §8.6.4 below. See [449, Appendix I] (based on [428, pp. 48-50]) for a discussion of equation-error IIR filter design and a derivation of a fast equation-error method based on the Fast Fourier Transform (FFT) (used in invfreqz).

  • Conversion to real-valued approximation: For example, power spectrum matching, i.e., minimization of $ \vert\vert\,\vert H(e^{j\omega})\vert^2-\vert{\hat H}(e^{j\omega})\vert^2\,\vert\vert $, is possible using the Chebyshev or $ L-infinity$ norm [428]. Similarly, linear-phase filter design can be carried out with some guarantees, since again the problem reduces to real-valued approximation on the unit circle. The essence of these methods is that the phase-response is eliminated from the error measure, as in the norm of the ratio error, in order to convert a complex approximation problem into a real one. Real rational approximation of a continuous curve appears to be solved in principle only under the $ L-infinity$ norm [373,374].

  • Decoupling poles and zeros: An effective example of this approach is Kopec's method [428] which consists of using ratio error to find the poles, computing the error spectrum $ E=H/{\hat H}$, inverting it, and fitting poles again (to $ 1/E(e^{j\omega})$). There is a wide variety of methods which first fit poles and then zeros. None of these methods produce optimum filters, however, in any normal sense.

In addition to these modifications, sometimes it is necessary to reformulate the problem in order to achieve a different goal. For example, in some audio applications, it is desirable to minimize the log-magnitude frequency-response error. This is due to the way we hear spectral distortions in many circumstances. A technique which accomplishes this objective to the first order in the $ L-infinity$ norm is described in [428].

Sometimes the most important spectral structure is confined to an interval of the frequency domain. A question arises as to how this structure can be accurately modeled while obtaining a cruder fit elsewhere. The usual technique is a weighting function versus frequency. An alternative, however, is to frequency-warp the problem using a first-order conformal map. It turns out a first-order conformal map can be made to approximate very well frequency-resolution scales of human hearing such as the Bark scale or ERB scale [459]. Frequency-warping is especially valuable for providing an effective weighting function connection for filter-design methods, such as the Hankel-norm method, that are intrinsically do not offer choice of a weighted norm for the frequency-response error.

There are several methods which produce $ {\hat H}(z){\hat H}(z^{-1})$ instead of $ {\hat H}(z)$ directly. A fast spectral factorization technique is useful in conjunction with methods of this category [428]. Roughly speaking, a size $ 2{\hat{N}_a}$ polynomial factorization is replaced by an FFT and a size $ {\hat{N}_a}$ system of linear equations.

Digital Differentiator Design

We saw the ideal digital differentiator frequency response in Fig.8.1, where it was noted that the discontinuity in the response at $ \omega=\pm\pi$ made an ideal design unrealizable (infinite order). Fortunately, such a design is not even needed in practice, since there is invariably a guard band between the highest supported frequency $ f_{\mbox{\tiny max}}$ and half the sampling rate $ f_s/2$.

Figure 8.2: FIR differentiator designed by the matlab function invfreqz (Octave). Top: Overlay of the ideal amplitude response ($ \vert j2\pi f\vert$), fitted filter amplitude response, and guard-band limit (at 20 kHz). Bottom: Overlay of ideal phase response ($ \pi /2$ radians), fitted filter phase response, and guard-band limit (20 kHz).
\includegraphics[width=\twidth]{eps/iirdiff-mag-phs-N0-M10}

Figure 8.2 illustrates a more practical design specification for the digital differentiator as well as the performance of a tenth-order FIR fit using invfreqz (which minimizes equation error) in Octave.9.12 The weight function passed to invfreqz was 1 from 0 to 20 kHz, and zero from 20 kHz to half the sampling rate (24 kHz). Notice how, as a result, the amplitude response follows that of the ideal differentiator until 20 kHz, after which it rolls down to a gain of 0 at 24 kHz, as it must (see Fig.8.1). Higher order fits yield better results. Using poles can further improve the results, but the filter should be checked for stability since invfreqz designs filters in the frequency domain and does not enforce stability.9.13

Fitting Filters to Measured Amplitude Responses

The preceding filter-design example digitized an ideal differentiator, which is an example of converting an LTI lumped modeling element into a digital filter while maximally preserving its frequency response over the audio band. Another situation that commonly arises is the need for a digital filter that matches a measured frequency response over the audio band.

Measured Amplitude Response

Figure 8.3 shows a plot of simulated amplitude-response measurements at 10 frequencies equally spread out between 100 Hz and 3 kHz on a log frequency scale. The ``measurements'' are indicated by circles. Each circle plots, for example, the output amplitude divided by the input amplitude for a sinusoidal input signal at that frequency [449]. These ten data points are then extended to dc and half the sampling rate, interpolated, and resampled to a uniform frequency grid (solid line in Fig.8.3), as needed for FFT processing. The details of these computations are listed in Fig.8.4. We will fit a four-pole, one-zero, digital-filter frequency-response to these data.9.14

Figure 8.3: Example measured amplitude-response samples at 10 exponentially spaced frequencies. Circles: Measured amplitude-response points. Solid: Extrapolated, spline-interpolated, and uniformly resampled amplitude response, ready for ifft.
\includegraphics[width=\twidth]{eps/tmps2-G}

Figure: Script (matlab) for simulating a measured amplitude response at 10 exponentially spaced frequencies and extrapolating/interpolating/resampling to obtain a complete, nonparametric amplitude response, uniformly sampled at FFT frequencies. This script generated Fig.8.3.
  
NZ = 1;      % number of ZEROS in the filter to be designed
NP = 4;      % number of POLES in the filter to be designed
NG = 10;     % number of gain measurements
fmin = 100;  % lowest measurement frequency (Hz)
fmax = 3000; % highest measurement frequency (Hz)
fs = 10000;  % discrete-time sampling rate
Nfft = 512;  % FFT size to use
df = (fmax/fmin)^(1/(NG-1)); % uniform log-freq spacing
f = fmin * df .^ (0:NG-1);   % measurement frequency axis

% Gain measurements (synthetic example = triangular amp response):
Gdb = 10*[1:NG/2,NG/2:-1:1]/(NG/2); % between 0 and 10 dB gain

% Must decide on a dc value.
% Either use what is known to be true or pick something "maximally
% smooth".  Here we do a simple linear extrapolation:
dc_amp = Gdb(1) - f(1)*(Gdb(2)-Gdb(1))/(f(2)-f(1));

% Must also decide on a value at half the sampling rate.
% Use either a realistic estimate or something "maximally smooth".
% Here we do a simple linear extrapolation. While zeroing it
% is appealing, we do not want any zeros on the unit circle here.
Gdb_last_slope = (Gdb(NG) - Gdb(NG-1)) / (f(NG) - f(NG-1));
nyq_amp = Gdb(NG) + Gdb_last_slope * (fs/2 - f(NG));

Gdbe = [dc_amp, Gdb, nyq_amp];
fe = [0,f,fs/2];
NGe = NG+2;

% Resample to a uniform frequency grid, as required by ifft.
% We do this by fitting cubic splines evaluated on the fft grid:
Gdbei = spline(fe,Gdbe); % say `help spline'
fk = fs*[0:Nfft/2]/Nfft; % fft frequency grid (nonneg freqs)
Gdbfk = ppval(Gdbei,fk); % Uniformly resampled amp-resp

figure(1);
semilogx(fk(2:end-1),Gdbfk(2:end-1),'-k'); grid('on');
axis([fmin/2 fmax*2 -3 11]);
hold('on'); semilogx(f,Gdb,'ok');
xlabel('Frequency (Hz)');   ylabel('Magnitude (dB)');
title(['Measured and Extrapolated/Interpolated/Resampled ',...
       'Amplitude Response']);

Desired Impulse Response

It is good to check that the desired impulse response is not overly aliased in the time domain. The impulse-response for this example is plotted in Fig.8.5. We see that it appears quite short compared with the inverse FFT used to compute it. The script in Fig.8.6 gives the details of this computation, and also prints out a measure of ``time-limitedness'' defined as the $ L2$ norm of the outermost 20% of the impulse response divided by its total $ L2$ norm--this measure was reported to be $ 0.02$% for this example.

Figure: Desired impulse response obtained from a length 512 inverse FFT of the interpolated measured amplitude response in Fig.8.3.
\includegraphics[width=\twidth]{eps/tmps2-ir}

Note also that the desired impulse response is noncausal. In fact, it is zero phase [449]. This is of course expected because the desired frequency response was real (and nonnegative).

Figure 8.6: Script (matlab) for checking the desired impulse-response for time aliasing. If it wraps around in the time buffer, one can increase the inverse FFT length (Nfft) and/or smooth the desired amplitude response (Sdb).
  
Ns = length(Gdbfk); if Ns~=Nfft/2+1, error("confusion"); end
Sdb = [Gdbfk,Gdbfk(Ns-1:-1:2)]; % install negative-frequencies

S = 10 .^ (Sdb/20); % convert to linear magnitude
s = ifft(S); % desired impulse response
s = real(s); % any imaginary part is quantization noise
tlerr = 100*norm(s(round(0.9*Ns:1.1*Ns)))/norm(s);
disp(sprintf(['Time-limitedness check: Outer 20%% of impulse ' ...
              'response is %0.2f %% of total rms'],tlerr));
% = 0.02 percent
if tlerr>1.0 % arbitrarily set 1% as the upper limit allowed
  error('Increase Nfft and/or smooth Sdb');
end

figure(2);
plot(s,'-k'); grid('on');   title('Impulse Response');
xlabel('Time (samples)');   ylabel('Amplitude');

Converting the Desired Amplitude Response to Minimum Phase

Phase-sensitive filter-design methods such as the equation-error method implemented in invfreqz are normally constrained to produce filters with causal impulse responses.9.15 In cases such as this (phase-sensitive filter design when we don't care about phase--or don't have it), it is best to compute the minimum phase corresponding to the desired amplitude response [449].

As detailed in Fig.8.8, the minimum phase is constructed by the cepstral method [449].9.16

The four-pole, one-zero filter fit using invfreqz is shown in Fig.8.7.

Figure 8.7: Overlay of desired amplitude response (solid) and that of a fourth-order filter fit (dashed) using invfreqz.
\includegraphics[width=\twidth]{eps/tmps2-Hh}

Figure: Script (matlab) for converting the (real) desired amplitude response to minimum-phase form for invfreqz. This script generated Fig.8.7.
  
c = ifft(Sdb); % compute real cepstrum from log magnitude spectrum
% Check aliasing of cepstrum (in theory there is always some):
caliaserr = 100*norm(c(round(Ns*0.9:Ns*1.1)))/norm(c);
disp(sprintf(['Cepstral time-aliasing check: Outer 20%% of ' ...
    'cepstrum holds %0.2f %% of total rms'],caliaserr));
% = 0.09 percent
if caliaserr>1.0 % arbitrary limit
  error('Increase Nfft and/or smooth Sdb to shorten cepstrum');
end
% Fold cepstrum to reflect non-min-phase zeros inside unit circle:
% If complex:
% cf=[c(1),c(2:Ns-1)+conj(c(Nfft:-1:Ns+1)),c(Ns),zeros(1,Nfft-Ns)];
cf = [c(1), c(2:Ns-1)+c(Nfft:-1:Ns+1), c(Ns), zeros(1,Nfft-Ns)];
Cf = fft(cf); % = dB_magnitude + j * minimum_phase
Smp = 10 .^ (Cf/20); % minimum-phase spectrum

Smpp = Smp(1:Ns); % nonnegative-frequency portion
wt = 1 ./ (fk+1); % typical weight fn for audio
wk = 2*pi*fk/fs;
[B,A] = invfreqz(Smpp,wk,NZ,NP,wt);
Hh = freqz(B,A,Ns);

figure(3);
plot(fk,db([Smpp(:),Hh(:)])); grid('on');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Magnitude Frequency Response');
% legend('Desired','Filter');

Further Reading on Digital Filter Design

This section provided only a ``surface scratch'' into the large topic of digital filter design based on an arbitrary frequency response. The main goal here was to provide a high-level orientation and to underscore the high value of such an approach for encapsulating linear, time-invariant subsystems in a computationally efficient yet accurate form. Applied examples will appear in later chapters. We close this section with some pointers for further reading in the area of digital filter design.

Some good books on digital filter design in general include [343,362,289]. Also take a look at the various references in the help/type info for Matlab/Octave functions pertaining to filter design. Methods for FIR filter design (used in conjunction with FFT convolution) are discussed in Book IV [456], and the equation-error method for IIR filter design was introduced in Book II [449]. See [281,282] for related techniques applied to guitar modeling. See [454] for examples of using matlab functions invfreqz and invfreqs to fit filters to measured frequency-response data (specifically the wah pedal design example). Other filter-design tools can be found in the same website area.

The overview of methods in §8.6.2 above is elaborated in [428], including further method details, application to violin modeling, and literature pointers regarding the methods addressed. Some of this material was included in [449, Appendix I].

In Octave or Matlab, say lookfor filter to obtain a list of filter-related functions. Matlab has a dedicated filter-design toolbox (say doc filterdesign in Matlab). In many matlab functions (both Octave and Matlab), there are literature citations in the source code. For example, type invfreqz in Octave provides a URL to a Web page (from [449]) describing the FFT method for equation-error filter design.

Next Section:
Commuted Synthesis
Previous Section:
Modal Expansion