1. Experimental purpose
1. Be familiar with the principles and methods of designing IIR digital filters using bilinear transformation method.
2. Master the computer simulation method of digital filters.
3. Obtain perceptual knowledge of digital filtering by observing the filtering effect on actual electrocardiogram signals.
2. Experimental content
1. Use bilinear transformation method to design a Butterworth low-pass IIR digital filter. The design index parameters are: when the frequency in the passband is lower than 0.2Π, the maximum attenuation is less than 1dB; in the stopband [0.3Π, Π] frequency interval, the minimum attenuation is greater than 15dB.
2. With T=1s as the sampling interval, print out the amplitude-frequency response characteristic curve of the digital filter in the frequency interval [0, Π/2].
3. Use the designed filter to perform simulated filtering on the actual ECG signal sampling sequence (given later in this experiment), and print out the ECG signal waveforms before and after filtering to observe and summarize the filtering effects and effects.
3. Experimental steps
1. Review the content about Butterworth analog filter design and the use of bilinear transformation method to design IIR digital filters, and use the bilinear transformation method to design the digital filter system function H(z).
The digital filter function in the textbook that meets the requirements of this experiment 
In the formula: 
,k=1,2,3
A=0.09036,
B1=1.2686, C1=-0.7051
B2=1.0106, C2=-0.3583
B3=0.9044, C3=-0.2155
According to the design indicators, H(z) can also be obtained by calling the MATLAB signal processing tool functions buttonrd and butter.
It can be seen from the function of the filter that the filter H(z)
It consists of three second-order filters H1(z), H2(z), H3(z) cascaded.
2. Write a filter simulation program to calculate H(z) for the ECG signal sampling sequence 
response sequence y(n).
yk(n) is the output sequence of the Kth-order filter Hk(z), Yk-1(n) is the input sequence, and the difference equation can be obtained according to the composition of the filter H(z).
When K =1 when
So the total response sequence y(n) of H(z) to x(n) can be obtained using a sequential iteration algorithm. That is, the difference equation is solved for k=1, 2, and 3 in sequence, and finally y3(n)=y(n) is obtained. The simulation program is a general program that implements the above-mentioned solution of the difference equation and the sequential iteration algorithm. The MATLAB filter function can also be directly called to implement the simulation. .
3. Run the simulation filter program on a general-purpose computer to complete the experimental contents (2) and (3).
Human ECG signals are often interfered by industrial high-frequency interference during the measurement process, so they must be processed by low-pass filtering before they can be used as useful information for judging heart function. An actual electrocardiogram signal sampling sequence sample x(n) is given below, in which high-frequency interference exists. In the experiment, x(n) is used as the input sequence to filter out the interference components. x(n) ={-4, -2, 0, -4, -6, -4, -2, -4, -6, -6, – 4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, – 4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0}
The code looks like this:
%%
clear all;
close all;
T=1; % can be chosen arbitrarily, usually 1
Fs=1/T;
Wp=0.2*pi;
Ws=0.3*pi;
Rp=1;
As=15;
% gives the digital filter index, which is first converted into the prototype analog filter index.
Wap=2/T*tan(Wp/2);
Was=2/T*tan(Ws/2);
[N,Wn]=buttord(Wap,Was,Rp,As,'s'); % order estimation function (in line with the order and cutoff frequency of the index)
[z0,p0,k0]=buttap(N); %Normalized Butterworth low-pass simulation prototype function (z0, p0, k0 are zero point, pole and gain respectively)
[Bap,Aap]=zp2tf(z0,p0,k0); %Convert the zero and pole points of the system function into the coefficients of the general form of the system function, and convert the zero points and poles into the transfer function representation of the simulated low-pass filter
[bs,as]=lp2lp(Bap,Aap,Wn); % frequency conversion function to convert the analog normalized prototype low-pass filter into the required de-normalized analog low-pass filter
[bz,az]=bilinear(bs,as,Fs); % Bilinear transformation method for analog to digital filter transformation
%Print out the numerator and denominator polynomial coefficients of H(z)
bz;
az;
%With 0.02pi as the sampling interval, the amplitude-frequency response characteristic curve in the frequency interval [0,pi/2]
[H,W]=freqz(bz,az); % Calculate frequency response freqz digital filter frequency response
figure(1); %Draw the amplitude-frequency response curve
plot(W,abs(H));
figure(2);
title('Logarithmic amplitude-frequency characteristic curve graph');
plot(W,20*log10(abs(H)));
%signal filtering
xn=[-4, -2, 0, -4, -6, -4, -2, -4, -6, -6,-4, -4, -6, -6, -2, 6, 12 , 8, 0, -16,-38,-60, -84, -90, -66, -32, -4, -2, -4, 8,12, 12, 10, 6, 6, 6, 4 , 0, 0, 0,0, 0, -2, -4, 0, 0, 0, -2, -2, 0,0, -2, -2, -2, -2, 0];
figure(3);
subplot(1,1,1);
stem(xn);
y=filter(bz,az,xn);% performs filtering, x is which signal is filtered, y is the filter output result
figure(4);
subplot(1,1,1);
stem(y);% filtering produced by y
%%
clear all;
close all;
k=1;
x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12 ,8, 0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4 ,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0,0];
figure(1)
subplot(2,1,1)
n=0:56;
stem(n,x,'o');% draw discrete sequence data
axis([0 56 -100 80]);
hold on;
n=0:60;
m=zeros(61);%Create a 61×61 matrix composed of zero values.
ylabel('x(n)');
title('ECG signal sampling sequence x(n)');
B=[0.09036 2*0.09036 0.09036];
A=[1.2686 -0.7051];
plot(n,m);
xlabel('n');
A1=[1.0106 -0.3583];
A2=[0.9044 -0.2155];
while(k<=3)
y=filter(B,A,x);
x=y;
if k==2
A=A1;
end
if k==3
A=A2;
end
k=k + 1;
end
subplot(2,1,2)
n=0:56;
stem(n,y,'*');
axis([0 56 -15 15]);
hold on;
n=0:60;
m=zeros(61);
plot(n,m);
xlabel('n');
ylabel('y(n)');
title('Bilinear transformation method to design IIR filtered electrocardiogram signal');
%Bilinear transformation method to design ButterWorth digital filter amplitude-frequency characteristics
A=[0.09036,0.1872,0.09036];
B1=[1,-1.2686,0.7051];
B2=[1,-1.0106,0.3583];
B3=[1,-0.9044,0.2155];
[H1,~]=freqz(A,B1,100);% digital filter frequency response
[H2,~]=freqz(A,B2,100);
[H3,w]=freqz(A,B3,100);
H4=H1.*(H2);
H=H4.*(H3);
mag=abs(H);% Calculate the absolute value and complex number of H
db=20*log10((mag + eps)/max(mag));
figure(2);
subplot(2,1,1)% divides the current figure into a 2×1 grid and draws the figure at the position specified in the first row.
plot(w/pi,db);
axis([0,1.2,-160,10]);%Change the coordinate axis range so that the X-axis ranges from 0 to 1.2 and the y-axis ranges from -160 to 10.
xlabel('w/π');%x-axis coordinate is w/Π
ylabel('20lg|H(ejw|');%Y-axis coordinate axis is 20lg|H(ejw|
title('Amplitude-frequency response curve of IIR filter designed by bilinear transformation method');%title
% Impulse response invariant method to design ButterWorth digital filter amplitude-frequency characteristics
C1=[0.2871 -0.4466];
C2=[-2.1428 1.1454];
C3=[1.8558 -0.6304];
D1=[1 -0.1297 0.6949];
D2=[1 -1.0691 0.3699];
D3=[1 -0.9972 0.2570];
[E1,~]=freqz(C1,D1,100); % Calculate frequency response freqz digital filter frequency response
[E2,~]=freqz(C2,D2,100); % Calculate frequency response freqz digital filter frequency response
[E3,w]=freqz(C3,D3,100); % Calculate frequency response freqz digital filter frequency response
E4=E1 + E2;
E=E4 + E3;
mag=abs(E);% Calculate the absolute value of E
db=20*log10((mag + eps)/max(mag));
subplot(2,1,2)% divides the current figure into a 2×1 grid and draws the figure at the position specified in the second line.
plot(w/pi,db);% make a two-dimensional line graph in which x is w/pi and y is db.
axis([0,1,-50,10]);%Change the coordinate axis range so that the X-axis ranges from 0 to 1 and the y-axis ranges from -50 to 10.
xlabel('w/π');
ylabel('20lg|H(ejw|');
The knowledge points of the article match the official knowledge files, and you can further learn related knowledge. Algorithm skill tree Home page Overview 57581 people are learning the system