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


Good description and comments!
9 years ago
Sorry, you need javascript enabled to post any comments.
Rick Lyons
Hello JOS, Thanks for posting this very interesting code. I'm going to experiment with it. It looks like the line: Ns = length(Gdbfk); if Ns!=Nfft/2+1, error("confusion"); end should be: Ns = length(Gdbfk); if Ns~=Nfft/2+1, error('confusion'); end (I imagine those three "typos" occurred because of some sort of software text-translation problem when when your code was posted as a blog.) Thanks again JOS, [-Rick-]
9 years ago
Sorry, you need javascript enabled to post any comments.
Fixed - thanks! - jos
9 years ago
Sorry, you need javascript enabled to post any comments.
Dr. Smith or Rick (Hi), Unstable impulse reponse seems to be what I am getting, with a largely growing plot with impz of the transfer function derived with invfreqz? If this is true, and I think it is; then I am making a mistake? I am using the mag and phase data file from a newtwork/spectrum analyzer on a passive filter. I am trying to derive the transer function, then develop the inverse filter. It is another story why I am trying to get the inverse filter. Matlab code is below, if it is not too much work for you too review it. I have a data file that I use with it. I have researched and found this paper, but the Matlab code would be involved to write for the IIR case: % grab_lab_network_analyzer_data.m % Author: Brian Farmer % Purpose: Grabs data, parses, and derives transfer function using IIR. % Clean up the environment clc clear all close all format compact [filename, pathname, filterindex] = uigetfile( ... '*.txt','Text file (*.txt)'); fprintf('Opening and analyzing the file "%s" ...', filename); rawfile = importdata(filename); fprintf('done.\n') % If first line is added in the data file: 0 -70 0 180 % Then the gain response drops a little more, but not flat % Amplitude is in dB (column 2), and phase is in degres (column 6) % Note: Degrees unwraps in the data file, but wraps in the graphical % plot. The graphical plot equals the data plot with a visual approximation % of (in this LPF test case): % -180 + -360 + -144 degrees is -684 degrees. -676 degrees final in file. DUT_frequency =,1); % column 1, frequency DUT_magnitude_dB =,2); % column 2, gain or amplitude (db) %column_3 =,3); column 3, toss 'NAN' text %column_4 =,4); column 4, toss 'NAN' text %DUT_frequency =,5);% column 5, frequency (same as column 1) DUT_phase =,6); % column 6, phase (degrees) % Convert data values in preparation for invfreqz % Don't know if this is power gain or voltage gain. % Power gain = 10*log10(voltage out/voltage in), this according to scope? % Voltage gain = 20*log10(voltage out/voltage in) for n = 1:length(DUT_magnitude_dB) DUT_magnitude(n) = 10^(DUT_magnitude_dB(n)/10); end for n = 1:length(DUT_phase) while DUT_phase(n) < -180 % loop to wrap DUT_phase(n) = DUT_phase(n) + 360; end DUT_rad(n) = DUT_phase(n)*(pi/180); x(n) = DUT_magnitude(n)*cos(DUT_rad(n)); y(n) = DUT_magnitude(n)*sin(DUT_rad(n)); end DUT_h = x + i*y; subplot(2,1,1) %semilogx(DUT_f,10*log10(abs(DUT_h))) %plot(DUT_f,10*log10(abs(DUT_h))) semilogx(DUT_frequency,10*log10(abs(DUT_h))) grid on ylabel('Magnitude (dB)') title('Bode Plot of Raw Lab Data') subplot(2,1,2) %semilogx(DUT_f,angle(DUT_h)*(180/pi)) %plot(DUT_f,angle(DUT_h)*(180/pi)) semilogx(DUT_frequency,angle(DUT_h)*(180/pi)) grid on ylabel('Phase (deg))') xlabel('Frequency (Rad/s)') pause DUT_h_inv = 1./DUT_h; subplot(2,1,1) semilogx(DUT_frequency,10*log10(abs(DUT_h_inv))) grid on ylabel('Magnitude (dB)') title('Bode Plot Inverse of Raw Lab Data') subplot(2,1,2) semilogx(DUT_frequency,angle(DUT_h_inv)*(180/pi)) grid on ylabel('Phase (deg))') xlabel('Frequency (Rad/s)') pause % Prepare for inverse % Remove any amplification that is beyound the amplitude resolution % of the oscilloscope's 8-bit ADC. magnitude_ADC = 20*log10(1/2^8); % 10dB or 20dB? % normalize to get the original waveform shape, not the orginal amplitude. DUT_magnitude_dB_norm = DUT_magnitude_dB - max(DUT_magnitude_dB); for n = 1:length(DUT_magnitude_dB) if DUT_magnitude_dB_norm(n) <= magnitude_ADC DUT_magnitude_ADC(n)= 1; % gain = 1 (no gain). reduces noise %DUT_magnitude_ADC(n)= DUT_magnitude_dB_min; % not same as passband max else DUT_magnitude_ADC(n) = 10^(DUT_magnitude_dB_norm(n)/10); end end for n = 1:length(DUT_phase) x_ADC(n) = DUT_magnitude_ADC(n)*cos(DUT_rad(n)); y_ADC(n) = DUT_magnitude_ADC(n)*sin(DUT_rad(n)); end DUT_h_ADC = x_ADC + i*y_ADC; % subplot(2,1,1) % semilogx(DUT_frequency,10*log10(abs(DUT_h_ADC))) % grid on % ylabel('Magnitude (dB)') % title('Bode Plot with ADC Amplitude Quantization Limit (48dB)') % subplot(2,1,2) % semilogx(DUT_frequency,angle(DUT_h_ADC)*(180/pi)) % grid on % ylabel('Phase (deg))') % xlabel('Frequency (Rad/s)') % % pause DUT_h_ADC_inv = 1./DUT_h_ADC; subplot(2,1,1) semilogx(DUT_frequency,10*log10(abs(DUT_h_ADC_inv))) grid on ylabel('Magnitude (dB)') title('Bode Plot Inverse Normalized and ADC Quantization Limited (48dB)') subplot(2,1,2) semilogx(DUT_frequency,angle(DUT_h_ADC_inv)*(180/pi)) grid on ylabel('Phase (deg))') xlabel('Frequency (Rad/s)') %break pause nyquist_frequency = DUT_frequency(end); sample_frequency = nyquist_frequency*2; DUT_f = DUT_frequency/nyquist_frequency; DUT_w = DUT_f*pi; %[bb,aa] = invfreqz(DUT_h,DUT_w,105,115); % for BPF %[bb,aa] = invfreqz(DUT_h,DUT_w,405,395); % for BPF [bb,aa] = invfreqz(DUT_h,DUT_w,404,426); % for BPF %[bb,aa] = invfreqz(DUT_h,DUT_w,41,51); % for LPF [FIT_h_eq_spaced,FIT_w_eq_spaced] = freqz(bb,aa,4096); % reverse (inverts) subplot(2,1,1) %plot((FIT_w_eq_spaced/pi)*nyquist_frequency,10*log10(abs(FIT_h_eq_spaced))) semilogx((FIT_w_eq_spaced/pi)*nyquist_frequency,10*log10(abs(FIT_h_eq_spaced))) grid on ylabel('Magnitude (dB)') title('Bode Plot') subplot(2,1,2) %plot((FIT_w_eq_spaced/pi)*nyquist_frequency,angle(FIT_h_eq_spaced)*(180/pi)) semilogx((FIT_w_eq_spaced/pi)*nyquist_frequency,angle(FIT_h_eq_spaced)*(180/pi)) grid on ylabel('Phase (deg))') xlabel('Frequency (Hz)') %instead of Rad/s pause %break %[FIT_h,FIT_f] = freqz(bb,aa,DUT_frequency,sample_frequency); %FIT_f and DUT_frequency are the same % FIT_h = freqz(bb,aa,DUT_w); % % subplot(2,1,1) % semilogx(FIT_f,10*log10(abs(FIT_h))) % grid on % ylabel('Magnitude (dB)') % title('Bode Plot') % subplot(2,1,2) % semilogx(FIT_f,angle(FIT_h)*(180/pi)) % grid on % ylabel('Phase (deg))') % xlabel('Frequency (Hz)') xx = zeros(2048); xx(64) = 1; yy = filter(bb,aa,xx); clf impz(bb,aa); %plot(abs(yy),'o') break pause [bbb,aaa] = invfreqz(DUT_h_ADC_inv,DUT_w,426,404); [h,w] = freqz(bbb,aaa,1024); %inverse filter just reverse (inverts) %y_inv = filter(aa,bb,yy); y_inv = filter(bbb,aaa,yy); plot(y_inv) break %w = 1350000*w; % to convert to Hz plot(w/pi,10*log10(abs(h))) pause plot(w/pi,angle(h)*180/pi) title('Frequency Response Magnitude') break subplot(2,1,1) %semilogx(w/pi,10*log10(abs(h))) plot(w/pi,10*log10(abs(h))) grid on ylabel('Magnitude (dB)') title('Bode Plot') subplot(2,1,2) %semilogx(w/pi,angle(h)*(180/pi)) plot(w/pi,angle(h)*(180/pi)) grid on ylabel('Phase (deg))') xlabel('Frequency (Rad/s)') break diary test.txt %or mat file? disp('data') disp(' ') for n = 1:1024 disp(['10#' int2str(DUT_amplitude(n))]) end diary off
8 years ago
Sorry, you need javascript enabled to post any comments.
Sorry I just now saw your comment (hopefully I will catch comment notifications from now on).

It is generally true that a growing impulse response will yield an unstable filter. What I normally do is convert to minimum phase before using invfreqz. That usually takes care of it.
4 years ago
Sorry, you need javascript enabled to post any comments.
Here is my implementation of db() in case anyone needs it:

function [xdb] = db(x,mindb)
% [y] = db(x,mindb)
% convert x to db, optionally clipped to mindb - see also dbn()

ax = abs(x);
xmax = max(ax(;
if nargin<2, xmin = xmax*eps; else xmin = 10^(mindb/20); end
xdb = 20*log10(max(abs(x),xmin*ones(size(x))));
3 years ago
Sorry, you need javascript enabled to post any comments.
I'm wondering if there's any advice here! I'm designing a filter to model the frequency response of an old tape. I have that impulse response data, H(z), as a '.wav' file, and am trying to design a filter from it. From what I can tell, while the program runs, the issue or uncertainty I'm having is what order of zeros (NB) and poles (NA) to use for invfreqz()... [B,A]=invfreqz(H,W,NB,NA,Wt,ITER,TOL). Although you don't have the audio file I'm using, the code it as follows:

close all

[y, fs] = audioread('Metal.wav');

y = y - mean(y);
y = detrend(y);
y1 = abs(fft(y));

[p,f] = pspectrum(y,fs,'FrequencyResolution',1000); %*((2*pi)/fs)

grid on
title('Noise Power Spectrum - Metal');
xlabel('Frequency (Hz)')
ylabel('Power Spectrum (dB)')
set(gca, 'XScale', 'log')

[b,a] = invfreqz(p(:,1),f,20,22, [], 10000, 0.00000001);

1 year ago
Sorry, you need javascript enabled to post any comments.
Sorry, you need javascript enabled to post any comments.