Filter Implementation

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Filter Implementation

shall689
This post was updated on .
Hello,

I have been struggling with trying to implement a low pass filter in excel using filter coefficients generated from octave (see m file text below). The coefficients are
b0 = 0.0074833989
b1 = 0.0074833989
a1 = -0.9850332

The excel file is computing the equation:y[n] = b0*x[n] + b1*x[n-1] + a1*y[n-1]
and filtering the signal: y=2048*sin(2*pi*10000*t)

Excel computes a filtered sinewave having a peak value of 15.6 without any phase shift.

I expected a sinewave with a peak value of 123 with a -89 phase shift.

Is the above equation correct?  Is there something special I have to do?

Thanks,
Stephen

m file
------
pkg load control
pkg load signal
format long
fc = 600;
SampleRate = 250000
fin = 10000;
Tc = 1/(2*pi*fc);
s = tf ('s');
sysC = 1/(1+s*Tc) % continuous transfer functions
sysD = c2d (sysC, 1/SampleRate, 'tustin') % discrete transfer function
[b,a] = tfdata(sysD,'v') % extract the coefficents
InputSignalLengthInSeconds = 10/fin;   % input signal length in seconds
t =linspace(0,InputSignalLengthInSeconds,InputSignalLengthInSeconds*250000);
wt = 2*pi*t*fin;
y=2048*sin(wt);
plot(t,y);
hold on;
y2 = filter(b,a,y);
plot(t,y2);
%print ("FilteredSignal.jpg")
%bode(sysC);
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Filter Implementation

Doug Stewart-4


On Mon, Jul 3, 2017 at 4:39 PM, shall689 <[hidden email]> wrote:
Hello,

I have been struggling with trying to implement a low pass filter in excel
using filter coefficients generated from octave (see m file text below). The
coefficients are
b0 = 0.0074833989
b1 = 0.0074833989
a1 = -0.9850332

The excel file is computing the equation:y[n] = b0*x[n] + b1*x[n-1] +
a1*y[n-1]
and filtering the signal: y=2048*sin(2*pi*10000*t)

Excel computes a filtered sinewave having a peak value of 15.6 without any
phase shift.

I expected a sinewave with a peak value of 123 with a -89 phase shift.

Is the above equation correct?  Is there some special I have to do?

Thanks,
Stephen

m file
------
pkg load control
pkg load signal
format long
fc = 600;
SampleRate = 250000
fin = 10000;
Tc = 1/(2*pi*fc);
s = tf ('s');
sysC = 1/(1+s*Tc) % continuous transfer functions
sysD = c2d (sysC, 1/SampleRate, 'tustin') % discrete transfer function
[b,a] = tfdata(sysD,'v') % extract the coefficents
InputSignalLengthInSeconds = 10/fin;   % input signal length in seconds
t =linspace(0,InputSignalLengthInSeconds,InputSignalLengthInSeconds*250000);
wt = 2*pi*t*fin;
y=2048*sin(wt);
plot(t,y);
hold on;
y2 = filter(b,a,y);
plot(t,y2);
%print ("FilteredSignal.jpg")
%bode(sysC);




I am working on it give my 20 min.


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Filter Implementation

Doug Stewart-4


On Mon, Jul 3, 2017 at 7:09 PM, Doug Stewart <[hidden email]> wrote:


On Mon, Jul 3, 2017 at 4:39 PM, shall689 <[hidden email]> wrote:
Hello,

I have been struggling with trying to implement a low pass filter in excel
using filter coefficients generated from octave (see m file text below). The
coefficients are
b0 = 0.0074833989
b1 = 0.0074833989
a1 = -0.9850332

The excel file is computing the equation:y[n] = b0*x[n] + b1*x[n-1] +
a1*y[n-1]
and filtering the signal: y=2048*sin(2*pi*10000*t)

Excel computes a filtered sinewave having a peak value of 15.6 without any
phase shift.

I expected a sinewave with a peak value of 123 with a -89 phase shift.

Is the above equation correct?  Is there some special I have to do?

Thanks,
Stephen

m file
------
pkg load control
pkg load signal
format long
fc = 600;
SampleRate = 250000
fin = 10000;
Tc = 1/(2*pi*fc);
s = tf ('s');
sysC = 1/(1+s*Tc) % continuous transfer functions
sysD = c2d (sysC, 1/SampleRate, 'tustin') % discrete transfer function
[b,a] = tfdata(sysD,'v') % extract the coefficents
InputSignalLengthInSeconds = 10/fin;   % input signal length in seconds
t =linspace(0,InputSignalLengthInSeconds,InputSignalLengthInSeconds*250000);
wt = 2*pi*t*fin;
y=2048*sin(wt);
plot(t,y);
hold on;
y2 = filter(b,a,y);
plot(t,y2);
%print ("FilteredSignal.jpg")
%bode(sysC);




I am working on it give my 20 min.




--
DAS

ok there must be something wrong in what you are doing in excel

I added a
figure(1) at the top of you M file and
this at the end of you M file.

 t=0:0.001/1000:0.001;
x=2048*sin(2*pi*10000*t);
y(1) = b(1)*x(1) 
for n=2:1001
y(n) = b(1)*x(n) + b(2)*x(n-1) + a(1)*y(n-1);
endfor

figure(2)
plot(t,x,t,y)


using this I see the phase shift etc.

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Filter Implementation

shall689
This post was updated on .
Thanks Doug,


Your work is getting me going in the right direction.  I don't know why I ever used excel in the first place.

You used a(1) (=1) for the y(n-1) term.   I was thinking that a(2) (=-0.985033202259162) should be used for that term, i.e. y(n) = b(1)*x(n) + b(2)*x(n-1) + a(2)*y(n-1).  Is that correct?

If a(2) is used for the y(n-1) term, the plot look like my excel plot, so there is something I am misunderstanding about implementing the filter equation.

Also, the filtered output has an offset if a(1) is used with y(n-1).

Stephen
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Filter Implementation

Doug Stewart-4


On Mon, Jul 3, 2017 at 8:05 PM, shall689 <[hidden email]> wrote:
Thanks Doug,


Your work is getting me going in the right direction.  I don't know why I
was ever used excel in the first place.

You used a(1) (=1) for the y(n-1) term.   I was thinking that a(2)
(=-0.985033202259162) should be used for that term, i.e. y(n) = b(1)*x(n) +
b(2)*x(n-1) + a(2)*y(n-1).  Is that correct?

If a(2) is used for the y(n-1) term, the plot look like my excel plot, so
there is something I am misunderstanding about implementing the filter
equation.

Also, the filtered output has an offset if a(1) is used with y(n-1).

Stephen



you have

y                   b0 *Z +b1
----      =    -----------------------
x                   a0*Z  - a1

y                   b0 *Z +b1                  z^-1
----      =    -----------------------    *    -------
x                   a0*Z  - a1                  z^-1


y                   b0  +b1* z^-1
----      =    -----------------------
x                   a0  - a1 *z^-1



z^-1  means 1 time ago

y*(a0  - a1 *z^-1)= x*( b0  +b1* z^-1)

a0*y=  x*( b0  +b1* z^-1)   + y*a1*z^-1

since octave cant index at 0 we rename them and get my program.



--
DAS


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Filter Implementation

shall689
Fantastic!

It gives the correct response for various frequencies, If I use a1 = 0.985033202259162 instead of -0.985033202259162 in your code.

Thanks alot,
Stephen
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Filter Implementation

Doug Stewart-4


On Mon, Jul 3, 2017 at 8:30 PM, shall689 <[hidden email]> wrote:
Fantastic!

It gives the correct response for various frequencies, If I use a1 =
0.985033202259162 instead of -0.985033202259162 in your code.

Thanks alot,
Stephen



--
View this message in context: http://octave.1599824.n4.nabble.com/Filter-Implementation-tp4683979p4683992.html
Sent from the Octave - General mailing list archive at Nabble.com.

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


Yes  you are correct. I had incorporated the minus sign into my code.

the correct code is


 t=0:0.001/1000:0.001;
x=2048*sin(2*pi*10000*t);
y(1) = b(1)*x(1) ;
for n=2:1001
y(n) = b(1)*x(n) + b(2)*x(n-1) - a(1)*y(n-1);
endfor

figure(2)
plot(t,x,t,y)

--
DASCertificate for 206392


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Loading...