Reputation: 31
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.
Upvotes: 3
Views: 1074
Reputation: 3177
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);
Upvotes: 1
Reputation: 14579
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 b
s are the coefficients of the x
sequence and a
s are the coefficients of the y
sequence, such that the b
are in the numerator and the a
s 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');
Upvotes: 0
Reputation: 1331
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.
Upvotes: 0