Frequency response of a continuous-time (analog) filter

Summary

A matlab example that shows how to calculate the frequency response of a continuous-time (s-domain, “analog”) filter.

Introduction

Matlab has a built-in function (freqs) to calculate and plot the frequency response of a continuous-time filter.
This example shows how it is done step-by-step.

S-domain filter

The frequency response of a continuous-time filter has the following structure:

$H(s)=\frac{b_{m}s^{m}+b_{m-1}s^{m-1}+\ldots+b_{1}s+b_{0}}{a_{n}s^{n}+a_{n-1}s^{n-1}+\ldots+a_{1}s+a_{0}}$

$s=\sigma+i\omega$ is the argument of the Laplace transform, which simplifies with $\sigma=0$ to the continous Fourier transform.
The frequency response is found by evaluating

at a given angular frequency $\omega$ with $s=i\omega$ .

The resulting $H(i\omega)$ gives the frequency response as a complex number for each frequency.
Its amplitude can be converted to dB using $20\log_{10}|H(i\omega)|$ .

Example program

The program generates an example filter, calculates the frequency response on a logarithmic axis, and plots it (figure 1, 2).

Figure 1: Amplitude response

Figure 2: Phase response

Download

% ***********************************************************************
% Continuous-time ("analog") filter example
% © Markus Nentwig 2007
% This program is provided without any warranty.
% ***********************************************************************
close all; clear all;

% Filter design parameters
order=5;
fCutoff=2000;

% A fifth-order polynomial of frequency needs to be evaluated.
% Unless frequencies are normalized,  f^1 and f^5 are of very different
% order of magnitude, and numerical accuracy would be insufficient. 
% Therefore, divide all frequencies by fCutoff.

fNorm=2.3456789*fCutoff; % (almost) arbitrary choice!
fCutoff=fCutoff/fNorm;

% ***********************************************************************
% Design a filter
% ***********************************************************************
if 0 
  [b,a]=butter(order, 2*pi*fCutoff, 's');
elseif 1
  ripple_dB=0.5;
  [b,a]=cheby1(order, ripple_dB, 2*pi*fCutoff, 's');
elseif 0
  ripple_dB=0.5;
  stopband_dB=60;
  [b,a]=ellip(order, ripple_dB, stopband_dB, 2*pi*fCutoff, 's');
end

% ***********************************************************************
% Build a logarithmic frequency axis
% ***********************************************************************
f1=10;
f2=48000;
nPoints=3000;
fbase=exp(linspace(log(f1), log(f2), nPoints));
% The filter requires normalized frequency
fbase=fbase/fNorm;

% ***********************************************************************
% Numerator and denominator are polynomials in s=i*omega+alpha (see 
% Laplace transform). Alpha is 0.
% ***********************************************************************
H_num=polyval(b, i*2*pi*fbase);
H_denom=polyval(a, i*2*pi*fbase);
% H is the complex gain of the filter at any given frequency
H=H_num ./ H_denom;

% ***********************************************************************
% - Convert to dB. H could be a voltage or current => use 20 log |H|
% - Obtain phase in degrees
% - plot
% ***********************************************************************
gain_dB=20*log10(abs(H)+1e-15); % don't show below -300 dB
phase_deg=unwrap(arg(H))/pi*180;
figure(); semilogx(fbase*fNorm, gain_dB, '+'); xlabel('f/Hz'); ylabel('|H|/dB'); grid on;
figure(); semilogx(fbase*fNorm, phase_deg, '+'); xlabel('f/Hz'); ylabel('phase/deg'); grid on;

© Markus Nentwig 2007-2008
The content of this page is provided without any warranty and may not be reproduced without permission.

Comments? Questions?

Please send me a mail! mnentwig@elisanet.fi