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); |
On Mon, Jul 3, 2017 at 4:39 PM, shall689 <[hidden email]> wrote:
Hello, _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
On Mon, Jul 3, 2017 at 7:09 PM, Doug Stewart <[hidden email]> wrote:
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 |
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 |
On Mon, Jul 3, 2017 at 8:05 PM, shall689 <[hidden email]> wrote: you haveThanks Doug, 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. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
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 |
On Mon, Jul 3, 2017 at 8:30 PM, shall689 <[hidden email]> wrote: Fantastic! 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) _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave |
Free forum by Nabble | Edit this page |