Help implementing a digital filter

I'm trying to implement a digital filter. I have the coefficients. I make the calculation in a for-loop and it's not working. The results are aberrant = 10^307.

The strange thing is that when I use the BUTTER function and use the returned coefficients, everything works correctly. So, the algorithm is correct.

If I just take the coefficients returned by butter() from workspace and use them in another *.m file as normal variables, it's not working. Why?

If I want to calculate the coefficients and check them in this way, why I cannot do it?

[Todd Flanagan]

Hi. Posting a code snippet would be very helpful.

[Monica]

Hi!

This is the way it's working. If I uncomment lines 62-63 (give the coefficients a,b it's not working - this is what i

I need, coz I have to compute my coefficients depending on some electronic components)

Thanks!

code:

clear all;

clc;

magnitude_LL =[ 400

400

400];

A12=magnitude_LL(1)*sqrt(2)/400;

A23=magnitude_LL(2)*sqrt(2)/400;

A31=magnitude_LL(3)*sqrt(2)/400;

arg_LL =[ 30

-90

150];

fi12_degrees=arg_LL(1);

fi12=fi12_degrees*pi/180; % radians

fi23_degrees=arg_LL(2);

fi23=fi23_degrees*pi/180; % radians

fi31_degrees=arg_LL(3);

fi31=fi31_degrees*pi/180; % radians

%A12=1;

N12=2; %harmonic

B12=0; %amplitude of the N12th harmonic

A31=0;

N31=5; %harmonic

B31=0.5; %amplitude of the N31th harmonic

B23=0;

%-------------------------------------------------------------------------

f=50;

T=1/f;

step=T/1024;

t=0:step:1.5; %simulation time

%input signals

U12=A12.*sin(2*pi*f.*t+fi12)+B12.*sin(2*pi*(f*N12).*t+fi12);

U31=A31.*sin(2*pi*f.*t+fi31)+B31.*sin(2*pi*(f*N31).*t+fi31);

U23=A23.*sin(2*pi*f.*t)+B23.*sin(2*pi*(f*N12).*t);

fc=60;

nr_samples=1024;

fc_normal=2*fc*step;

[b,a] = butter(5,fc_normal,'low')

figure;

freqz(b,a,5120,1/step)

a0= a(1);

a1= a(2);

a2= a(3);

a3= a(4);

a4= a(5);

a5= a(6);

b0= b(1);

b1= b(2);

b2= b(3);

b3= b(4);

b4= b(5);

b5= b(6);

% b = 1.0e-011 *[0.0668 0.3342 0.6683 0.6683 0.3342 0.0668]

% a =[ 1.0000 -4.9762 9.9050 -9.8579 4.9055 -0.9765]

out_PLL_12=cos(2*pi*f.*t);

s=U12.*out_PLL_12;

figure('name','test2');

plot(t,U12);hold on

plot(t,out_PLL_12,'r');

for i=1:length(s)

if i>=6

s_filtered(i)=(1/a0)*(b0*s(i)+b1*s(i-1)+b2*s(i-2)+b3*s(i-3)+b4*s(i-4)+b5*s(i-5)...

-a1*s_filtered(i-1)-a2*s_filtered(i-2)-a3*s_filtered(i-3)-a4*s_filtered(i-4)-a5*s_filtered(i-5));

else

s_filtered(i)=0;

end

end

figure;

plot(t,out_PLL_12);hold on;grid on;

plot(t,s,'g');hold on; grid on;

plot(t,s_filtered,'r');hold on; grid on;

Answer 1

[Vieniava]

The problem is that those lines

b = 1.0e-011 *[0.0668 0.3342 0.6683 0.6683 0.3342 0.0668]
a =[ 1.0000 -4.9762 9.9050 -9.8579 4.9055 -0.9765]

are only approximation (for display purpose) of values computed by butter() .

To store your coefficients you should after

[b,a] = butter(5,fc_normal,'low')

write this

save('MyFilter','a','b')

When you need restore those values in another m-file just load via

load MyFilter

[Monica]

Thanks a lot for the answer!

So, as I understand it, I won't be able to see the exact coefficients. I'm using MatLab just to test the filter, but in the end I must implement it on a DSP, so in the code I must give the filter coefficients. I don't have experience with DSPs so far, but I believe I must work with basic mathematical functions.

Many thanks again,

Monica

[joanna]

Pleasure to read so clear answers:)

[Bhaskar]

Use

>> format long

to see results with more number of decimal points.