bandpass - The reasons for filter explosion

Assume signal $x$, sampled at $f_s = 44100\; \mathrm{Hz}$. I tried to filter it using the Butterworth bandpass filter ( $30\; \mathrm{Hz} - 70\; \mathrm{Hz}$) of order $8$. However, as a result I get a vector with most elements being NaN (and some of them extremely small, approx. $-2.5 \cdot 10^{306}$`).

If I try the same filter of order $6$, I get results as expected. What could be possible reason for order $8$ filter to 'explode'?

Here is the MATLAB code, just in case I made an error which I don't see:

[b, a] = butter(4, [60 / fs, 140/fs]);

x_filtered = filter(b, a, x);

Answer

Basically you never want to use the Transfer Function representation (with b and a) and rather use the Zeros-Poles-Gain (z,p,k). This will allow you to avoid the numerical errors. In your case you might design your filter in following way:

fs = 44100;           % Sampling frequency
Wp = [30, 70]/(fs/2); % Pass band frequencies (as normalized frequency)
Ws = [20, 90]/(fs/2); % Stop band frequencies
Rp = 3;               % Ripple at pass band
Rs = 50;              % Ripple at stop band

[n, Wn] = buttord(Wp, Ws, Rp, Rs);     % Get order and omega vector

[z, p, k] = butter(n, Wn, 'bandpass'); % Design filter accordingly
[sos, g] = zp2sos(z, p, k);            % Convert to state matrix
Hd = dfilt.df2sos(sos, g);             % Create the filter object

Which for some dummy random signal:

x = rand(1,100000);
y = filter(Hd, x);

Will produce a stable output:

And here is the filter frequency response (everything looks as requested):