Farrow Filter implementation in Matlab

Hi,

for my diploma thesis I am studying different resampling techniques.
Just now I wrote a Farrow resampler in Matlab. The structure of the
resampler is quite easy to implement, but I struggle with the
polynomial coefficients. I tried to derive them from FIR filters, but
although the filter approximations look good in time domain, the
attenuation is quite low (20dB-40dB for an original filter designed
with nearly 80dB attenuation).

I currently use 10 polynoms of 4th order with a sufficiently (or what
I considered sufficiently) shaped filter. I tried different filter
shapes (equiripple (doesn't work so well), different windows, least-
squares).

Does anyone have recommendations for the practical filter
coefficients? Can you guess what I did wrong?

> Can you guess what I did wrong?

For example, trying to use a filter designed using a least-squares method
(such as Parks McClellan)?

Check the impulse response of your prototype, whether the outermost taps
fall on a "smooth" impulse response with the others.
Some methods give filters where the outermost coefficients are
"discontinuous", and that can't be approximated with a polynomial.

That's a guess, as you asked, so it may well be something completely
different.
-mn

Hi,

Thanks for the answer.
I think I figured it out now...

I also tried out Kaiser and other windows to design a filter -- and
then I decided to use the example coefficients from a paper which give
me 50dB attenuation. And I use the filter in a place in my
architecture, where this is enough... a filter guy from our department
said, that 50dB ar probably quite good already, because the
coefficients are only an approximation and small errors have big
effects in high attenuation.

-- Imke

>>and small errors have big
>effects in high attenuation.

yes, small errors have big effects on small numbers :)

You could compare the performance of your Farrows implementation with that
of the filter prototype. Or even easier: if the prototype alone doesn't
meet the specs you know where you have to look for.
Otherwise, it's the polynomials.

> Otherwise, it's the polynomials.

I found a suitable place in my processing architecture for the farrow
-- if I put it in front of a shaper FIR, that one can do the rest of
the attenuation.

not sure if I posted this reference to one of your posts.
http://yehar.com/dsp/deip.pdf

On Nov 16, 8:29 am, Imke Huismann <pikku.kor...@gmail.com> wrote:
> Hi,
>
> for my diploma thesis I am studying different resampling techniques.
> Just now I wrote a Farrow resampler in Matlab. The structure of the
> resampler is quite easy to implement, but I struggle with the
> polynomial coefficients. I tried to derive them from FIR filters, but
> although the filter approximations look good in time domain, the
> attenuation is quite low (20dB-40dB for an original filter designed
> with nearly 80dB attenuation).
>
> I currently use 10 polynoms of 4th order with a sufficiently (or what
> I considered sufficiently) shaped filter. I tried different filter
> shapes (equiripple (doesn't work so well), different windows, least-
> squares).
>
> Does anyone have recommendations for the practical filter
> coefficients? Can you guess what I did wrong?

You may want to take a look at Harris' book "Multirate Signal
Processing for Communication Systems" pgs 190-191 which discuss the
Farrow Filter. The following is essentially paraphrasing the book with
the numbers changed for your application:

Assuming out of band sidelobe levels are bounded by 5dB /bit. If you
need 80 dB attenuation then this is 16 bits. Then the precision needed
is 2/(2^16) ~= 3.05e-5. Thus the error in your polynomial
approximations need to be less than this.

He also warns about the Remez filters which have constant sidelobe
levels. You get an 6dB /octave falloff in the sidelobe's for every
continuous n'th derivative (this is an asymptotic value). Obviously
the usual Parks McCellan(Remez) filters don't have this, so there is
an impulse at the ends of the time domain filter, which subsequently
precludes using low order polynomials near this area.

Hope that helps.
Cheers,
David

On Nov 28, 6:02 pm, Dave <dspg...@netscape.net> wrote:
> On Nov 16, 8:29 am, Imke Huismann <pikku.kor...@gmail.com> wrote:
>
>
>
> > Hi,
>
> > for my diploma thesis I am studying different resampling techniques.
> > Just now I wrote aFarrowresampler in Matlab. The structure of the
> > resampler is quite easy to implement, but I struggle with the
> > polynomial coefficients. I tried to derive them from FIR filters, but
> > although thefilterapproximations look good in time domain, the
> > attenuation is quite low (20dB-40dB for an originalfilterdesigned
> > with nearly 80dB attenuation).
>
> > I currently use 10 polynoms of 4th order with a sufficiently (or what
> > I considered sufficiently) shapedfilter. I tried differentfilter
> > shapes (equiripple (doesn't work so well), different windows, least-
> > squares).
>
> > Does anyone have recommendations for the practicalfilter
> > coefficients? Can you guess what I did wrong?
>
> You may want to take a look at Harris' book "Multirate Signal
> Processing for Communication Systems" pgs 190-191 which discuss theFarrowFilter. The following is essentially paraphrasing the book with
> the numbers changed for your application:
>
> Assuming out of band sidelobe levels are bounded by 5dB /bit. If you
> need 80 dB attenuation then this is 16 bits. Then the precision needed
> is 2/(2^16) ~= 3.05e-5. Thus the error in your polynomial
> approximations need to be less than this.
>
> He also warns about the Remez filters which have constant sidelobe
> levels. You get an 6dB /octave falloff in the sidelobe's for every
> continuous n'th derivative (this is an asymptotic value). Obviously
> the usual Parks McCellan(Remez) filters don't have this, so there is
> an impulse at the ends of the time domainfilter, which subsequently
> precludes using low order polynomials near this area.
>
> Hope that helps.
> Cheers,
> David

Hi,

I read this book (or at least the resampling chapter) and also some of
Harris' papers. I tried several window types and also equiripple
(remez) filters as prototype filters -- but I don't believe that it is
possible to get much more than 50dB attenuation.

In the paper posted by mnentwig (http://yehar.com/dsp/deip.pdf) there
are several polynomial filtering methods reviewed and non of them got
more than about 50 dB. Although the pictures in  the paper
"Performance and design considerations of the Farrow filter when used
for arbitrary resampling of sampled time series" from Harris, F. look
more promising (60 dB). I used the coefficients and got 50dB.

I think it is not the right approach to try to get 80dB with a Farrow
resampler -- if it is possible it is probably a very cost intensive
implementation. Please correct me when I'm wrong.

Cheers
 Imke

On Nov 29, 5:28 am, Imke Huismann <pikku.kor...@gmail.com> wrote:
> On Nov 28, 6:02 pm, Dave <dspg...@netscape.net> wrote:
>
>
>
> > On Nov 16, 8:29 am, Imke Huismann <pikku.kor...@gmail.com> wrote:
>
> > > Hi,
>
> > > for my diploma thesis I am studying different resampling techniques.
> > > Just now I wrote aFarrowresampler in Matlab. The structure of the
> > > resampler is quite easy to implement, but I struggle with the
> > > polynomial coefficients. I tried to derive them from FIR filters, but
> > > although thefilterapproximations look good in time domain, the
> > > attenuation is quite low (20dB-40dB for an originalfilterdesigned
> > > with nearly 80dB attenuation).
>
> > > I currently use 10 polynoms of 4th order with a sufficiently (or what
> > > I considered sufficiently) shapedfilter. I tried differentfilter
> > > shapes (equiripple (doesn't work so well), different windows, least-
> > > squares).
>
> > > Does anyone have recommendations for the practicalfilter
> > > coefficients? Can you guess what I did wrong?
>
> > You may want to take a look at Harris' book "Multirate Signal
> > Processing for Communication Systems" pgs 190-191 which discuss theFarrowFilter. The following is essentially paraphrasing the book with
> > the numbers changed for your application:
>
> > Assuming out of band sidelobe levels are bounded by 5dB /bit. If you
> > need 80 dB attenuation then this is 16 bits. Then the precision needed
> > is 2/(2^16) ~= 3.05e-5. Thus the error in your polynomial
> > approximations need to be less than this.
>
> > He also warns about the Remez filters which have constant sidelobe
> > levels. You get an 6dB /octave falloff in the sidelobe's for every
> > continuous n'th derivative (this is an asymptotic value). Obviously
> > the usual Parks McCellan(Remez) filters don't have this, so there is
> > an impulse at the ends of the time domainfilter, which subsequently
> > precludes using low order polynomials near this area.
>
> > Hope that helps.
> > Cheers,
> > David
>
> Hi,
>
> I read this book (or at least the resampling chapter) and also some of
> Harris' papers. I tried several window types and also equiripple
> (remez) filters as prototype filters -- but I don't believe that it is
> possible to get much more than 50dB attenuation.
>
> In the paper posted by mnentwig (http://yehar.com/dsp/deip.pdf) there
> are several polynomial filtering methods reviewed and non of them got
> more than about 50 dB. Although the pictures in  the paper
> "Performance and design considerations of the Farrow filter when used
> for arbitrary resampling of sampled time series" from Harris, F. look
> more promising (60 dB). I used the coefficients and got 50dB.
>
> I think it is not the right approach to try to get 80dB with a Farrow
> resampler -- if it is possible it is probably a very cost intensive
> implementation. Please correct me when I'm wrong.
>
> Cheers
>  Imke

From Harris' description he uses a least squares polynomial fit, which
will be dependent on the prototype filter coefficients used. You also
have some control over the decomposition of the filter used in the
Farrow. The polynomials are calculated off line and fixed - they only
need to evaluated at different values. The computational cost is only
upfront and remains fixed during runtime.

This is fundamentally different than using fixed polynomials which try
to fit the input of the filter. The fit and resulting error will be
dependent on the actual input data - this is different  than in the
Farrow, where the fit is predetermined by the filter. I don't think
the results in the paper posted by mnentwig are applicable to the
results of the Farrow. I believe there is a paper by some scientists
in FMV (Swedish Defence Research) which essentially says that it is
best to upsample and filter first before applying polynomial
approximations for data fits. Even the paper posted hints at that.
Consider if the -60 dB stop band already sits on top of sidelobes 80
dB down - you still end up with sidelobes that are atleast 80 dB down
- so the polynomial filter isn't what is driving the performance.

That said, I'll qualify - this isn't my exact area of expertise so
what I'm saying could be wrong. Also, your result could be correct -
it may be difficult to achieve the 80 dB attenuation via a Farrow
implementation, the goal of my post was simply to provide the original
poster with some quantitative measure of what was needed in order to
achieve the desired attenuation levels - whether that is actually
achievable is another question.

Cheers,
David

I just wrote some matlab code to do a polynomial interpolation of the
filter coefficients. To make it work properly I had to handle the end
of each subfilter properly - essentially this means including the next
tap in the filter in the polynomial fit. In my case the original
filter had 72 dB attenuation. In the evaluations the resulting
interpolated filter does meet these requirements. See the line
b2=[b2;b2(1,2:end),0]; .

Note: I had to use a 5th order polynomial to do the fits.


Here's my matlab code:

% b = coefficients of a predesigned filter
% In my test case b has 300 taps
% Designed using 72 dB attenuation Kaiser window filter
Ntaps=length(b);
D=50;
Npoly=Ntaps/D;
b2= reshape(b,D,Npoly);
%For a good polynomial fit we need to handle the end points in a
smooth
%fashion
b2=[b2;b2(1,2:end),0];
x=(0:D)';
p_ord = 5; %Order of polynomials for fitting
p=zeros(Npoly,p_ord+1);
for k =1:Npoly
    k
    p(k,:) =polyfit(x,b2(:,k),p_ord);
    %This plots the error in the polynomial evaluation
    %plot(b2(:,k)-polyval(p(k,:),x))
    %pause
end

%Lets evaluate the resulting filter at some other sample point.
x2 = x(1:(length(x)-1))+0.5;
h2=[];
for k =1:Npoly
    h2=[h2,polyval(p(k,:),x2)];
end
plot(20*log10(abs(fft(h2(:),2048))))

Hope that helps.
Cheers,
David