How to generate frequency response given b,a coefficients of the system?

Question

I have the following system, specified by the set of coefficients:

b = [1 2 3];
a = [1 .5 .25];

In the Z-Domain, such function will have the following transfer function:

 H(Z) = Y(Z)/X(Z)

So the frequency response will be just the unit circle, where:

 H(e^jw) = Y(e^jw)/X(e^jw)

Do I just substitute in the e^jw for 'Z' in my transfer function to obtain the frequency response of the system mathematically, on paper? Seems a bit ridiculous from my (a student's) point of view.

Answer 1

Have you tried freqz()? It returns the frequency response vector, h, and the corresponding angular frequency vector, w, for the digital filter with numerator and denominator polynomial coefficients stored in b and a, respectively.

In your case, simply follow the help:

[h,w]=freqz(b,a);

Answer 2

You do sub in e^jw for Z. This isn't ridiculous. Then you just sweep w from -pi to pi. Your freq response will be the absolute value of the result.

As Alessiox mentioned, freqz is the command you want to use in matlab.

Answer 3

I would indeed be as simple as substituting exp(j*w) in your transfer function. There are of course different ways to implement this with Matlab. For the purpose of illustration, I will be assuming bs are the coefficients of the x sequence and as are the coefficients of the y sequence, such that the b are in the numerator and the as are in the denominator:

A direct evaluation with Matlab could be done with:

b = [1 2 3];
a = [1 .5 .25];

N = 513; % number of points at which to evaluate the transfer function
w = linspace(0,2*pi,N);
num = 0;
for i=1:length(b)
  num = num + b(i) * exp(-j*i*w);
end
den = 0;
for i=1:length(a)
  den = den + a(i) * exp(-j*i*w);
end
H = num ./ den;

This would be equivalent to the following which makes use of the builtin polyval:

N = 513; % number of points at which to evaluate the transfer function
w = linspace(0,2*pi,N);
H = polyval(fliplr(b),exp(-j*w))./polyval(fliplr(a),exp(-j*w));

Also, this is really evaluating the transfer function at discrete equally spaced angular frequencies w = 2*pi*k/N which corresponds to the Discrete Fourier Transform (DFT). As such it could also be done with:

N = 512;
H = fft(b,N) ./ fft(a,N);

Incidentally this is what freqz does, so you could also get the same result with:

N = 512;
H = freqz(b,a,N,'whole');