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..
Convolution versus LPFiltering
Started by ●February 22, 2004
Reply by ●February 22, 20042004-02-22
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
Reply by ●February 23, 20042004-02-23
In article 103iilcfglcrre4@corp.supernews.com, One Usenet Poster at me@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
Reply by ●February 23, 20042004-02-23
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
Reply by ●February 23, 20042004-02-23
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. > > XEDFor 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. �����������������������������������������������������������������������
Reply by ●February 23, 20042004-02-23
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. > > JerryFWIW, 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
Reply by ●February 24, 20042004-02-24
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 BorgerdingIn 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. �����������������������������������������������������������������������
Reply by ●February 24, 20042004-02-24
>> >> 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. > > JerryI'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))
Reply by ●February 24, 20042004-02-24
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. �����������������������������������������������������������������������
Reply by ●February 24, 20042004-02-24
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
