Digital filtering in matlab not giving expected results

Question

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!

Answer 1

I found the problem, which is that the numden() function of matlab doesn't give the right numerators and denominators. thanks to the people who took time to help me!

the following code gives a wrong result:

z=sym('z');
Pi=sym('Pi');
Ts=sym('Ts');
vc=sym('vc');

Fz=1/((Pi^2*Ts^2*vc^2)/(z^2 - 2*z + 1) + (2^(1/2)*Pi*Ts*vc)/(z - 1) + (2*Pi^2*Ts^2*vc^2*z)/(z^2 - 2*z + 1) + (Pi^2*Ts^2*vc^2*z^2)/(z^2 - 2*z + 1) + (2^(1/2)*Pi*Ts*vc*z)/(z - 1) + 1)

[num, den]=numden(Fz);
cn=coeffs(num,z)
cd=coeffs(den,z)

matlab gives the coefficients:

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]

which are obvisoulsy wrong

thanks!

Answer 2

The result is probably correct. Your sine wave frequency is 500 and your sampling frequency is 10000. Your normalized frequency is 500/10000=0.05 which is filtered. So you see the highly attenuated output. You should increase your sine wave frequency and see whether you get your sine wave back after filtering.