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

Question

Abhinit 2025-8-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2179494-create-a-code-for-butterworth-4th-order-bandpass-filter-without-signal-processing-toolbox

评论： William Rose 2025-8-23

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2025-8-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2179494-create-a-code-for-butterworth-4th-order-bandpass-filter-without-signal-processing-toolbox#answer_1569690

在 MATLAB Online 中打开

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.

.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

William Rose 2025-8-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2179494-create-a-code-for-butterworth-4th-order-bandpass-filter-without-signal-processing-toolbox#answer_1569687

编辑：Torsten 2025-8-23

在 MATLAB Online 中打开

@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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

William Rose 2025-8-23

@Torsten, thank you.

William Rose 2025-8-23

在 MATLAB Online 中打开

emg+force.txt

@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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论