I'm trying to design and test a simple digital high pass filter in matlab.
I have two scripts: the first one is to design the filter, the second one is the implementation fo the recursion algorithm
1st script: desgn of the filter:
s=sym('s');
z=sym('z');
w=sym('w');
p=sym('p');
f=sym('f');
x=sym('x');
Pi=sym('Pi');
Ts=sym('Ts');
vc=sym('vc');
n=2;
Fspb=1/(s^2+sqrt(2)*s+1); % simple low pass butterworth
Fs=subs(Fspb, s, 2*Pi*vc/s); % tranform to high pass butterworth with vc cutting freq
Fz=subs(Fs, s, 2/Ts*(z-1)/(z+1)); %bilinear transformation
Fp=subs(Fz, z, exp(Ts*p))
Fw=subs(Fp, p, 1i*w);
Ff=subs(Fw, w, 2*Pi*f);
vmax=5000;
vc_val=1000; %1000 Hz
vccont=1/Pi/Ts*tan(Pi*Ts*vc_val);
Ts_val=0.0001
fval=0:0.1:vmax;
pretty(expand(Fz))
[num, den]=numden(Fz);
cn=coeffs(num,z)
cd=coeffs(den,z)
it gives me the coefficients (for z) and the frequency response of the filter:
cn =
[ 1, -2, 1]
cd =
[ Pi^2*Ts^2*vc^2 - 2^(1/2)*Pi*Ts*vc + 1, 2*Pi^2*Ts^2*vc^2 - 2, Pi^2*Ts^2*vc^2 + 2^(1/2)*Pi*Ts*vc + 1]
here is the mfilter displayed with the freqz function:
and my second script, implementing the filter and testing on a simple sinus function:
Pi=3.14;
v=500;
Ts=0.0001;
x=0:Ts:100000*Ts;
y=sin(x*2*Pi*v); % sinus, freq=v
figure;
plot(x,y);
%% filter def
vc_val=1; %1000 Hz
vccont=1/Pi/Ts*tan(Pi*Ts*vc_val);
a=Pi^2*Ts^2*vccont^2;
b=sqrt(2)*Pi*Ts*vccont;
%% filtering
yf=zeros(1,size(y,2));
y1=0; y2=0; y3=0; y4=0; x1=0; x2=0; x3=0; x4=0;
for i=3:size(y,2)
if (i>1)
y1=yf(i-1);
x1=xyi-1);
if (i>2)
y2=yf(i-2);
x2=y(i-2);
end
end
yf(i)=1/(a-b+1)*(-(2*a-2)*y1-(a+b+1)*y2+y(i)-2*x1+x2);
end
figure;
plot(x,yf);
But it doesn't give me the result i was expecting (left: the original sinus, right: the filtered result at 1/2 of the cutting frequency):
I checked that my filter is stable, I can't see what the problem may be. Have you got any hint?
thanks!