Fitting Filters to Measured Amplitude Response Data Using invfreqz in Matlab

Julius Orion Smith III November 8, 20107 comments Coded in Matlab

In the 1980s I wrote the matlab function invfreqz that implements a fast (FFT-based) equation-error method for recursive filter design given samples of the desired frequency response. This method designs causal filters, fitting both the phase and magnitude of the desired (complex) frequency response. As a result, one must be careful about the specified phase response, especially when it doesn't matter!  This article gives some pointers on this topic.

The case I see most often is where there are no phase data at all. For example, a filter response may be specified simply as a desired (real) gain at various frequencies, i.e., only the amplitude response is specified. The easiest choice of phase response in this case is zero, but that would correspond to an impulse response that is symmetric about time zero, and causal filters cannot produce any response before time zero. In this case, phase-sensitive filter-design methods such as invfreqz will generally give their best results for minimum phase in place of zero phase. In other words, we need to convert our desired amplitude response to the corresponding minimum-phase frequency response (whose magnitude equals the original desired amplitude response).

A commonly used method for computing the minimum-phase frequency response from its magnitude is the cepstral method.  The steps are as follows:

  1. Interpolate the amplitude response samples from 0 to half the sampling rate, if necessary, and resample to a uniform "FFT frequency axis", if necessary.  Denote the real, sampled amplitude response by S(k).
  2. Perform an inverse FFT of log(S) to obtain the real cepstrum of s, denoted by c(n).
  3. Fold the noncausal portion of c(n) onto its causal portion.
  4. Perform a forward FFT, followed by exponentiation to obtain the minimum phase frequency response Sm(k), where now sm(n) is causal, and |Sm(k)|=S(k).

These steps are illustrated by the matlab example below (tested in Octave, version 3.2.4, on Mac OS X).  Copy-paste it and give it a try!

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

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']);

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');

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

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');
% 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);

plot(fk,db([Smpp(:),Hh(:)])); grid('on');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Magnitude Frequency Response');