Time domain convolution in a real-time situation

Hi group,

I'm still working on my software BFSK-modem project.
More specifically, i want to use better filters, so I'm experimenting 
with a windowed sinc filter (Blackman window).
I can succesfully calculate the impulse response H[]

My problem is that when I apply the filtering algorithm ( as found in 
ch. 16 of The Scientist & Engineer's Guide to DSP ) , that the first M 
output datapoints in Y[] are zero. ( M = length of H[] )

Quote example code:
  N=4999, M=100

for j := 100 to 4999 do
begin
   Y[j] := 0;
   for i:= 0 to 100 do
       Y[j] := Y[j] + X[j-i] * H[i];
end;

If this were only a transient startup effect, it wouldn't be a problem,
but if I apply the algorithm on each incoming datablock, then I will 
obviously "loose" data.

I understand why the first M points must be zero in a "statical" 
situation, but can't figure out how to handle this in a real-time 
environment, so that I get a continuous filtered output.

I have read about OLA/OLS, but all examples seemed to apply only to FFT 
filters.

Please bear with me, my mathematical background leaves a lot to be 
desired :0(

On Mar 6, 6:01 am, DSP-Newbie <N...@way.invalid> wrote:
> Hi group,
>
> I'm still working on my software BFSK-modem project.
> More specifically, i want to use better filters, so I'm experimenting
> with a windowed sinc filter (Blackman window).
> I can succesfully calculate the impulse response H[]
>
> My problem is that when I apply the filtering algorithm ( as found in
> ch. 16 of The Scientist & Engineer's Guide to DSP ) , that the first M
> output datapoints in Y[] are zero. ( M = length of H[] )
>
> Quote example code:
>   N=4999, M=100
>
> for j := 100 to 4999 do
> begin
>    Y[j] := 0;
>    for i:= 0 to 100 do
>        Y[j] := Y[j] + X[j-i] * H[i];
> end;
>
> If this were only a transient startup effect, it wouldn't be a problem,
> but if I apply the algorithm on each incoming datablock, then I will
> obviously "loose" data.
>
> I understand why the first M points must be zero in a "statical"
> situation, but can't figure out how to handle this in a real-time
> environment, so that I get a continuous filtered output.
>
> I have read about OLA/OLS, but all examples seemed to apply only to FFT
> filters.
>
> Please bear with me, my mathematical background leaves a lot to be
> desired :0(

The comment given at the top of the table that this code came out of
says:

"This program filters 5000 samples with a 101 point windowed-sinc
filter, resulting in 4900 samples of filtered data."

This statement tells you that the author is choosing to overlook the
transients that would occur at the start and end of the output. In
general, if you convolve an M length sequence with an N length
sequence, the output will have length M + N - 1. The author in this
case has ignored the data that would occur on the edges and just
returns the non-transient portion of the output sequence.

The fact that the start of the sequence is ignored is shown by the
fact that the loop index j (which indexes the output array) starts at
100; no values are loaded into the output array Y for indices smaller
than 100. You could certainly change the code to calculate the entire
sequence, starting at index 0, but you would need some additional
logic to handle this case, as in these cases you will not have enough
input samples to evaluate every term of the sum. You would probably
instead assume that the input signal is causal and is therefore zero
for indices less than zero.

Jason

cincydsp@gmail.com wrote:

> You could certainly change the code to calculate the entire
> sequence, starting at index 0, but you would need some additional
> logic to handle this case, as in these cases you will not have enough
> input samples to evaluate every term of the sum. You would probably
> instead assume that the input signal is causal and is therefore zero
> for indices less than zero.

Thanks for the input Jason,

By causal I am assuming I should take the input signal "as is"?

Here's how I do it now:

- I add zeroes (number = length of H[]) at the beginning and end of 
X[].

- I have changed the original algorithm to:

 if j-i < 0 then
   Y[j] := Y[j] + Tx[i] * H[i]
 else
   Y[j] := Y[j] + Tx[j-i] * H[i];

Y[] now becomes something like:
http://users.pandora.be/dirk.claessens2/DSP/bandfilter.JPG

My problem now was, which part of Y[] should I extract and use?
After much and trial & error and some hair-pulling, I finally found a 
intuitive solution that actually works:

As starting index I use 1.5 times the length of H[],
For the extraction length I use the original length of X[].

I've tested it for a wide range of bandwidths, and it works.
I wonder if you or someone else can comment on whether this is:
-  a plain hack ? :-)
-  a mathematically correct method?

Thanks - Dirk








   // 75=900 125=500
   Y := Copy(y, Round(Length(H)*1.5) , Length(X));

On 6 Mar, 12:01, DSP-Newbie <N...@way.invalid> wrote:

--- start up effects snipped ---

> If this were only a transient startup effect, it wouldn't be a problem,
> but if I apply the algorithm on each incoming datablock, then I will
> obviously "loose" data.
>
> I understand why the first M points must be zero in a "statical"
> situation, but can't figure out how to handle this in a real-time
> environment, so that I get a continuous filtered output.

The issue with missing data points only appears at the very start
of the data sequence. Once the filter starts running, it always
contains valid data, you don't clear the buffer between frames.

Rune

> By causal I am assuming I should take the input signal "as is"?

By causal, I meant that the signal is zero for all instances of time
before time zero.

> Here's how I do it now:
>
> - I add zeroes (number = length of H[]) at the beginning and end of
> X[].
>
> - I have changed the original algorithm to:
>
>  if j-i < 0 then
>    Y[j] := Y[j] + Tx[i] * H[i]
>  else
>    Y[j] := Y[j] + Tx[j-i] * H[i];
>
> Y[] now becomes something like:http://users.pandora.be/dirk.claessens2/DSP/bandfilter.JPG

Explicitly padding the input signal with zeros like this will work.
However, note that you're essentially delaying the input signal by
length(H) samples, so your output signal (even before you see any
transients) will also be delayed by length(H) samples.

> My problem now was, which part of Y[] should I extract and use?
> After much and trial & error and some hair-pulling, I finally found a
> intuitive solution that actually works:

-snipped-

I forgot to ask before: why do you want to discard the transients? The
filter "startup" time is valid data, and occurs as your input signal
gradually fills up the delay lines that feed each filter tap. If you
don't want to see any of the startup transients, then you would just
have to discard all of the samples before all of the filter's taps
have been filled. This doesn't occur until length(H) input samples
have entered the filter. You'll also see the same effect on the way
out; after your input signal ends (and you start passing the zeros
that you padded onto the end into the filter), you'll have length(H)
instances of time where some of the input data stored in the filter's
delay lines are the zeros that you used for padding. If you're doing
the filtering on a block basis and you're going to continue to feed
the filter with new data as time goes on, then you need to maintain
these transients to get correct results.

Remember, all you're trying to do here is implement the convolution
sum, nothing more. The process of convolution causes the transients to
appear, and they are a valid part of its output.

Jason

cincydsp@gmail.com wrote:

>
> Explicitly padding the input signal with zeros like this will work.
> However, note that you're essentially delaying the input signal by
> length(H) samples, so your output signal (even before you see any
> transients) will also be delayed by length(H) samples.
>

I am aware of that, but its is not a problem; the delay is the same for 
all samples.



>> My problem now was, which part of Y[] should I extract and use?
>> After much and trial & error and some hair-pulling, I finally found a
>> intuitive solution that actually works:
>
> -snipped-
>
> I forgot to ask before: why do you want to discard the transients? The
> filter "startup" time is valid data, and occurs as your input signal
> gradually fills up the delay lines that feed each filter tap. 

Perhaps I have not explained my problem clearly.
I have set the soundcard to sample at 11.025 KHz, and to report the 
samples in blocks of 2048 samples, so I have to process about 5 
samples/sec.

Each of these samples - after filtering, envelope detection, PLL bit 
sampler etc... - represents a number of bits, not necessarily a *whole* 
number of bits, which I have to cobble together block by block to 
reconstruct the original synchronous bitstream.

Now, if I apply the algorithm as found in the book for each incoming 
datablock, I will have a startup transient for *each* block. There is 
no history from previous blocks. This results in repeating bit errors.

But, as I mentioned in my previous reply, the zero-padding etc. now 
made it working correctly.
BTW:  the filter performs **tremendously** well; I can now decode 
almost inaudible signals totally buried in the noise.


[snip]
> Remember, all you're trying to do here is implement the convolution
> sum, nothing more. 

Sure, but I've had a hard time with it >:|

Thanks again for the explanations Jason, it is appreciated!

Rune Allnor wrote:
> On 6 Mar, 12:01, DSP-Newbie <N...@way.invalid> wrote:
>
> --- start up effects snipped ---
>
>> If this were only a transient startup effect, it wouldn't be a problem,
>> but if I apply the algorithm on each incoming datablock, then I will
>> obviously "loose" data.
>> 
>> I understand why the first M points must be zero in a "statical"
>> situation, but can't figure out how to handle this in a real-time
>> environment, so that I get a continuous filtered output.
>
> The issue with missing data points only appears at the very start
> of the data sequence. Once the filter starts running, it always
> contains valid data, you don't clear the buffer between frames.
>
> Rune

Thanks for the reaction Rune, but see my reply to Jason.

DSP-Newbie wrote:
> cincydsp@gmail.com wrote:
> 
>>
>> Explicitly padding the input signal with zeros like this will work.
>> However, note that you're essentially delaying the input signal by
>> length(H) samples, so your output signal (even before you see any
>> transients) will also be delayed by length(H) samples.
>>
> 
> I am aware of that, but its is not a problem; the delay is the same for 
> all samples.
> 
> 
> 
>>> My problem now was, which part of Y[] should I extract and use?
>>> After much and trial & error and some hair-pulling, I finally found a
>>> intuitive solution that actually works:
>>
>> -snipped-
>>
>> I forgot to ask before: why do you want to discard the transients? The
>> filter "startup" time is valid data, and occurs as your input signal
>> gradually fills up the delay lines that feed each filter tap. 
> 
> Perhaps I have not explained my problem clearly.
> I have set the soundcard to sample at 11.025 KHz, and to report the 
> samples in blocks of 2048 samples, so I have to process about 5 
> samples/sec.
> 
> Each of these samples - after filtering, envelope detection, PLL bit 
> sampler etc... - represents a number of bits, not necessarily a *whole* 
> number of bits, which I have to cobble together block by block to 
> reconstruct the original synchronous bitstream.
> 
> Now, if I apply the algorithm as found in the book for each incoming 
> datablock, I will have a startup transient for *each* block. There is no 
> history from previous blocks. This results in repeating bit errors.

Then it is clearly the wrong algorithm. Look for overlap-add in the index.

> But, as I mentioned in my previous reply, the zero-padding etc. now made 
> it working correctly.
> BTW:  the filter performs **tremendously** well; I can now decode almost 
> inaudible signals totally buried in the noise.
> 
> 
> [snip]
>> Remember, all you're trying to do here is implement the convolution
>> sum, nothing more. 
> 
> Sure, but I've had a hard time with it >:|

You probably ended up implementing one of the overlap methods.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

"DSP-Newbie" <No@way.invalid> wrote in message 
news:mn.35547d7371047ae5.59994@way.invalid...
> Rune Allnor wrote:
>> On 6 Mar, 12:01, DSP-Newbie <N...@way.invalid> wrote:
>>
>> --- start up effects snipped ---
>>
>>> If this were only a transient startup effect, it wouldn't be a problem,
>>> but if I apply the algorithm on each incoming datablock, then I will
>>> obviously "loose" data.
>>>
>>> I understand why the first M points must be zero in a "statical"
>>> situation, but can't figure out how to handle this in a real-time
>>> environment, so that I get a continuous filtered output.
>>
>> The issue with missing data points only appears at the very start
>> of the data sequence. Once the filter starts running, it always
>> contains valid data, you don't clear the buffer between frames.
>>
>> Rune
>
> Thanks for the reaction Rune, but see my reply to Jason.
>

Nowhere did you seem to pick up on Rune's suggestion that you *not* throw 
out the filter contents when moving to the next block....

Fred 

Fred Marshall wrote:

   ...

> Nowhere did you seem to pick up on Rune's suggestion that you *not* throw 
> out the filter contents when moving to the next block....

It seems like a sad case of canned code that he doesn't fully 
understand, but as long as it seems to work, what the hey! Sort of how I 
write my web pages. I hope I'm wrong.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯