create a code for butterworth 4th order bandpass filter without Signal Processing Toolbox

I am working with Surface EMG data and need to use a 4th order butterworth bandpass filter but don't have the signal processing toolbox. The cutoff frequency range is from 20 to 450 hz. I am also trying to understand how to come up with the coefficients using this bandpass filter and where I can apply that if I'm using that filtered signal to measure the MVIC of the data.

Thank you

Answer 1

[Star Strider]

A relatively straightforward presentation stats with (for a fourth-order prototype lowpass filter) --

where is a normalising constant and --

and then do a lowpass-to bandpass transformation --

where and are the high and low corner frequencies, resspectively.

Coding it --

syms K_0 s t 
lpf =  20;
hpf = 450;
omega_L = 2*pi*lpf

omega_L = 125.6637

omega_H = 2*pi*hpf

omega_H = 2.8274e+03

sbp = @(s) (s^2 + omega_H/omega_L) / (s*(omega_H - omega_L))

sbp = function_handle with value:
    @(s)(s^2+omega_H/omega_L)/(s*(omega_H-omega_L))

n = 4;
sexp = @(k,n) exp((1j*pi/2)*(1+(2*k*1)/n))

sexp = function_handle with value:
    @(k,n)exp((1j*pi/2)*(1+(2*k*1)/n))

for k = 1:n
    st(k) = s - sexp(k,n);
end
H = vpa(expand(1/prod(st)), 7)

H = 

[n,d] = numden(H)

n = 

1.0e+48

d = 

H = vpa(simplify(expand(subs(H,s,sbp(s))), 500), 7)

H = 

figure
fplot(abs(H), [1, 500])
grid
set(gca, XScale='log', YScale='log')        % The 'mag2db' function doesn't work here for some reason

h = vpa(simplify(ilaplace(H, s, t), 500), 4)

h = 

●

Using the Signal Processing Toolbox --

[b,a] = butter(4, [lpf, hpf]/500, 'bandpass');

figure
freqs(b, a, 2^16)

I doubt that you will get the bandpass result you want with a 4th order Butterworth design.

To use this with a sampled signal, you would then need to use the bilinear transform.

I would use an elliptic filter, such as that designed by the bandpass function with the 'iir' option --

Fs = 1000;
L = 1;
t = linspace(0, Fs*L, Fs*L+1)/Fs;
s = randn(size(t));
[sf,df] = bandpass(s, [lpf, hpf], Fs, ImpulseResponse='iir');
df

df = 
 digitalFilter with properties:

   Coefficients:
               Numerator: [8×3 double]
             Denominator: [8×3 double]

   Specifications:
       FrequencyResponse: 'bandpass'
         ImpulseResponse: 'iir'
              SampleRate: 1000
      PassbandFrequency2: 450
      StopbandFrequency2: 457.8500
    StopbandAttenuation1: 60
    StopbandAttenuation2: 60
          PassbandRipple: 0.1000
      StopbandFrequency1: 16.8600
      PassbandFrequency1: 20
            DesignMethod: 'ellip'

 Use filterAnalyzer to visualize filter
 Use designfilt to edit filter
 Use filter to filter data

figure
freqz(df, 2^16, Fs)

% sgtitle('Efficient Elliptic Filter Designed By the ''bandpass'' Function')

If you are doing biomedical signal processing, you really need to get the Signal Processing Toolbox.

.

Answer 2

[William Rose]

@Abhinit,

[Edit:

Make corections to formulas for a1 and a2 in the lowpass filter.

Many of the equations I entrered with the LaTeX editor look OK when I enter them, but after I submit my answer, they do not display correctly on my monitor. I hope this is not happening to you also. It happens to me with Firefox and with Edge browsers.]

Since you want a 4th order Butterworth bandpass, make two second order Butterworths: low pass at 450 Hz plus high pass at 20 Hz.

In theory, the order (low then high, versus high then low) does not matter. Whether this is true in practice depends on the data and the choice of filters and cutoff frequencies.

x=input, y1=output from filter 1=input to filter 2, y2=output from filter 2

We have two second order digital filters with transfer functions

where the coefficients ak and bk are different for the two filters.

Let fs=sampling rate and T=1/fs=sampling interval.

Lowpass flter

Define fcL=lowpass cutoff frequency (Hz); =lowpass cutoff frequency (radian/s);

.

betaL = tan(wcL*T/2) = tan(pi*fcL/fs) (same equation as above but without LaTeX)

Then, for the 2nd order Butterworth lowpass filter,

, ,

,

Define bL=[b0,b1,b2] and aL=[1,a1,a2] for the lowpass filter with the formulas above. Then apply the lowpass filter to signal x by doing

y1=filter(bL,aL,x);

Function filter() is part of standard Matlab - not in a toolbox.

Highpass filter

Define fcH=highpass cutoff frequency (Hz); =highpass cutoff frequency (radian/s);

Then, for the 2nd order Butterworth highpass filter,

, ,

,

Define bH(1:3)=[b0,b1,b2] and aH(1:3)=[1,a1,a2] with the highpass formulas above. Then apply the highpass filter to signal y1 by doing

y2=filter(bH,aH,y1);

y2 is the bandpass-filtered signal that you want.

Good luck with your analysis.

[Torsten]

@William Rose

Opening and submitting your answer made the formulae look correct.

[William Rose]

@Torsten, thank you.

[William Rose]

@Abhinit,

Here is an example. The data file has EMGs from four leg muscles (columns 1-4), and ground reaction force (column 5, Newtons), all recorded at 1 kHz. The subject does a vertical jump on a forceplate. Column 1= vastus lateralis muscle EMG.

data=load('emg+force.txt');
x=data(:,1);        % vastus lateralis EMG
fs=1000;            % sampling rate (Hz)

Compute coefficients for digital 2nd order lowpass Butterworth filter at 450 Hz, using formulas from my answer above:

fcL=450;            % lowpass cutoff frequency (Hz)
Wn=2*fcL/fs;        % normalized frequency 
beta=tan(Wn*pi/2);
b0=beta^2/(1+sqrt(2)*beta+beta^2);
b1=2*b0;
b2=b0;
a1=2*(beta^2-1)/(1+sqrt(2)*beta+beta^2);
a2=(1-sqrt(2)*beta+beta^2)/(1+sqrt(2)*beta+beta^2);
b=[b0,b1,b2];       % numerator coefficients for H(z)
a=[ 1,a1,a2];       % denominator coefficients for H(z)

Apply the lowpass filter to the signal:

y1=filter(b,a,x);

Compute coefficients for digital 2nd order highpass Butterworth filter at 20 Hz:

fcH=20;             % highpass cutoff frequency (Hz)
Wn=2*fcH/fs;        % normalized frequency
beta=tan(Wn*pi/2);
b0=1/(1+sqrt(2)*beta+beta^2);
b1=-2*b0;
b2=b0;
a1=2*(beta^2-1)/(1+sqrt(2)*beta+beta^2);
a2=(1-sqrt(2)*beta+beta^2)/(1+sqrt(2)*beta+beta^2);
b=[b0,b1,b2];       % numerator coefficients for H(z)
a=[ 1,a1,a2];       % denominator coefficients for H(z)

Apply the highpass filter to the signal that has been lowpass-filtered:

y2=filter(b,a,y1);

Display original signal and filtered signal.

t=(0:length(x)-1)/fs;       % time vector
plot(t,x,'-b',t,y2,'-r');
grid on; xlabel('Time (s)'); ylabel('EMG');
title('Vastus Lateralis Muscle EMG During A Jump');
legend('Unfiltered','Filtered');

●

The formulas above produce the same results as if you use the toolbox functions butter(2,Wn) and butter(2,Wn,'high') to get the coefficients for the lowpass and highpass filters.