FPGA Central

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

[email protected] 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.claesse...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.claesse...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

[email protected] 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:
> [email protected] 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" <[email protected]> wrote in message
news:[email protected]..
> 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.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

"DSP-Newbie" <[email protected]> wrote in message
news:[email protected]..
> [email protected] 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's a problem you don't need to carry along. Don't do what the book
did. Do continuous processing - if even the continous process has to wait
between availability of data. Something like this:

Get 2048 data points:
Compute 2048 output points. Initially this will include the M-1 startup
transient points.
Wait for another 2048 block.
Continue where you left off using M-1 values of the last input block still
in the filter and the one next new input point.
There will be no transients in the second and succeeding blocks.

The code won't look all that much different from what you're using.

Take a look at the circular code at:
http://www.dspguru.com/sw/lib/fir_algs_1-0.c

Fred

On Mar 6, 10:05 pm, "Fred Marshall" <fmarshallx@remove_the_x.acm.org>
wrote:
> "DSP-Newbie" <N...@way.invalid> wrote in message
>
> news:[email protected]..
>
>
>
> > cincy...@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's a problem you don't need to carry along. Don't do what the book
> did. Do continuous processing - if even the continous process has to wait
> between availability of data. Something like this:
>
> Get 2048 data points:
> Compute 2048 output points. Initially this will include the M-1 startup
> transient points.
> Wait for another 2048 block.
> Continue where you left off using M-1 values of the last input block still
> in the filter and the one next new input point.
> There will be no transients in the second and succeeding blocks.
>
> The code won't look all that much different from what you're using.
>
> Take a look at the circular code at:http://www.dspguru.com/sw/lib/fir_algs_1-0.c
>
> Fred

I'm amazed that it's even working at all. Or maybe he really is
keeping the filter history and doesn't realize it. I can't see how you
would get anything meaningful if the filter history was flushed after
every block; you'd have to be really lucky, I guess.

Jason

Fred Marshall wrote:
> The code won't look all that much different from what you're using.
>
> Take a look at the circular code at:
> http://www.dspguru.com/sw/lib/fir_algs_1-0.c
>
> Fred

Many thanks Fred, this is *exactly* what I needed to implement
continuous block-to-block filtering!

printf("Testing fir_circular:\n ");
clear(NTAPS, z);
state = 0;
for (ii = 0; ii < IMP_SIZE; ii++) {
output = fir_circular(imp[ii], NTAPS, h, z, &state);
printf("%3.1lf ", (double) output);
}

Regards - Dirk

DSP-Newbie wrote:
> Fred Marshall wrote:
>> The code won't look all that much different from what you're using.
>>
>> Take a look at the circular code at:
>> http://www.dspguru.com/sw/lib/fir_algs_1-0.c
>>
>> Fred


And I finally got it working...
Thanks to all responders!

<http://users.pandora.be/dirk.claessens2/DSP/scrshot.JPG>

DSP-Newbie <[email protected]> writes:

> DSP-Newbie wrote:
>> Fred Marshall wrote:
>>> The code won't look all that much different from what you're using.
>>>
>>> Take a look at the circular code at:
>>> http://www.dspguru.com/sw/lib/fir_algs_1-0.c
>>>
>>> Fred
>
>
> And I finally got it working...
> Thanks to all responders!
>
> <http://users.pandora.be/dirk.claessens2/DSP/scrshot.JPG>

Very Cool!
--
% Randy Yates % "The dreamer, the unwoken fool -
%% Fuquay-Varina, NC % in dreams, no pain will kiss the brow..."
%%% 919-577-9882 %
%%%% <[email protected]> % 'Eldorado Overture', *Eldorado*, ELO
http://home.earthlink.net/~yatescr