Sign in

Not a member? | Forgot your Password?

Search code

Search tips

Free PDF Downloads

A Quadrature Signals Tutorial: Complex, But Not Complicated

Understanding the 'Phasing Method' of Single Sideband Demodulation

Complex Digital Signal Processing in Telecommunications

Introduction to Sound Processing

C++ Tutorial

Introduction of C Programming for DSP Applications

Fixed-Point Arithmetic: An Introduction

Cascaded Integrator-Comb (CIC) Filter Introduction

FIR Filter Design Software

See Also

Embedded SystemsFPGA

DSP Code Sharing > Fitting Filters to Measured Amplitude Response Data Using invfreqz in Matlab

Fitting Filters to Measured Amplitude Response Data Using invfreqz in Matlab

Language: Matlab

Processor: Not Relevant

Submitted by Julius Orion Smith III on Nov 8 2010

Licensed under a Creative Commons Attribution 3.0 Unported License

Fitting Filters to Measured Amplitude Response Data Using invfreqz 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

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

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

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

 
 
Rate this code snippet:
4.75
Rating: 4.75 | Votes: 4
 
   
 
posted by Julius Orion Smith III
Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and Associate Professor (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.


Comments


 

cfelton wrote:

11/9/2010
 
Good description and comments!
 

Rick Lyons wrote:

12/10/2010
 
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-]
 

JOS wrote:

12/10/2010
 
Fixed - thanks! - jos
 

brianwfarmer wrote:

7/18/2011
 
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:

http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1084468



% 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 = rawfile.data(:,1); % column 1, frequency
DUT_magnitude_dB = rawfile.data(:,2); % column 2, gain or amplitude (db)
%column_3 = rawfile.data(:,3);       column 3, toss 'NAN' text
%column_4 = rawfile.data(:,4);       column 4, toss 'NAN' text
%DUT_frequency = rawfile.data(:,5);% column 5, frequency (same as column 1)
DUT_phase     = rawfile.data(:,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

Add a Comment
You need to login before you can post a comment (best way to prevent spam). ( Not a member? )