Filter question: Problem with time-varying filter

Hello,

I am going to implement a filter on a DSP. Before I do that I want to check 
that the filter
works fine in matlab.

So I have made a filter function in MATLAB:

function 
[y,new_DelayLine1,new_DelayLine2]=filtx(x,num_b,den_a,old_DelayLine1,old_DelayLine2)

where x is a row vector of samples.

num_b and den_a are the numerator and denominator filter coefficients.

If I generate a vector x with random data and use MATLABs filter and my own 
function I get the same
output. As soon as I change the filter coefficients the output from MATLABs 
function and my own is no
longer the same.

How do I solve this?

Here is the function:

function 
[y,new_DelayLine1,new_DelayLine2]=filtx(x,num_b,den_a,old_DelayLine1,old_DelayLine2)

K=size(x,2);
M=size(num_b,2);
N=size(den_a,2);

for k=1:K

% Filter (FIR)
old_DelayLine1(1,1)=x(1,k);
sum1=x(1,k)*num_b(1,1);
for m=2:M
    sum1=sum1+old_DelayLine1(1,m)*num_b(1,m);
end
% Update FIR delay line
for m=M-1:-1:1
    old_DelayLine1(1,m+1)=old_DelayLine1(1,m);
end


% Filter (IIR)
sum2=0;
for n=2:N
    sum2=sum2+old_DelayLine2(1,n)*den_a(1,n);
end
% Filter output
y(1,k)=(sum1-sum2)/den_a(1,1);

% Update IIR delay line
old_DelayLine2(1,1)=y(1,k);
for n=N-1:-1:1
    old_DelayLine2(1,n+1)=old_DelayLine2(1,n);
end

end
new_DelayLine1=old_DelayLine1;
new_DelayLine2=old_DelayLine2;



And here is the test program:

clc
close all
clear

% generate new samples
x=randn(1,10);
% define filter
[num_b,den_a]=butter(4,0.2);

% reset delay lines used in my function
new_DelayLine1=zeros(1,size(num_b,2));
new_DelayLine2=zeros(1,size(den_a,2));

% reset filter state for MATLABs filter function
s=zeros(1,4);

% Use my function to calculate filter output
[y,new_DelayLine1,new_DelayLine2]=filtx(x,num_b,den_a,new_DelayLine1,new_DelayLine2);
% Use MATLABs filter function
[z,s]=filter(num_b,den_a,x,s);

y % OK
z % OK

% Check that functions perform the same way when x is varying

% new samples
x=randn(1,10);

% Use my function to calculate filter output
[y,new_DelayLine1,new_DelayLine2]=filtx(x,num_b,den_a,new_DelayLine1,new_DelayLine2);
% Use MATLABs filter function
[z,s]=filter(num_b,den_a,x,s);

y % OK
z % OK


% Check that functions perform the same way when x is varying AND filter is 
varying

% new samples
x=randn(1,10);

% update filter coefficients
[num_b,den_a]=butter(4,0.6);


% filter
[y,new_DelayLine1,new_DelayLine2]=filtx(x,num_b,den_a,new_DelayLine1,new_DelayLine2);
[z,s]=filter(num_b,den_a,x,s);


% Output y and z is not the same anymore
y
z

Hmmm...I guess it takes time for the filter to converge towards zero. 
Therfore, I must add
the extra non-zero samples from the filter-output to the filter-output for 
the next block of
samples. Right?

"James" <james@nadaspam.thanx> wrote in message 
news:43eb619e$0$78279$157c6196@dreader1.cybercity.dk...
> Hmmm...I guess it takes time for the filter to converge towards zero. 
> Therfore, I must add
> the extra non-zero samples from the filter-output to the filter-output for 
> the next block of
> samples. Right?
>
James,

That's a heck of a lot of code to deal with.

What are you doing *in English*?  You might get more responses.

Fred 

>
> That's a heck of a lot of code to deal with.
>
> What are you doing *in English*?  You might get more responses.
>
> Fred


Sorry about that.

Let me try in english then...

I am trying to implement a filter in MATLAB:

H[z]=(b[0]+b[1]*z^-1+.....+b[N]*z^-N) /  (a[0]+a[1]*z^-1+.....+a[M]*z^-M)

that will act like MATLABs own filter-function:

[y,s]=filter(num,den,x,s)

I have transfered the above function from the Z-domain to the time-domain 
and
implemented the resulting expression as a function in MATLAB.

It works fine as long as the filter-coefficients don't change from one block 
of
samples to the next block of samples.

However, as soon as I change the filter-coefficients from one block to 
another
I get one result using MATLABs filter function and another result using my
own function.

I suspect it is because I have forgotten to implement some code that adds
the transient signal from the filtering of block j to the filter output when
processing block j+1....but I am not sure...

"James" <james@nadaspam.thanx> wrote in message 
news:43eb7b4e$0$78287$157c6196@dreader1.cybercity.dk...
>
>>
>> That's a heck of a lot of code to deal with.
>>
>> What are you doing *in English*?  You might get more responses.
>>
>> Fred
>
>
> Sorry about that.
>
> Let me try in english then...
>
> I am trying to implement a filter in MATLAB:
>
> H[z]=(b[0]+b[1]*z^-1+.....+b[N]*z^-N) /  (a[0]+a[1]*z^-1+.....+a[M]*z^-M)
>
> that will act like MATLABs own filter-function:
>
> [y,s]=filter(num,den,x,s)
>
> I have transfered the above function from the Z-domain to the time-domain 
> and
> implemented the resulting expression as a function in MATLAB.
>
> It works fine as long as the filter-coefficients don't change from one 
> block of
> samples to the next block of samples.
>
> However, as soon as I change the filter-coefficients from one block to 
> another
> I get one result using MATLABs filter function and another result using my
> own function.
>
> I suspect it is because I have forgotten to implement some code that adds
> the transient signal from the filtering of block j to the filter output 
> when
> processing block j+1....but I am not sure...
>

James,

I think you have the right idea.

Note that the matlab function returns the initial transient (and the ending 
transient if you ask for it).  You don't need the initial transient part at 
all if you're working with a "streaming" set of values for x.

It looks like this:

X1        X2         X3         X4         X5         X6         X7 
X8         X9
0 0 0 0 1  0 0 0 1 1  0 0 1 1 1  0 1 1 1 1  1 1 1 1 1  1 1 1 1 0  1 1 1 0 0 
1 1 0 0 0  1 0 0 0 0

Using a purely non-recursive structure so a=1 and let's have b=[1 1 1]

Y1        Y2         Y3         Y4         Y5         Y6         Y7 
Y8         Y9
0 0 0 0 1  0 0 0 1 2  0 0 1 2 3  0 1 2 3 3  1 2 3 3 3  1 2 3 3 2  1 2 3 2 1 
1 2 2 1 0  1 1 1 0 0
YF        YF         YF         YF         YF         YF         YF 
YF         YF
1 1        2 1        2 1        2 1        2 1        1 0        0 0 
0 0        0 0

Here, each block for x represents a single sample shift in x so each block 
in y represents a single sample shift in y.  I've also limited the sequence 
in x to a total length of 9.  You could also look at x and y like this:
b=[1 1 1]
x=[0 0 0 0 1 1 1 1 1 0 0 0 0]
y=[0 0 0 0 1 2 3 3 3 2 1 0 0]

I think this means the following if you're using shorter blocks of x:
Throw out the startup transient in all the Y vectors.
Now we have:
Y1        Y2         Y3         Y4         Y5         Y6         Y7 
Y8
0 0 1      0 1 2      1 2 3      2 3 3      3 3 3      3 3 2      3 2 1 
1 0 0
y=[1 2 3 3 3 2 1 0 0]

This type of filter implementation is probably fine for dealing with long 
arrays of numbers.  It's obviously complicated when dealing with a sequence 
of short blocks - for whatever reason you choose to do that.

You also might ask, what if the blocks aren't overlapping but, rather, 
contiguous:

X1        X2         X3
0 0 0 0 1  1 1 1 1 0  0 0 0 0 0
for x=[0 0 0 0 1 1 1 1 1 0 0 0 0 0 0]

Take a look at the Ys and YFs:
Y1  00001 11  *   > 001
Y2  00012 21
Y3  00123 21
Y4  01233 21  *   > 233
Y5  12333 21
Y6  12332 10
Y7  12321 00  *   > 321
Y8  12210 00
Y9  11100 00
Y10 11000 00  *   > 000
Y11 10000 00
Which parts of which ones should/can be thrown out:

Throw out all the startup transients and compute an output for every N-L+1 
(5-3+1=3) input block meaning that the overlap in input is the same as the 
length of the startup transient L-1=2 here.
Now, that makes sense because there has to be enough extra input data to 
just make up the startup transient that will be thrown out.

Fred 

Thanks Fred.  I am not sure if it answers my question though.
However, I would be very happy if you could tell
me how matlab calculates the final conditions?

Using MATLAB notation:

b=1
a=[0.5 0.2 0.1]
s=[0 0]
x=[1 2 3 4]
[y,s]=filter(b,a,x,s)

Result:

y=2 3.2 4.32 5.632
s=-3.1168 -1.1264

The difference equation for the above filter is:

y(n)=[x(n)-0.2y(n-1)-0.1y(n-2)]/0.5

Verification shows:
y(1)=2
y(2)=3.2
y(3)=4.32
y(4)=5.632

Final Conditions:
y(5)=-3.1168
y(6)=[x(6)-0.2y(5)-0.1y(4)]/0.5=(0-0.2*(-3.1168)-0.1*5.632)/0.5=0.1203

As you can see 0.1203 is not equal to the number -1.1264 which MATLAB 
returns in the
variable s and I really don't understand why...???

Another question is also why stop at n=6. An IIR-filter has an infinite 
response so why
choose n to stop at max(length(a),length(b))-1 as MATLAB does?



Please help. Thank you.

"James" <james@nadaspam.thanx> wrote in message 
news:43ebf95c$0$67258$157c6196@dreader2.cybercity.dk...
> Thanks Fred.  I am not sure if it answers my question though.
> However, I would be very happy if you could tell
> me how matlab calculates the final conditions?

Well you're now asking about the hows and whys of matlab.

I think if you Google on the functions you might learn about it.
I don't see that I can add much for questions like this.

But, a couple of comments:

- read the definition of "zf" in "filter".  Maybe it seems weird but at 
least it seems to be well defined.

- for zf would you expect an array of infinite length?  or, perhaps an array 
of length sufficient to run out of numeric significance?  or, perhaps that 
the word length be increased to continue the sequence for more samples? 
Except the first one, the other two might make sense for somebody.

I really don't know how people normally *use* this function.  Maybe someone 
with experience with matlab and specifically "filter" for IIR might comment.

***!
I didn't think of it before because I didn't pay that much attention but it 
seems pretty cool now that I do think about it:
The zf provide you the initial conditions for the next block!

Something like this:
[Y1,ZF1]=filter(a,b,X1,ZF0)
[Y2,ZF2]=filter(a,b,X2,ZF1)
etc.
where it appears the Xi can be contiguous, nonoverlapping blocks (made 
possible by retaining the initial conditions of the filter).

It appears this helps a whole lot regarding the transients!  It's sure a lot 
simpler. So, I learned something.  Hope you did too.

Fred 

> ***!
> I didn't think of it before because I didn't pay that much attention but 
> it seems pretty cool now that I do think about it:
> The zf provide you the initial conditions for the next block!
>
> Something like this:
> [Y1,ZF1]=filter(a,b,X1,ZF0)
> [Y2,ZF2]=filter(a,b,X2,ZF1)
> etc.
> where it appears the Xi can be contiguous, nonoverlapping blocks (made 
> possible by retaining the initial conditions of the filter).
>
> It appears this helps a whole lot regarding the transients!  It's sure a 
> lot simpler. So, I learned something.  Hope you did too.
>
> Fred

That's exactly why I want to find out how the algorithm works because
I have to implement it in C.

"James" <james@nadaspam.thanx> wrote in message 
news:43ed2bf5$0$78281$157c6196@dreader1.cybercity.dk...
>> ***!
>> I didn't think of it before because I didn't pay that much attention but 
>> it seems pretty cool now that I do think about it:
>> The zf provide you the initial conditions for the next block!
>>
>> Something like this:
>> [Y1,ZF1]=filter(a,b,X1,ZF0)
>> [Y2,ZF2]=filter(a,b,X2,ZF1)
>> etc.
>> where it appears the Xi can be contiguous, nonoverlapping blocks (made 
>> possible by retaining the initial conditions of the filter).
>>
>> It appears this helps a whole lot regarding the transients!  It's sure a 
>> lot simpler. So, I learned something.  Hope you did too.
>>
>> Fred
>
> That's exactly why I want to find out how the algorithm works because
> I have to implement it in C.
>
>
>
Oh, is that all... in that case:

Note what it does.  I think that's all the specification you need:

It computes the output of a filter under two different situations:
- Initial conditions given.
- Initial conditions not given - assumed to be zero.

If the initial conditions aren't given then it computes all of the output 
samples including the startup transient and the decay transient.
The startup transient is included in the standard "output" and the filter 
state in a separate array.
If the initial conditions *are* given, then it computes all of the output 
samples with no startup transient (assuming that the next input has some 
relation to the previous input).
The provision of the initial conditions avoid a startup transient so inputs 
can be in contiguous blocks.

The equations don't change, only the initial conditions.

                               /---\
   x[k]----------+------->-----| + |---->-------+------------> y[k]
                 |      b0     \---/            |
                 |               |              |
                 |               ^              |
                 |               |              |
                 |            +--+--+           |
                 V            |     |           V
                 |            |z^-1 |           |
                 |            +--+--+           |
                 |               |              |
                 |               ^              |
                 |               |              |
                 |             /---\            |
                 +------->-----| + |----<-------+
                 |       b1    \---/   a1       |
                 |               |              |
                 |               ^              |
                 |               |              |
                 |            +--+--+           |
                 V            |     |           V
                 |            |z^-1 |           |
                 |            +--+--+           |
                 |               |              |
                 |               ^              |
                 |               |              |
                 |             /---\            |
                 +------->-----| + |----<-------+
                 |       b2    \---/   a2       |
                 |               |              |
                 V               ^              V


                                 ^
                 V               |              V
                 |            +--+--+           |
                 |            |     |           |
                 |            |z^-1 |           |
                 |            +--+--+           |
                 |               |              |
                 |               ^              |
                 |               |              |
                 |             /---\            |
                 +------->-----| + |----<-------+
                         bn    \---/   an



       Direct Form II Transposed Filter

        y[k]=b0*x[k] + b1*x[k-1] + ... + a1*y[k-1] + a2*y[k-2] - ''''


       H(z) =[b0 + b1*z^-1 +b2*z^-2...]/ [1 - a1*z^-1 - a2*z^-2 ...]



I *think* this is it:

The initial conditions give you the output of each of the delays:
i1 = {b1*x[k-1]+a1*y[k-1]}
i2 = {b2*x[k-2]+a2*y[k-2]}
i3 = {b3*x[k-3]+a3*y[k-3]}

and, looking at the diagram, you compute the new outputs:
y[k]=b0*x[k] + i1
y[k+1]=b0*x[k+1] + b1*x[k] + a1*y[k] + i2
y[k+2]=b0*x[k+2] + b1*x[k+1] + a1*y[k+1] + i3

So it appears to be the same as computing the outputs with zero initial 
conditions and adding the initial conditions to the results.

I would test it!

Fred

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message 
news:j46dnZfl2b87yHDenZ2dnUVZ_tadnZ2d@centurytel.net...

Also:

Because there are no multiplies on the outputs of the delays, the initial 
conditions aren't modified from block to block as the filter coefficients 
change in your time-verying filter.

Fred