# Attempting DSP in Scilab -- filtering using convol()

Started by July 9, 2008
```I'm attempting to use convol() to implement a simple filter.
Scilab's Help entry is terse to say the least.

I suspect I should be using the overlap and add variant, but Scilab
gives no example.

I've experimented and wrote the following code. Calculating the
convolution by both methods. The results can only be described as
"similar". What am I missing?

TIA

t=[0:.00001:1];
s=size(t,2);
h= exp(-t(1:101)/.001);                     // impulse response

f1=123;   / arbitrary test frequencies
f2=5678;
f3=76767;
signal=sin(2*%pi*f1*t) +sin(2*%pi*f2*t) +sin(2*%pi*f3*t); // test signal

rco=convol(h,signal);
sizerco=size(rco);
rco=rco(1:s)/max(rco);

N=int(s/100);
for i=[1:N]
k=i-1;
x(i,:)=signal((1+100*k):(101+100*k));   // I was guessing that h and
x(i)
// should be same size
end

xx=x(1,:);
[yy,ee] =convol(h,xx);         // calc y(1)
y=yy;

for k=[2:N-1]
xx=x(k,:);
[yy,ee] =convol(h,xx,ee);   // calc y(2)->y(N-1)
y=[y,yy];
end

k=k+1;
xx=x(k,:);
yy=convol(h,xx,ee);           // calc y(N)
y=[y,yy];

sizey=size(y)
sizerco

yin=signal(1:s)/max(signal(1:s));
yout=y(1:s)/max(y(1:s));
// plot2d(t(1:s)',[yin(1:s)' yout(1:s)'])
plot2d(t(1:s)',[rco(1:s)' yout(1:s)'])

```
```On Jul 9, 12:52&#4294967295;pm, Richard Owlett <rowl...@atlascomm.net> wrote:
> I'm attempting to use convol() to implement a simple filter.
> Scilab's Help entry is terse to say the least.
>
> I suspect I should be using the overlap and add variant, but Scilab
> gives no example.
>
> I've experimented and wrote the following code. Calculating the
> convolution by both methods. The results can only be described as
> "similar". What am I missing?
>
> TIA
>
> t=[0:.00001:1];
> s=size(t,2);
> h= exp(-t(1:101)/.001); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; // impulse response
>
> f1=123; &#4294967295; / arbitrary test frequencies
> f2=5678;
> f3=76767;
> signal=sin(2*%pi*f1*t) +sin(2*%pi*f2*t) +sin(2*%pi*f3*t); // test signal
>
> rco=convol(h,signal);
> sizerco=size(rco);
> rco=rco(1:s)/max(rco);
>
> N=int(s/100);
> for i=[1:N]
> &#4294967295; &#4294967295;k=i-1;
> &#4294967295; &#4294967295;x(i,:)=signal((1+100*k):(101+100*k)); &#4294967295; // I was guessing that h and
> x(i)
> &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;// should be same size
> end
>
> xx=x(1,:);
> [yy,ee] =convol(h,xx); &#4294967295; &#4294967295; &#4294967295; &#4294967295; // calc y(1)
> y=yy;
>
> for k=[2:N-1]
> &#4294967295; &#4294967295;xx=x(k,:);
> &#4294967295; &#4294967295;[yy,ee] =convol(h,xx,ee); &#4294967295; // calc y(2)->y(N-1)
> &#4294967295; &#4294967295;y=[y,yy];
> end
>
> k=k+1;
> xx=x(k,:);
> yy=convol(h,xx,ee); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; // calc y(N)
> y=[y,yy];
>
> sizey=size(y)
> sizerco
>
> yin=signal(1:s)/max(signal(1:s));
> yout=y(1:s)/max(y(1:s));
> // plot2d(t(1:s)',[yin(1:s)' yout(1:s)'])
> plot2d(t(1:s)',[rco(1:s)' yout(1:s)'])

Richard,

I only looked at the first few lines of what you have. It looks like
you have a sample rate of 100Ksps and one input frequency of
76.767KHz, which is > 1/2 the sample rate, so it will alias, possibly
into the passband of the filter (which you didn't specify).  That
might account for some conflict between your expectations and results.

Dirk.
```
```On Jul 9, 2:40&#4294967295;pm, dbell <bellda2...@cox.net> wrote:
> On Jul 9, 12:52&#4294967295;pm, Richard Owlett <rowl...@atlascomm.net> wrote:
>
>
>
>
>
> > I'm attempting to use convol() to implement a simple filter.
> > Scilab's Help entry is terse to say the least.
>
> > I suspect I should be using the overlap and add variant, but Scilab
> > gives no example.
>
> > I've experimented and wrote the following code. Calculating the
> > convolution by both methods. The results can only be described as
> > "similar". What am I missing?
>
> > TIA
>
> > t=[0:.00001:1];
> > s=size(t,2);
> > h= exp(-t(1:101)/.001); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; // impulse response
>
> > f1=123; &#4294967295; / arbitrary test frequencies
> > f2=5678;
> > f3=76767;
> > signal=sin(2*%pi*f1*t) +sin(2*%pi*f2*t) +sin(2*%pi*f3*t); // test signal
>
> > rco=convol(h,signal);
> > sizerco=size(rco);
> > rco=rco(1:s)/max(rco);
>
> > N=int(s/100);
> > for i=[1:N]
> > &#4294967295; &#4294967295;k=i-1;
> > &#4294967295; &#4294967295;x(i,:)=signal((1+100*k):(101+100*k)); &#4294967295; // I was guessing that h and
> > x(i)
> > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;// should be same size
> > end
>
> > xx=x(1,:);
> > [yy,ee] =convol(h,xx); &#4294967295; &#4294967295; &#4294967295; &#4294967295; // calc y(1)
> > y=yy;
>
> > for k=[2:N-1]
> > &#4294967295; &#4294967295;xx=x(k,:);
> > &#4294967295; &#4294967295;[yy,ee] =convol(h,xx,ee); &#4294967295; // calc y(2)->y(N-1)
> > &#4294967295; &#4294967295;y=[y,yy];
> > end
>
> > k=k+1;
> > xx=x(k,:);
> > yy=convol(h,xx,ee); &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; // calc y(N)
> > y=[y,yy];
>
> > sizey=size(y)
> > sizerco
>
> > yin=signal(1:s)/max(signal(1:s));
> > yout=y(1:s)/max(y(1:s));
> > // plot2d(t(1:s)',[yin(1:s)' yout(1:s)'])
> > plot2d(t(1:s)',[rco(1:s)' yout(1:s)'])
>
> Richard,
>
> I only looked at the first few lines of what you have. It looks like
> you have a sample rate of 100Ksps and one input frequency of
> 76.767KHz, which is > 1/2 the sample rate, so it will alias, possibly
> into the passband of the filter (which you didn't specify). &#4294967295;That
> might account for some conflict between your expectations and results.
>
> Dirk.- Hide quoted text -
>
> - Show quoted text -

Richard,

I see that you actually did specify your filter, so I will correct
myself. Aliasing comment still correct.

Dirk
```
```Richard Owlett wrote:
> I'm attempting to use convol() to implement a simple filter.
> Scilab's Help entry is terse to say the least.
>
> I suspect I should be using the overlap and add variant, but Scilab
> gives no example.
>
> I've experimented and wrote the following code. Calculating the
> convolution by both methods. The results can only be described as
> "similar". What am I missing?
>
> TIA
>
(code snipped)

What are you trying to accomplish?  I tested it out on a few simple
cases; it appears to do a nicely padded-out convolution, so if you want
to do 'straight' convolution you just do the call without the

I really don't know what the overlap-add is about, unless it's to more
efficiently convolve a pair of sequences with widely disparate lengths;
I'd be tempted to just do it all in one step for starters.

You could always examine the code to figure out it's functionality, and
file a bug against the help entry with a suggested wording.

--

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" gives you just what it says.
See details at http://www.wescottdesign.com/actfes/actfes.html
```
```dbell wrote:
> On Jul 9, 2:40 pm, dbell <bellda2...@cox.net> wrote:
>
>>On Jul 9, 12:52 pm, Richard Owlett <rowl...@atlascomm.net> wrote:
>>
>>
>>
>>
>>
>>
>>>I'm attempting to use convol() to implement a simple filter.
>>>Scilab's Help entry is terse to say the least.
>>
>>>I suspect I should be using the overlap and add variant, but Scilab
>>>gives no example.
>>
>>>I've experimented and wrote the following code. Calculating the
>>>convolution by both methods. The results can only be described as
>>>"similar". What am I missing?
>>
>>>TIA
>>[snip code]
>>Richard,
>>
>>I only looked at the first few lines of what you have. It looks like
>>you have a sample rate of 100Ksps and one input frequency of
>>76.767KHz, which is > 1/2 the sample rate, so it will alias, possibly
>>into the passband of the filter (which you didn't specify).  That
>>might account for some conflict between your expectations and results.
>>
>>Dirk.- Hide quoted text -
>>
>>- Show quoted text -
>
>
> Richard,
>
> I see that you actually did specify your filter, so I will correct
> myself. Aliasing comment still correct.
>
> Dirk

Alias folded into "harmless" location. Just gave "another" noise source.
Moved it down to ~36 kHz. Thanks.
```
```Tim Wescott wrote:

> Richard Owlett wrote:
>
>> I'm attempting to use convol() to implement a simple filter.
>> Scilab's Help entry is terse to say the least.
>>
>> I suspect I should be using the overlap and add variant, but Scilab
>> gives no example.
>>
>> I've experimented and wrote the following code. Calculating the
>> convolution by both methods. The results can only be described as
>> "similar". What am I missing?
>>
>> TIA
>>
> (code snipped)
>
> What are you trying to accomplish?

Immediate goal was to figure out how to use Scilab's convol() routine.
Motivation was need to filter a signal with code speed taking a distant
second to code correctness. Not being expert in programming NOR DSP that
indicated using a canned routine.

> I tested it out on a few simple
> cases; it appears to do a nicely padded-out convolution, so if you want
> to do 'straight' convolution you just do the call without the
>
> I really don't know what the overlap-add is about, unless it's to more
> efficiently convolve a pair of sequences with widely disparate lengths;
> I'd be tempted to just do it all in one step for starters.

If it's mathematically correct to convolve sequences with widely
differing lengths, I'll do it.

>
> You could always examine the code to figure out it's functionality, and
> file a bug against the help entry with a suggested wording.
>

First I'd have to figure out DSP ;)

Thanks

```