Notes regarding invfreqsmethod listing in Fig.19:

• The basic method is to use invfreqs to find the coefficients of a second-order analog filter having a frequency response close to what we measured. Then, we calculate the $Q$ and resonance frequency $\omega_r$ from the simple formulas relating these quantities to the transfer-function coefficients:
  
  $$H(s) = \frac{s-\xi}{\left(\frac{s}{\omega_r}\right)^2 + \frac{2}{Q} \left(\frac{s}{\omega_r}\right) + 1}$$
  
  
  where $\xi$ is an arbitrary finite zero location near dc. Finally, using the impulse invariant method,²⁰we map the estimated $Q$ and $\omega_r$ to digital biquad coefficients, as shown in the code.

• A useful quantitative measure of filter approximation error can be defined as follows (insert after the call to freqs):

    err = norm(wt(:) .* (db(Hp(:))-db(Hph(:))))/norm(wt(:) .* Hp(:));
    disp(sprintf(['Relative weighted L2 norm of frequency ',
                  'db-magnitude response error = %f'],err));
  

• The minphaseir function for converting a spectrum to its minimum-phase counterpart is listed in Fig.20 and discussed in [12].²¹This is a time-domain version of mps.m from the RTFMT matlab code.²² The support routines clipdb and fold are similarly listed and discussed in [12] (and included with the RTFMT matlab code).

  Conversion to minimum phase is an important preprocessing step for phase-sensitive filter-design methods when the desired frequency response is given as the FFT of a measured impulse response with an unknown excess delay. Otherwise, the leading zeros can be trimmed away manually. Further discussion on this point appears in [11].

  Note that minphaseir will give poor results (in the form of time-aliasing) if there are poles or zeros too close to the unit circle in the $z$ plane. To address this, the spectrum of h can be smoothed to eliminate any excessively sharp peaks or nulls. The requirement is that the inverse-FFT of the log magnitude spectrum H = fft(h) must not time-alias appreciably at the FFT size used (which is determined by the smoothness and the amount of zero padding used in the FFT).

• The approximation $z = e^{sT} \approx 1 + sT$ assumes the highest resonance frequency (measured to be 2.2 kHz here--see Fig.16) is much less than the sampling rate $f_s$. Since we will use sampling rates no lower than $f_s=44.1$ kHz, and since the resonance frequency does not need to be exact, this approximation is adequate for wah-pedal simulation.

Figure 20: Listing of matlab function minphaseir for converting an impulse response to its minimum-phase counterpart.

 

function [hmp] = minphaseir(h) 
%
% MINPHASEIR - Convert a real impulse response to its 
%              minimum phase counterpart
% USAGE:
%       [hmp] = minphaseir(h) 
% where
%   h   = impulse response (any length - will be zero-padded)
%   hmp = min-phase impulse response (at zero-padded length)

nh = length(h);
nfft = 2^nextpow2(5*nh);
Hzp = fft(h,nfft);
Hmpzp = exp( fft( fold( ifft( log( clipdb(Hzp,-100) )))));
hmpzp = ifft(Hmpzp);
hmp = real(hmpzp(1:nh));

In summary, a second-order analog transfer-function was fit to each RTFMT-measured frequency response using invfreqs in Octave. Closed-form expressions relating the returned coefficients to $Q$, peak-frequency, and peak-gain were used to obtain these parameters.