Why does the phase response get cutoff beyond a frequency?

Hi, I am designing some butterworth lowpass filters of different orders and cutoff frequencies to test on some ECG data sampled at 1000 Hz. I am designing them as:

1. fc = 10Hz; order = 2
2. fc = 20Hz; order = 8
3. fc = 40Hz; order = 8
4. fc=70Hz; order = 8

I am using the same structure for all the filter designs, with just a change in variable names, like so:

fs = 1000;
fc1 = 10;
N1 = 2;
[b1,a1] = butter(N1,fc1/(fs/2),'low'); %low indicates LPF
figure('Name','LPF 20 Hz');
freqz(b1,a1,L,fs);

However I get a difference in the phase response for the 10Hz order 2 filter and the 20Hz order 8 filter. Here is the phase response of the 10 Hz filter:

and the phase response of the 20 Hz filter:

Can someone please explain why the second phase response is getting cut off abruptly?

Answer 1

[Star Strider]

The most likely reason is that the filter numerator coefficients are close to the limit of floating-point approximation error. This is a problem with transfer-function realisations of digital filters in general, and the reason that the second-order-section realisation is preferred.

I did the same filter with a second-order-section realisation (and initiial zero-pole-gain output from butter), and it works as desired —

L = 2^16;
fs = 1000;
fc1 = 20;
N1 = 8;
[b1,a1] = butter(N1,fc1/(fs/2),'low'); %low indicates LPF
format longG
ND = [b1; a1]

ND = 2×9

1.77837938807345e-10 1.42270351045876e-09 4.97946228660565e-09 9.95892457321131e-09 1.24486557165141e-08 9.95892457321131e-09 4.97946228660565e-09 1.42270351045876e-09 1.77837938807345e-10 1 -7.35591423914845 23.6970393591454 -43.6658296202337 50.3366193367395 -37.1713473565871 17.171275470867 -4.53667256190413 0.524829656648063

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

format short
figure('Name','LPF 20 Hz');
freqz(b1,a1,L,fs);
h1 = subplot(2,1,1);
h1.YLim = [min(h1.Children.YData) 50];

[z,p,k] = butter(N1,fc1/(fs/2),'low'); %low indicates LPF
[sos,g] = zp2sos(z,p,k);
SOS = sos

SOS = 4×6

1.0000 2.0000 1.0000 1.0000 -1.7670 0.7811 1.0000 2.0000 1.0000 1.0000 -1.7970 0.8112 1.0000 2.0000 1.0000 1.0000 -1.8551 0.8698 1.0000 2.0000 1.0000 1.0000 -1.9369 0.9523

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

G = g

G = 1.7784e-10

figure('Name','LPF 20 Hz');
freqz(sos,L,fs);
h1 = subplot(2,1,1);
% get(h1)
h1.Children.YData = h1.Children.YData - h1.Children.YData(1);
h1.YLim = [min(h1.Children.YData) 50];

.

[Shruti Jeyaraman]

Thank you so much, I will keep this in mind as I do filter designs in the future. i wonder what's causing the FP error though. Nevertheless, I tried the SOS implementation and it worked out great.

[Star Strider]

As always, my pleasure!

Thee floatiing-point error is siimply due to the values of the coefficients (on the order of ) so apparently beyond about 280 Hz the real and iimaginary parts of the Fourier transforms of the filter output are both 0 or close to it, and since the atan2 function that is used in the angle calculation requires the imag/real division to be performed, that produces a 0/0 result. A 0/0 result equates to NaN, and NaN values do not plot. (This iis my assumption. I did not actually explore that by calculating it.)

Answer 2

[William Rose]

Check the values returned by freqz. Answering on my cellphone so I’m limited in what I can do. Look for NaNs at high frequencies in the computed 8th order response. The [b,a]=butter() method of specifying a filter can be unstable even at relatively low order. Try z,p,k or second order sections instead.

[Shruti Jeyaraman]

Thanks a lot for responding, i will analyze what freqz returns as you suggest. As the other answer suggested, I implemented the design with an SOS and that worked out