Return to 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 stable signal:

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

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):

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 stable signal:

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 stable signal:

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

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; % Pass band frequencies
Ws = [20, 90]/fs; % 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 stable signal:

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 stable signal: