I'm starting with chapters 5 and 6 of the 2nd edition of Rick's _Understanding Digital Signal Processing_ . I have no problem with explanations using summation operator ( CAP SIGMA ). How do I filter a digitized signal? I *think* I'm looking for a matrix representation. I'm using Scilab. I have an extremely large linear vector as input. How do I apply filter? Size input vector >>>>> size matrix of filter coefficients. Or am I asking "wrong question"? Any pointers appreciated.
[NEWBIE question] How do I filter a digitized signal
Started by ●August 13, 2004
Reply by ●August 13, 20042004-08-13
Richard Owlett <rowlett@atlascomm.net> writes:> I'm starting with chapters 5 and 6 of the 2nd edition of Rick's > _Understanding Digital Signal Processing_ . > > > I have no problem with explanations using summation operator > ( CAP SIGMA ). > > How do I filter a digitized signal? > I *think* I'm looking for a matrix representation. > > I'm using Scilab. > I have an extremely large linear vector as input. > How do I apply filter? > Size input vector >>>>> size matrix of filter coefficients. > > Or am I asking "wrong question"? > > Any pointers appreciated.Hi Richard, The linear-algebra view of an FIR filter is simply this: you form an "inner product" of the "data vector" at time n (call it X[n]) with the "coefficient vector" (call it H). Let's say both of these are column vectors of length N, where N is the number of coefficients. So they're both Nx1 vectors. Then X[n] = [x[n] x[n-1] x[n-2] ... x[n - N + 1]]' (where "'" denotes transpose) and H = [h[0] h[1] ... h[N-1]]'. Then the output of the filter at time n, y[n], is simply the inner product of X[n] and H, i.e., y[n] = X[n]' * H. Note the dimensions agree - X[n]' is 1xN, H is Nx1, and the product is 1x1. Is that what you were looking for? -- Randy Yates Sony Ericsson Mobile Communications Research Triangle Park, NC, USA randy.yates@sonyericsson.com, 919-472-1124
Reply by ●August 13, 20042004-08-13
Richard Owlett wrote:> I'm starting with chapters 5 and 6 of the 2nd edition of Rick's > _Understanding Digital Signal Processing_ . > > I have no problem with explanations using summation operator > ( CAP SIGMA ). > > How do I filter a digitized signal? > I *think* I'm looking for a matrix representation. > > I'm using Scilab. > I have an extremely large linear vector as input. > How do I apply filter? > Size input vector >>>>> size matrix of filter coefficients. > > Or am I asking "wrong question"? > > Any pointers appreciated. >In Scilab flts works nicely for an IIR, and would work (albeit inefficiently) for an FIR -- just describe it as the filter vector over z^N, and it will treat it like any other transfer function. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
Reply by ●August 13, 20042004-08-13
Randy Yates wrote:> Richard Owlett <rowlett@atlascomm.net> writes: > > >>I'm starting with chapters 5 and 6 of the 2nd edition of Rick's >>_Understanding Digital Signal Processing_ . >> >> >>I have no problem with explanations using summation operator >>( CAP SIGMA ). >> >>How do I filter a digitized signal? >>I *think* I'm looking for a matrix representation. >> >>I'm using Scilab. >>I have an extremely large linear vector as input. >>How do I apply filter? >>Size input vector >>>>> size matrix of filter coefficients. >> >>Or am I asking "wrong question"? >> >>Any pointers appreciated. > > > Hi Richard, > > The linear-algebra view of an FIR filter is simply this: you form an > "inner product" of the "data vector" at time n (call it X[n]) with the > "coefficient vector" (call it H). Let's say both of these are column > vectors of length N, where N is the number of coefficients. So they're > both Nx1 vectors. > > Then X[n] = [x[n] x[n-1] x[n-2] ... x[n - N + 1]]' (where "'" denotes > transpose) and H = [h[0] h[1] ... h[N-1]]'. Then the output of the > filter at time n, y[n], is simply the inner product of X[n] and H, i.e., > > y[n] = X[n]' * H. > > Note the dimensions agree - X[n]' is 1xN, H is Nx1, and the > product is 1x1. Is that what you were looking for?I think I was going in that direction. input(1)*h(1) + input(2)*h(2) + ... + input(n+0)*h(n) = OUTPUT(1) input(2)*h(1) + input(3)*h(2) + ... + input(n+1)*h(n) = OUTPUT(2) input(3)*h(1) + input(4)*h(2) + ... + input(n+2)*h(n) = OUTPUT(3) | repeat M-N times ( M=length input vector ) That handles FIR case. What to do with IIR where h is N x2 ?
Reply by ●August 13, 20042004-08-13
Tim Wescott wrote:> Richard Owlett wrote: > >> I'm starting with chapters 5 and 6 of the 2nd edition of Rick's >> _Understanding Digital Signal Processing_ . >> >> I have no problem with explanations using summation operator >> ( CAP SIGMA ). >> >> How do I filter a digitized signal? >> I *think* I'm looking for a matrix representation. >> >> I'm using Scilab. >> I have an extremely large linear vector as input. >> How do I apply filter? >> Size input vector >>>>> size matrix of filter coefficients. >> >> Or am I asking "wrong question"? >> >> Any pointers appreciated. >> > In Scilab flts works nicely for an IIR, and would work (albeit > inefficiently) for an FIR -- just describe it as the filter vector over > z^N, and it will treat it like any other transfer function. >Sorry, but that lost me. Suspect I'm missing something obvious.
Reply by ●August 13, 20042004-08-13
Richard Owlett wrote:> Randy Yates wrote: > >> Richard Owlett <rowlett@atlascomm.net> writes: >> >> >>> I'm starting with chapters 5 and 6 of the 2nd edition of Rick's >>> _Understanding Digital Signal Processing_ . >>> >>> >>> I have no problem with explanations using summation operator >>> ( CAP SIGMA ). >>> >>> How do I filter a digitized signal? >>> I *think* I'm looking for a matrix representation. >>> >>> I'm using Scilab. >>> I have an extremely large linear vector as input. >>> How do I apply filter? >>> Size input vector >>>>> size matrix of filter coefficients. >>> >>> Or am I asking "wrong question"? >>> >>> Any pointers appreciated. >> >> >> >> Hi Richard, >> >> The linear-algebra view of an FIR filter is simply this: you form an >> "inner product" of the "data vector" at time n (call it X[n]) with the >> "coefficient vector" (call it H). Let's say both of these are column >> vectors of length N, where N is the number of coefficients. So they're >> both Nx1 vectors. >> >> Then X[n] = [x[n] x[n-1] x[n-2] ... x[n - N + 1]]' (where "'" denotes >> transpose) and H = [h[0] h[1] ... h[N-1]]'. Then the output of the >> filter at time n, y[n], is simply the inner product of X[n] and H, i.e., >> >> y[n] = X[n]' * H. >> >> Note the dimensions agree - X[n]' is 1xN, H is Nx1, and the >> product is 1x1. Is that what you were looking for? > > > I think I was going in that direction. > > input(1)*h(1) + input(2)*h(2) + ... + input(n+0)*h(n) = OUTPUT(1) > input(2)*h(1) + input(3)*h(2) + ... + input(n+1)*h(n) = OUTPUT(2) > input(3)*h(1) + input(4)*h(2) + ... + input(n+2)*h(n) = OUTPUT(3) > | > repeat M-N times ( M=length input vector ) > > That handles FIR case. > What to do with IIR where h is N x2 ? > >I'll have to get a copy of the book, but in the meantime have you gotten the concept of a transfer function yet? If not see http://www.wescottdesign.com/articles/zTransform/z-transforms.html or read ahead to the appropriate part of Rick's book. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
Reply by ●August 13, 20042004-08-13
Tim Wescott wrote:> Richard Owlett wrote: > >> Randy Yates wrote: >> >>> Richard Owlett <rowlett@atlascomm.net> writes: >>> >>> >>>> I'm starting with chapters 5 and 6 of the 2nd edition of Rick's >>>> _Understanding Digital Signal Processing_ . >>>> >>>> >>>> I have no problem with explanations using summation operator >>>> ( CAP SIGMA ). >>>> >>>> How do I filter a digitized signal? >>>> I *think* I'm looking for a matrix representation. >>>> >>>> I'm using Scilab. >>>> I have an extremely large linear vector as input. >>>> How do I apply filter? >>>> Size input vector >>>>> size matrix of filter coefficients. >>>> >>>> Or am I asking "wrong question"? >>>> >>>> Any pointers appreciated. >>> >>> >>> >>> >>> Hi Richard, >>> >>> The linear-algebra view of an FIR filter is simply this: you form an >>> "inner product" of the "data vector" at time n (call it X[n]) with the >>> "coefficient vector" (call it H). Let's say both of these are column >>> vectors of length N, where N is the number of coefficients. So they're >>> both Nx1 vectors. >>> >>> Then X[n] = [x[n] x[n-1] x[n-2] ... x[n - N + 1]]' (where "'" denotes >>> transpose) and H = [h[0] h[1] ... h[N-1]]'. Then the output of the >>> filter at time n, y[n], is simply the inner product of X[n] and H, i.e., >>> >>> y[n] = X[n]' * H. >>> >>> Note the dimensions agree - X[n]' is 1xN, H is Nx1, and the >>> product is 1x1. Is that what you were looking for? >> >> >> >> I think I was going in that direction. >> >> input(1)*h(1) + input(2)*h(2) + ... + input(n+0)*h(n) = OUTPUT(1) >> input(2)*h(1) + input(3)*h(2) + ... + input(n+1)*h(n) = OUTPUT(2) >> input(3)*h(1) + input(4)*h(2) + ... + input(n+2)*h(n) = OUTPUT(3) >> | >> repeat M-N times ( M=length input vector ) >> >> That handles FIR case. >> What to do with IIR where h is N x2 ? >> >> > I'll have to get a copy of the book, but in the meantime have you gotten > the concept of a transfer function yet?That's not a problem. I completed ~3 yrs towards a BSEE back in 60's Now I have to convert from analog to digital domain > If not see> http://www.wescottdesign.com/articles/zTransform/z-transforms.html or > read ahead to the appropriate part of Rick's book. >I've read ahead and I'll take a look at your page. I think I've been there before. But I know that I do not get everything the first time ;]
Reply by ●August 14, 20042004-08-14
Richard Owlett <rowlett@atlascomm.net> wrote in message news:<10hprc8pu4va64e@corp.supernews.com>...> I'm starting with chapters 5 and 6 of the 2nd edition of Rick's > _Understanding Digital Signal Processing_ . > > I have no problem with explanations using summation operator > ( CAP SIGMA ). > > How do I filter a digitized signal? > I *think* I'm looking for a matrix representation.Well, FIR filters are easy to implement in terms of linear algebra. Implementing IIR filters this way will definately be a lot more cumbersome, if at all possible.> I'm using Scilab. > I have an extremely large linear vector as input. > How do I apply filter? > Size input vector >>>>> size matrix of filter coefficients. > > Or am I asking "wrong question"?Nope, definately the right question. I'll run through the matrix implementation for the FIR filter, using matlab notation, as always. %%%% Start Matlab Pseudo Code %%%%%%%%%%%%%%%%%%%%%%%%% % % Preliminaries: % % x - Data column vector, Nx1 elements % h - Filter impulse response, Mx1 column vector % Prepare space for data matrix: xx=zeros(N+M-1,M); % Insert the data vector in each column of xx. Delay one sample % between each column. for m=1:M xx(m:m+N-1,m)=x; end % Filter data y = xx*h; %%%% End Matlab Pseudo Code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Here y is a vector of length N+M-1, containing your filtered data. Note that the matrix xx has more than NxM elements, which can be a very large number. Now, an off-line filtering method that goes easier on the memory on your computer, is to work in frequency domain. For FIR filters, this is very easy, since the filter coefficients are also the impulse response. For IIR filters, you need to compute the impulse response of the filter, i.e. to implement eq. (6.2) in Rick's book with an impulse sequence as input. This might be a bit cumbersome but isn't really difficult. %%%% Start Matlab Pseudo Code %%%%%%%%%%%%%%%%%%%%%%%%% % % Preliminaries: % % x - Data column vector, Nx1 elements % h - Filter impulse response, Mx1 column vector %%%% For IIR filters, compute the impulse response f here. % Take FFT of x and h. In case of IIR filters, M is the "efficient % length" of the impulse response. Zero-pad h and x to length N+M-1 % to make space for end transients. In practice, use the lowest % power of 2 > N+M-1. H=fft(h,N+M-1); X=fft(X,N+M-1); % Filter data as spectrum multiplication Y=X*H; % Find resulting time series as inverse FFT y = real(ifft(Y)); %%%% End Matlab Pseudo Code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Rune
Reply by ●August 14, 20042004-08-14
On Fri, 13 Aug 2004 16:43:51 -0500, Richard Owlett <rowlett@atlascomm.net> wrote: (snipped)> >I think I was going in that direction. > >input(1)*h(1) + input(2)*h(2) + ... + input(n+0)*h(n) = OUTPUT(1) >input(2)*h(1) + input(3)*h(2) + ... + input(n+1)*h(n) = OUTPUT(2) >input(3)*h(1) + input(4)*h(2) + ... + input(n+2)*h(n) = OUTPUT(3) >| >repeat M-N times ( M=length input vector ) > >That handles FIR case. >What to do with IIR where h is N x2 ?Hi Richard, I'm not sure I can answer your question, but I'll give it a try. My Eq. (5-4), on page 158, gives the necessary computations to implement a tapped-delay line FIR filter. That equation is the infamous "convolution equation" as it pertains to nonrecursive (tapped-delay line) FIR filters. I tried to introduce that equation as slowly, and gently, as I could. For various values of time-index n, Eq. (5-4) for a four-tap filter goes like this: n=3: y(3) = h3*x(3) + h2*x(2) + h1*x(1) + h0*x(0) n=4: y(4) = h3*x(4) + h2*x(3) + h1*x(2) + h0*x(1) n=5: y(5) = h3*x(5) + h2*x(4) + h1*x(3) + h0*x(2) n=6: y(6) = h3*x(6) + h2*x(5) + h1*x(4) + h0*x(3) etc. So, to implement nonrecursive (tapped-delay line) FIR filters, your computations are sum the products to compute an output sample shift input sequence relative to filter coeffs sum the products to compute an output sample shift input sequence relative to filter coeffs sum the products to compute an output sample shift input sequence relative to filter coeffs sum the products to compute an output sample etc. I don't know about Scilab, but MATLAB has a "filter" command that's used to filter some time-domain "x" variable (a list of input sample values) using a filter that's defined by its feedforward and feedback coefficients. For nonrecursive FIR filters, the feedforward coefficients needed by MATLAB are the FIR filter's coefficients (h), and the feedback coefficients needed by MATLAB turn out to be a simple single value of one. In MATLAB language, filtering an input sequence, x, with an FIR filter looks like: Output = filter([h0,h1,h2,h3],1,x); Maybe Scilab is similar to the above, I don't know. That above MATLAB command implements my Eq. (6-1) on page 213. [Please note that there is a "typo" in Eq. (6-1) that I'm desperately trying to get the publisher to correct. In Eq. (6-1) the term "h(4)x(n-4)" should NOT be there. Please cross out that h(4)x(n-4) term.] For an IIR filter, my Eq. (6-2), on page 214, gives the necessary computations to implement the IIR filter in Figure 6-3. In MATLAB language, filtering an input sequence, x, with the Figure 6-3 IIR filter looks like: Output = filter([b(0),b(1),b(2),b(3)],[1,-a(1),-a(2),-a(3)],x); No, the minus signs for the above a(k) feedback coefficients are not a mistake. The minus signs are necessary because of the way MATLAB expects the feedback coefficients to be defined. If you want to implement IIR filtering using some sort of "for" loop, you need to define variables for the instantaneous values at the output of each delay element in Figure 6-3. Then those variables are updated in preparation for computing a new filter output sample. Here's an example referring to my Figure 6-3: Let's call the x(n-1) value by the name "x1", call the x(n-2) value by the name "x2", call the x(n-3) value by the name "x3", and let's call the y(n-1) value by the name "y1", call the y(n-2) value by the name "y2", call the y(n-3) value by the name "y3". OK, next we set those six variables all to zero. That's our "initial condition" of the IIR filter. Now we dive into a "computational loop", like this, to compute our first IIR filter output using my Eq. (6-2): y(0) = b0*x(0) + b1*x1 + b2*x2 + b3*x3 + a1*y1 + a2*y2 + a3*y3. That y(0) is our 1st filter output value. Next we prepare for a new computation by updating the variables representing the outputs of the filter's delay elements. That "updating" is: x1 = x(0), x2 = x1, x3 = x2, y1 = y(0), y2 = y1, y3 = y2. Now we're ready to use the next input sample, x(1), to compute the next output sample y(1). That computation is a repeat (we're in a computational loop, right?) of the above, which we write as: y(1) = b0*x(1) + b1*x1 + b2*x2 + b3*x3 + a1*y1 + a2*y2 + a3*y3. Using the time-domain index n, we can say our filter computation inside the "for" loop is: y(n) = b0*x(n) + b1*x1 + b2*x2 + b3*x3 + a1*y1 + a2*y2 + a3*y3. Right after computing a y(n) output sample value, we update the x1, x2, x3, y1, y2, & y3 values in preparation for computing the next output y(n+1). We stop looping when we run out of x(n) inputs samples to filter. Sheece! I've rambled on like crazy here. Am hoping the above helps ya' a little Richard. [-Rick-] Where am I? In the Village. What do you want? Information. Whose side are you on? That would be telling ... We want Information. You won't get it. By hook or by crook ... we will. Who are you? The new Number Two. Who is Number One? You are Number Six. I am not a number ... I am a free man! (Mocking laughter)
Reply by ●August 14, 20042004-08-14
Rick Lyons wrote: ...> Hi Richard,...> For various values of time-index n, Eq. (5-4) for a > four-tap filter goes like this: > > n=3: y(3) = h3*x(3) + h2*x(2) + h1*x(1) + h0*x(0) > n=4: y(4) = h3*x(4) + h2*x(3) + h1*x(2) + h0*x(1) > n=5: y(5) = h3*x(5) + h2*x(4) + h1*x(3) + h0*x(2) > n=6: y(6) = h3*x(6) + h2*x(5) + h1*x(4) + h0*x(3) > etc....> Eq. (6-2): > > y(0) = b0*x(0) + b1*x1 + b2*x2 + b3*x3 > + a1*y1 + a2*y2 + a3*y3. > > That y(0) is our 1st filter output value. Next we > prepare for a new computation by updating the variables > representing the outputs of the filter's delay > elements. That "updating" is: > > x1 = x(0), x2 = x1, x3 = x2, > y1 = y(0), y2 = y1, y3 = y2. > > Now we're ready to use the next input sample, x(1), > to compute the next output sample y(1). That computation > is a repeat (we're in a computational loop, right?) > of the above, which we write as: > > y(1) = b0*x(1) + b1*x1 + b2*x2 + b3*x3 > + a1*y1 + a2*y2 + a3*y3. > > Using the time-domain index n, we can say > our filter computation inside the "for" loop is: > > y(n) = b0*x(n) + b1*x1 + b2*x2 + b3*x3 > + a1*y1 + a2*y2 + a3*y3. > > Right after computing a y(n) output sample value, > we update the x1, x2, x3, y1, y2, & y3 values in > preparation for computing the next output y(n+1).... Why does x(0) have parentheses that x1 etc. lacks? Is there a pattern? Why in the reference to Eq. (5-4) are there parentheses throughout? My best guess is that it's late and you're tired, but I have to ask just in case. I hope you don't mind. Jerry -- ... the worst possible design that just meets the specification - almost a definition of practical engineering. .. Chris Bore ������������������������������������������������������������������������
