Notes regarding the matlab script in Fig.18:

• The del parameter specifies the number of initial samples to skip over in the input sound (the round-trip delay of the measurement system). It can be a very critical parameter if the filter-design method is sensitive to frequency-response phase. Since invfreqsmethod converts the measured frequency response to minimum phase, it is not sensitive to del over a wide range.
• The frequency-zoom interval f1z to f2z should include the resonance peak. Since the measured response is noisier at very low and very high frequencies, frequency-zooming improves the resulting match.
• The estimated Q values printed at the end are
  
  $$Q = [9.4, 4.0, 1.9],$$
  
  
  and the estimated pole frequencies are
  
  $$f_p = [464, 838, 2252]\quad\mbox{Hz.}$$
  
  
  While these estimates could be improved by various means [8], they appear to be more than sufficiently accurate for the application at hand.
• The plot overlays in Figures 14 through 16 are plots of the returned frequency responses Hp (measured) and Hd (the model) versus sampled radian frequency w.

Figure 19: Listing of matlab code for estimating the Q and pole-frequency in the measured frequency response of a second-order resonator (wah pedal).

 

function [Q,wp,Hp,Hd,w] = invfreqsmethod(h,f1,f2,fs);
%
% INVFREQSMETHOD - use invfreqs to estimate Q and resonance 
%               frequency from a resonator impulse response.
% USAGE:
%       [Q,wp,w,Hp,Hd] = invfreqsmethod(h,f1,f2,fs);
% where
%   h = impulse response (power of 2 length preferred)
%   f1 = lowest frequency of interest (Hz)
%   f2 = highest frequency of interest (Hz)
%   fs = sampling rate (Hz)
%   Q  = estimated resonator Quality factor
%   wp = estimated pole frequency (rad/sec)
%
% invfreqs and invfreqz are VERY sensitive to where time 0 
% is defined.  Normalize this by converting the impulse 
% response to its minimum-phase counterpart:
h = minphaseir(h);      H = fft(h);

L0 = length(h);
k1  = round(L0*f1/fs);  k2  = round(L0*f2/fs);
Hp = H(k1:k2);          L0p = length(Hp)

w = 2*pi*[k1:k2]*fs/L0;
wt = 1 ./ w.^2; % Nominal weighting proportional to 1/freq
[Bh,Ah] = invfreqs(Hp,w,2,2,wt);
Hph = freqs(Bh,Ah,w);

% Denominator to canonical form A(s) = s^2 + (wp/Q) s + wp^2:
Ahn = Ah/Ah(1); wp = sqrt(Ahn(3)); Q = wp/Ahn(2);
fp = wp/(2*pi);
disp(sprintf('Pole frequency = %f Hz',fp));
disp(sprintf('Q = %f',Q));

kp = L0*fp/fs - k1 + 1; % bin number of pole (+1 for matlab)
k1z = max(1,floor(kp*0.5));
k2z = min(length(w),ceil(kp*1.5));
ndx = [k1z:k2z];

% BiQuad fit using z = exp(s T) ~ 1 + sT approximation:
frn = fp/fs; % Normalized pole frequency (cycles per sample)
R = 1 - pi*frn/Q; % pole radius
theta = 2*pi*frn; % pole angle
A = [1, -2*R*cos(theta), R*R];
B = [1 -1]; % zeros guessed by inspection of amp response
Hd = freqz(B,A,w/fs);
gain = max(abs(Hp))/max(abs(Hd))
B = B*gain;  Hd = Hd*gain;