freqz() magnitude plot behavior for estimated filter coefficients

I am working on a project involving estimating a digital filter representation of an experimentally-determined frequency response.

The system in question has a frequency response that, interestingly, has a slight increase in magnitude as frequency increases, like the following plot:

I've estimated FIR filter coefficients using the invfreqz() function with denominator n = 0, and numerator m = 200. This results in a set of coefficients, a and b, respectively.

Given a set of frequencies used in the initial experiment in which we obtained the frequency response, stored and normalized in the vector w -> [0,pi], validating these coefficients using the function freqz() in the following manner results in the following plot. It appears to be pretty consistent with the initial frequency response data.

h = freqz(b,a,w);
subplot(2,1,1); plot(w,abs(h)); grid on; ylabel('Magnitude'); xlabel('Frequency (Hz)');
subplot(2,1,2); plot(w,angle(h)); grid on; ylabel('Phase (rad)'); xlabel('Frequency (Hz)');

However, when I don't input the frequency vector, and instead let freqz() output the frequency bins, the resulting plot is quite different:

[h,w] = freqz(b,a);
subplot(2,1,1); plot(w,abs(h)); grid on; ylabel('Magnitude'); xlabel('Frequency (Hz)');
subplot(2,1,2); plot(w,angle(h)); grid on; ylabel('Phase (rad)'); xlabel('Frequency (Hz)');

The frequency vector from the experiment consists of only 101 points between 60 and 10060 Hz, so I suspect the freqz() function is interpolating to fill out the default size-512 frequency output vector in the latter case. But I'm not sure why that results in such a strong sinusoidal pattern, or why the phase appears to tend around 0 radians for the entire spectrum, rather than a slight linear decrease.

What could be going on here?

[Paul]

Hi Aaron,

The text states that "the vector w -> [0,pi]"

Then the code shown is:

h = freqz(b,a,w);

That call to freqz is fine. As long as w has two or more elements it is intrepreted as angular frequency in units of rad/sample.

Then we have the plotting commands:

subplot(2,1,1); plot(w,abs(h)); grid on; ylabel('Magnitude'); xlabel('Frequency (Hz)');
subplot(2,1,2); plot(w,angle(h)); grid on; ylabel('Phase (rad)'); xlabel('Frequency (Hz)');

Here, w is the independent variable but is xlabeled as Hz. But if it's the same w that was input to freqz it should be in rad/sample. However, and more importantly, w should only range from 0 to pi. But the plots are going from 0 to 10000.

Then, in this code

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

w is returned as angular frequency (rad/sample) ranging from 0 to pi, but the plotting commands are again labeled in Hz and the plots show that w ranges from 0 to 10000.

subplot(2,1,1); plot(w,abs(h)); grid on; ylabel('Magnitude'); xlabel('Frequency (Hz)');
subplot(2,1,2); plot(w,angle(h)); grid on; ylabel('Phase (rad)'); xlabel('Frequency (Hz)');

So there is a disconnect somewhere in this question between what's being coded and/or stated and what is being shown.

I suggest you save your a, b, and w vectors (from the first call to freqz) into a .mat file and then upload the .mat file using the paperclip icon on the Insert menu, either as an edit to the original question or in a response to this comment.

It might also be helpful to show all of your code. Not necessarily the code used to develop b and a, but everything after that to develop the plots.

[Aaron Wilson]

Thanks for the comment. You are right, I put the incorrect x-axis vector in the plotting function here on MathWorks, but when generating the plots I simply converted the w vector from rad/s to Hz.

[Paul]

Except I thought the w vector was in rad from 0->pi, not rad/s?

In any case, I suggest uploading your b, a, and w vectors in .mat file using the paper clip icon on the Insert menu.

Answer 1

[Star Strider]

The third argument to freqz (when the first two are the transfer function coefficients, different with a second-order-section matrix or a digital filter object) is the number of points at which to evalate the transfer function (in this instance). Ideally, that should be a scalar.

Experimenting with that idea —

b = fir1(200, 0.9, 'high');                                 % Prototype filter
a = 1;

w = linspace(60, 1E+4, 101)/(2*pi)                          % Radian Frequency Vector (Guess)

w = 1×101

9.5493 25.3693 41.1893 57.0093 72.8293 88.6493 104.4693 120.2893 136.1093 151.9293 167.7493 183.5693 199.3893 215.2093 231.0293 246.8493 262.6693 278.4893 294.3093 310.1293 325.9493 341.7693 357.5893 373.4093 389.2293 405.0493 420.8693 436.6893 452.5093 468.3293

figure
h = freqz(b,a,w);
subplot(2,1,1); plot(w,abs(h)); grid on; ylabel('Magnitude'); xlabel('Frequency (Hz)');
subplot(2,1,2); plot(w,angle(h)); grid on; ylabel('Phase (rad)'); xlabel('Frequency (Hz)');

figure
[h,w] = freqz(b,a);
subplot(2,1,1); plot(w,abs(h)); grid on; ylabel('Magnitude'); xlabel('Frequency (Hz)');
subplot(2,1,2); plot(w,angle(h)); grid on; ylabel('Phase (rad)'); xlabel('Frequency (Hz)');

figure
freqz(b, 1, 2^16, 10000)                                    % High-Resolution Version

What you are seeing may be a very low resolution version of the actual transfer function in the first instance, although I have no idea how freqz interprets a vector input for the third argument, since it expects a scalar, and apparently doesn’t check to be sure that it is, or that the value is an integer. When you leave freqz to its own devices, it produces a higher resolution version.

Without your code (or your ‘b’ vector), this is only a guess on my part, however it could explain the discrepancy.

The problem is likely not with freqz, and is instead what you are giving it.

.