Adaptive Filter | ADPCM | ADSP | ADSP-2181 | Aliasing | AMR | Anti-Aliasing | ARMA | Autocorrelation | AutoCovariance | Beamforming | Bessel | Blackfin | Butterworth | C6713 | CCS | Chebyshev | CIC Filter | Circular Convolution | Code Composer Studio | Comb Filter | Compression | Convolution | Cross Correlation | DCT | Decimation | Deconvolution | Demodulation | DM642 | DSP Boards | DSP/BIOS | DTMF | Echo Cancellation | Equalization | Equalizer | ETSI | EZLITE (Ez-kit Lite) | FFT | FFTW | FIR Filter | Fixed Point | FSK | G.711 | G.723 | G.729 | Gaussian Noise | Goertzel | GPIO | Hilbert Transform | IFFT | IIR Filter | Interpolation | Invariance | JTAG | Kalman | Laplace Transform | Levinson | LPC | McBSP | MIPS | Modulation | MPEG | Multirate | Notch Filter | Nyquist | OFDM | Oversampling | Pink Noise | Pitch | PLL | Polyphase | QAM | QDMA | Quantization | Quantizer | Radar | Random Noise | Reed Solomon | Remez | Resampling | RTDX | Sampling | Sharc | TI C6711 | Undersampling | Viterbi | Wavelets | White Noise | Wiener Filter | Windowing | XDS510PP | Z Transform

Hi everyone,

I'm trying to implement a digital filter in c++ (assignment work, so I can't
just use one that someone else has written (--bummer)) and along the way
I've discovered that my knowledge of (low pass) filtering isn't quite what
it should be. I designed my filtering algorithm around convolution of two
sequences, but when testing in Matlab I've come across an issue.

%example code
out=filter([coeffs],[1],[data]); %1 being there simply because the filter
needs both a num & den term
out1=conv(data,coeffs); %forgive if these are the wrong way round; this is
working from memory

Now, my understanding was that filtering was essentially the convolution of
the filter coefficients with the input data. But the fact that out != out1
seems to contradict this. My filter currently gives output agreeing with
conv() but not filter().

Can anyone give me an explanation or hints as to why this may be the case?

Cheers for any input,

XED

P.S I have a copy of Lyons (understanding dsp) and Ifeachor (dsp: a
practical approach) lying around if you can suggest page references.
However, please don't say "it's in ch4, go read it", since that'll take ages
and I may well miss the reasoning. I'm not scared to have a little peek in
the books though..

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

XED wrote:
> Hi everyone,
> 
> I'm trying to implement a digital filter in c++ (assignment work, so I can't
> just use one that someone else has written (--bummer)) and along the way
> I've discovered that my knowledge of (low pass) filtering isn't quite what
> it should be. I designed my filtering algorithm around convolution of two
> sequences, but when testing in Matlab I've come across an issue.
> 
> %example code
> out=filter([coeffs],[1],[data]); %1 being there simply because the filter
> needs both a num & den term
> out1=conv(data,coeffs); %forgive if these are the wrong way round; this is
> working from memory
> 
> Now, my understanding was that filtering was essentially the convolution of
> the filter coefficients with the input data. But the fact that out != out1
> seems to contradict this. My filter currently gives output agreeing with
> conv() but not filter().
> 
> Can anyone give me an explanation or hints as to why this may be the case?

You can filter a signal in either the time-domain or the frequency-domain.

In time-domain filtering, you have an initial transient period. The 
higher the filter order (more "taps"), the longer this transient.

In the frequency-domain using the Fourier transform, you will not see 
this. So if you're doing an ouput sample-by-sample comparison of the two 
techniques, they will *not* be the same. Also, note the frequency-domain 
technique is a "circular convolution".

Good luck,
OUP

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

In article 1...@corp.supernews.com, One Usenet Poster at
m...@my.computer.org wrote on 02/22/2004 19:34:

... 
> You can filter a signal in either the time-domain or the frequency-domain.
> 
> In time-domain filtering, you have an initial transient period. The
> higher the filter order (more "taps"), the longer this transient.
> 
> In the frequency-domain using the Fourier transform, you will not see
> this.

sure you will.  if they are ultimately equivalent, the same behavior has to
be observed in both cases assuming that you're doing it right.  if you are
frequency-domain filtering using overlap-add or overlap-save, you get the
same transient responses to an identical transient input.

> So if you're doing an ouput sample-by-sample comparison of the two
> techniques, they will *not* be the same.

sure they will.  at least to the extent that your bit precision allows.

> Also, note the frequency-domain technique is a "circular convolution".

overlap-add or overlap-save turns circular convolution into linear
convolution.

r b-j

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Thank-you both muchly.

OUP's answer was perfectly sufficient for my purposes, but thanks for your
corrections/additions, RBJ.

Incidentally RBJ, if that surfglobal address is your real email address, I
suggest you ditch it and get a new one. It'll be attracting hordes of spam,
in case you didn't realise.

XED

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

XED wrote:

> Thank-you both muchly.
> 
> OUP's answer was perfectly sufficient for my purposes, but thanks for your
> corrections/additions, RBJ.
> 
> Incidentally RBJ, if that surfglobal address is your real email address, I
> suggest you ditch it and get a new one. It'll be attracting hordes of spam,
> in case you didn't realise.
> 
> XED

For the record, the startup transient in overlap-add and overlap-save 
arises from there being no prior iteration.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Jerry Avins wrote:
> XED wrote:
> 
>> Thank-you both muchly.
>>
>> OUP's answer was perfectly sufficient for my purposes, but thanks for 
>> your
>> corrections/additions, RBJ.
>>
>> Incidentally RBJ, if that surfglobal address is your real email 
>> address, I
>> suggest you ditch it and get a new one. It'll be attracting hordes of 
>> spam,
>> in case you didn't realise.
>>
>> XED
> 
> 
> For the record, the startup transient in overlap-add and overlap-save 
> arises from there being no prior iteration.
> 
> Jerry

FWIW,
You can get rid of the startup transient (or rather shift it to the 
buffer end) by rotating the zero padded filter impulse response to the 
left by nh-1 samples.  This moves the "scrap" samples from the front of 
the output buffer to the back.

Mark Borgerding

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Mark Borgerding wrote:

> Jerry Avins wrote:
> 
>> XED wrote:
>>
>>> Thank-you both muchly.
>>>
>>> OUP's answer was perfectly sufficient for my purposes, but thanks for 
>>> your
>>> corrections/additions, RBJ.
>>>
>>> Incidentally RBJ, if that surfglobal address is your real email 
>>> address, I
>>> suggest you ditch it and get a new one. It'll be attracting hordes of 
>>> spam,
>>> in case you didn't realise.
>>>
>>> XED
>>
>>
>>
>> For the record, the startup transient in overlap-add and overlap-save 
>> arises from there being no prior iteration.
>>
>> Jerry
> 
> 
> FWIW,
> You can get rid of the startup transient (or rather shift it to the 
> buffer end) by rotating the zero padded filter impulse response to the 
> left by nh-1 samples.  This moves the "scrap" samples from the front of 
> the output buffer to the back.
> 
> Mark Borgerding

In real time? There are lots of neat games one can play off line.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

>>
>> FWIW,
>> You can get rid of the startup transient (or rather shift it to the 
>> buffer end) by rotating the zero padded filter impulse response to the 
>> left by nh-1 samples.  This moves the "scrap" samples from the front 
>> of the output buffer to the back.
>>
>> Mark Borgerding
> 
> 
> In real time? There are lots of neat games one can play off line.
> 
> Jerry

I'm not exactly sure what you are asking.

Any fast convolution filtering method introduces a delay due to the 
buffering of input samples.  ( Except for the ludicrous case with one 
fft->ifft operation for each sample. )

The overlap-scrap-tail method does not introduce any delay beyond that 
already imposed by fast-convolution filtering.  Does that answer your 
question?

The main advantage of this is that the final output buffer can be used 
as the output buffer for the ifft -- so long as you can write a few 
extra samples to it.  Those extra samples at the end are the transient 
response and will be overwritten @ the next ifft.

I posted a while ago asking if anyone had heard of this technique 
before.  I got deafening silence ( which proves nothing, of course ).

I will organize my thoughts and put together a pdf explanation.

In the mean time, there are working code examples in kissfft.
http://sourceforge.net/projects/kissfft
see kiss_fastfir.c and tailscrap.m


-- Mark Borgerding


Here's tailscrap.m
============================================
% test code for circular convolution with the scrapped portion
% at the tail of the buffer, rather than the front
%
% The idea is to rotate the zero-padded h (impulse response) buffer
% to the left nh-1 samples, rotating the junk samples as well.
% This could be very handy in avoiding buffer copies during fast filtering.
nh=10;
nfft=256;

h=rand(1,nh);
x=rand(1,nfft);

hpad=[ h(nh) zeros(1,nfft-nh) h(1:nh-1) ];

% baseline comparison
y1 = filter(h,1,x);
y1_notrans = y1(nh:nfft);

% fast convolution
y2 = ifft( fft(hpad) .* fft(x) );
y2_notrans=y2(1:nfft-nh+1);

maxabsdiff = max(abs(y2_notrans - y1_notrans))

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Mark Borgerding wrote:

>>>
>>> FWIW,
>>> You can get rid of the startup transient (or rather shift it to the 
>>> buffer end) by rotating the zero padded filter impulse response to 
>>> the left by nh-1 samples.  This moves the "scrap" samples from the 
>>> front of the output buffer to the back.
>>>
>>> Mark Borgerding
>>
>> In real time? There are lots of neat games one can play off line.
>>
>> Jerry
> 
> I'm not exactly sure what you are asking.
> 
> Any fast convolution filtering method introduces a delay due to the 
> buffering of input samples.  ( Except for the ludicrous case with one 
> fft->ifft operation for each sample. )
> 
> The overlap-scrap-tail method does not introduce any delay beyond that 
> already imposed by fast-convolution filtering.  Does that answer your 
> question?
> 
> The main advantage of this is that the final output buffer can be used 
> as the output buffer for the ifft -- so long as you can write a few 
> extra samples to it.  Those extra samples at the end are the transient 
> response and will be overwritten @ the next ifft.
> 
> I posted a while ago asking if anyone had heard of this technique 
> before.  I got deafening silence ( which proves nothing, of course ).
> 
> I will organize my thoughts and put together a pdf explanation.
> 
> In the mean time, there are working code examples in kissfft.
> http://sourceforge.net/projects/kissfft
> see kiss_fastfir.c and tailscrap.m
> 
> 
> -- Mark Borgerding
> 
> 
> Here's tailscrap.m
> ============================================
> % test code for circular convolution with the scrapped portion
> % at the tail of the buffer, rather than the front
> %
> % The idea is to rotate the zero-padded h (impulse response) buffer
> % to the left nh-1 samples, rotating the junk samples as well.
> % This could be very handy in avoiding buffer copies during fast filtering.
> nh=10;
> nfft=256;
> 
> h=rand(1,nh);
> x=rand(1,nfft);
> 
> hpad=[ h(nh) zeros(1,nfft-nh) h(1:nh-1) ];
> 
> % baseline comparison
> y1 = filter(h,1,x);
> y1_notrans = y1(nh:nfft);
> 
> % fast convolution
> y2 = ifft( fft(hpad) .* fft(x) );
> y2_notrans=y2(1:nfft-nh+1);
> 
> maxabsdiff = max(abs(y2_notrans - y1_notrans))

It looks interesting. The fourth conservation law, Conservation of
Inherent Limitation, makes me skeptical that a shorter start-up
transient than that imposed by a transversal filter is possible. I
hope the limitation isn't really inherent. I eagerly await the PDF.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Hi Jerry,

I had noticed this thread, but had not paid any attention until now.  And
now I can't see the earlier messages.  Well, I could google, but you know...

On Tue, 24 Feb 2004 12:03:12 -0500, Jerry Avins wrote:

> Mark Borgerding wrote:
> 
>>>>
>>>> FWIW,
>>>> You can get rid of the startup transient (or rather shift it to the 
>>>> buffer end) by rotating the zero padded filter impulse response to 
>>>> the left by nh-1 samples.  This moves the "scrap" samples from the 
>>>> front of the output buffer to the back.
>>>>
>>>> Mark Borgerding
>>>
>>> In real time? There are lots of neat games one can play off line.
>>>
>>> Jerry
>> 
>> I'm not exactly sure what you are asking.
>> 
>> Any fast convolution filtering method introduces a delay due to the 
>> buffering of input samples.  ( Except for the ludicrous case with one 
>> fft->ifft operation for each sample. )
>> 
>> The overlap-scrap-tail method does not introduce any delay beyond that 
>> already imposed by fast-convolution filtering.  Does that answer your 
>> question?
>> 
>> The main advantage of this is that the final output buffer can be used 
>> as the output buffer for the ifft -- so long as you can write a few 
>> extra samples to it.  Those extra samples at the end are the transient 
>> response and will be overwritten @ the next ifft.
>> 
>> I posted a while ago asking if anyone had heard of this technique 
>> before.  I got deafening silence ( which proves nothing, of course ).
>> 
>> I will organize my thoughts and put together a pdf explanation.
>> 
>> In the mean time, there are working code examples in kissfft.
>> http://sourceforge.net/projects/kissfft
>> see kiss_fastfir.c and tailscrap.m
>> 
>> 
>> -- Mark Borgerding
>> 
>> 
>> Here's tailscrap.m
>> ============================================
>> % test code for circular convolution with the scrapped portion
>> % at the tail of the buffer, rather than the front
>> %
>> % The idea is to rotate the zero-padded h (impulse response) buffer
>> % to the left nh-1 samples, rotating the junk samples as well.
>> % This could be very handy in avoiding buffer copies during fast filtering.
>> nh=10;
>> nfft=256;
>> 
>> h=rand(1,nh);
>> x=rand(1,nfft);
>> 
>> hpad=[ h(nh) zeros(1,nfft-nh) h(1:nh-1) ];
>> 
>> % baseline comparison
>> y1 = filter(h,1,x);
>> y1_notrans = y1(nh:nfft);
>> 
>> % fast convolution
>> y2 = ifft( fft(hpad) .* fft(x) );
>> y2_notrans=y2(1:nfft-nh+1);
>> 
>> maxabsdiff = max(abs(y2_notrans - y1_notrans))
> 
> It looks interesting. The fourth conservation law, Conservation of
> Inherent Limitation, makes me skeptical that a shorter start-up
> transient than that imposed by a transversal filter is possible. I
> hope the limitation isn't really inherent. I eagerly await the PDF.

By "transversal filter", you mean "Direct-form FIR filter", don't you?  So
what do you mean by "start-up transient"?  Isn't that a characteristic of
the response/coefficients, rather than the implementation?

Since a transversal, or direct-form FIR filter has no inherent limitations
other than causality, and Mark was not suggesting that you could
overlap-save with less delay than the blocking factor, I'm not sure what
the source of your skepticism is?  Could you elaborate?

Mark: how does changing where the "save" part of the result appears in the
ifft buffer "avoid a buffer copy"?  Presumably you have to at least copy
the data into an I/O unit or buffer at some stage, anyway.

Cheers,

-- 
Andrew

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.