Yeah, FIR filter design basically. There's a little bit of math in this talk. I have graphs
for mostly everything. So if you want to you can ignore the formulas and you can just
look at the graphs and also the last half of the talk is practical examples and applications.
So if you become lost or bored, that is a good point to rejoin. So what is there about
quick show of hands? Who knows about the parks, McLean, also known as Remes, filter design algorithm?
Okay, a few people. Next question, who does peep install PM Remes or you can substitute
peep install by whatever is your favorite way of installing pattern packages?
Marcus, maybe? Yeah. Okay, so the goal of these talk is to get more people racing hands for
both of these questions, especially the first one because the algorithm is good. The other one
is, I think, a good way of using it in Python, but you can use your favorite implementation.
So a quick recap about FIR design. If we want to design a Lopos filter, the filter is supposed
to be one which means let us in and through in the passband and zero exactly zero in the stopband,
which means kilo of the signal here. And for sure that's mathematically impossible. So the
best thing we can do is something like this, which is a practical example of a 64 tap FIR filter.
We have 64 numbers and this is the frequency response we get. This is the design with a humming
window using the fairway method. So this is something you might have done already. The two main
differences between these filters is a transition window because we cannot have this tip jump
from 1 to 0. And then also if we look at the passband and the stopband in more detail,
you can zoom in. And then in here we can see that there is a repo. So it's not exactly one.
It waves with an error. And in the stopband we can see it's normally plotted in decibels.
So we can see it reels about zero. And that means we are not canceling signals. We are
attenuating them by about 60 degrees. Okay, can we do better than something like this?
That is the question. So let's say I have a fixed number of taps because computation increases
with a number of taps. So let's say with the computation I can afford with 64 taps,
can I do any better than this window method? And the idea is, so I'm always going to have an error,
but what I want to do is with this fixed number of taps, I want to find the affair filter, which
has the lowest maximum error. What is this maximum error? So for example, in here for the blue one,
which is the window method, the maximum error is this one. Here's better, it has lower error.
But if the thing we want to optimize is the maximum error, the algorithm we need is these
parts, my clan, and algorithm. That's basically the goal of the algorithm. And we can see that
in orange, we basically get a much better filter with exactly the same number of taps. So that's the
fair thing I want to make clear in the talk. You should be using this. For most FIA design applications,
it's the right algorithm to use. You just get much better filters. As you can see in here,
there are no tricks. It's just the right number of taps. So, okay, this is kind of the mathy way of
stating what the algorithm does, it's basically what I said. It's finally the affair filter with
a fixed number of taps that minimizes the maximum error. And you can also put a weight function
for the error, which is going to be useful in the practical applications. We will see a couple examples.
Technical thingy, this is only supposed to work for real filters. It actually needs to have real
taps and the impulse response. So that's the taps needs to be either symmetric or anti-symmetric,
basically even symmetric or alt-symmetric. That is equivalent to the filter having linear phase,
which is the same as having constant group delay. So unless you are doing some dispersion
information or dispersion cancellation, these are the types of filters you are going to work with.
If you need to do some dispersion, then this is not the right algorithm. It just doesn't work.
So let's take a detour on polynomial approximation. Let's say I have this Gaussian function in blue.
And I want to find the polynomial, which approximates this function. For simplicity,
I want to find a degree five polynomial, just because if I had a higher degree,
the pot would get messier. A degree five polynomial has six coefficients. So with these six numbers,
I want to approximate these Gaussian function, which is an example. It could be whatever
continuous function we want to have here. And there is this theorem by ChevyChef,
which says that my polynomial is the best uniform approximation to this function. If there is this
equi oscillation property, and I'm going to explain what that means with a diagram,
because it sees it on the standard math, the math you have it up there. So to come up with
this polynomial, what I've done is I'm choosing interpolation notes in green. These are the six
green points. And then I compute the unique polynomial, which interpolates the function. It basically
has the same value as the function in the green points. And you can see that the two
blue and orange plots, they cross each other on the green points. That's how we come up with the
polynomial. And then we look at the local extreme for the error. So that's here, here,
here, these red points. That's when the orange and blue traces, they are like locally further apart.
Right, it's the maximum for them. A local maximum or local minimum for the error.
And what happens here, you might notice here the error is pretty large, here the error is pretty small,
right? So one, my thing, this is like intuitive, maybe I can wiggle a bit the coefficients of my
polynomial to try to make this thing go up and so the error decreases. And because nothing is free
in life when you are trying to optimize for something, then I would have to pay with something else,
which is probably the error here grows. So if I try to push this red thing up, these are the red
thing would go down. Okay, that's fair. But I mean, if I want to minimize the maximum error,
that's the right thing to do, right? I want to basically push this up. These will automatically
push down until both of them have this error. So that's kind of the idea of this theorem.
And then the question is, I'm mentioning I need to push these things, so how do I actually do it?
So the thing is basically we start with points on this interval in which we want to do the
approximation. They can be any points, you can take them equally spaced or or anything.
We solve this linear system. And here the thing is, so this thing is supers to be the error.
And then we have this plus and minus one because the error is supposed to alternate. So if I go back,
you can notice that because polynomials are like weekly things, they do not want to follow a Gaussian curve
or anything else. I am like above below and so on and so forth. So if I want these distances to
be all the same, what I want to do is basically something like this. I want to find these B
coefficients of a polynomial so that they go up and down the function, which is G here,
with respect to some error, which is E. So that is the linear system. I need to solve
there are M plus two equations, there are M plus two points. So there's a single solution. We solve
it. We have our polynomial, but then the thing will happen is I want this thing to be the error
extrema. But in general they are not going to be. So we need to find the local extrema of the error.
And then we replace the n plus two points that we started with by this local extrema.
And then we basically iterate. Okay, this is best explained with just the plots of the thing.
So I started with the thing I had before. And as I said, I want to push this thing up.
I want to push this thing down. So what I do is I compute these function,
which in the green nodes, it's supposed to go like up and down, up and down. And it's probably not
clear from this picture, but these distances on the green points,
there are exactly the same, because that's the linear program I'm solving for.
But then what happens is if I look at the error extrema, they are not exactly at the green points.
You can see them in red again. And for example, this one is actually at a green point.
But this one is slightly different and the error is slightly larger in here.
Okay, that's fine. So what we do now is we basically iterate. We put green points where we had red points.
We repeat the same thing. And now the next thing we get is basically the green points on the red points.
Now they are very, very similar. You cannot tell them apart in this plot. And so we can see now
all of these distances, they are the same. So that's the thing we wanted to have. We wanted to have
a five degree polynomial, which approximates best discussion curve in the uniform sense.
So this is kind of the recipe. But now you might be saying, hey, this is the DSP room.
I've come here to learn about FR filters. What does this thing have to do with FR filters at all?
Unfortunately, I don't have like a picture for this, but I want to convince you really quickly
that basically computing an FR filter and doing this polynomial approximation,
thinking they are kind of equivalent and because of some easy math. So I'm going to do it for the
old length FR with even symmetry. That's the easiest one, but the other ones are similar.
And then the thing is, so what is the frequency response of the FR filter? It is basically this thing.
It's a Fourier transform of the coefficients. And here there is some symmetry because I'm
saying that, for example, A0 is the same as A2n, right? The first coefficient is the same as the
last and so on and so forth. So what I can do to exploit this symmetry is to take out this thing as
a common factor. And surely I need to put it in there. I'm also writing my summation index
so that the symmetry between the coefficients is more obvious. So happens is the minus n is
emitted with the plus n and so on and so forth. And then what happens? For example, with the first one
and the last one, because they now have the same exponential, but with an opposite sign,
I can bring them together. I can use these allus formula and the thing I get is a cosine function.
So basically I get these B0 times the cosine of something, blah, blah. So I get a bunch of cosine
and these BK coefficients, which are basically two times the A's, except for the one in the
middle where you just have that one because you don't have two things to put together. So next thing,
we can do a change of variables to this cosine, we call it X. And then what happens is this function
here, can be written as a polynomial of degree n in X. And the thing you need to remember is
here actually, I don't have like powers of X. I don't have powers of cosine. I have cosine of
something with an n or with a k inside, but we can write that as powers of cosine function. For example,
for k equals two, so that's the cosine of twice the angle, you might remember this formula.
And then because I'm doing, so now I have powers of cosine. And now I can make the change,
like cosine, blah, blah, blah. And you can do this for free k. And so that's basically what happens.
Any a fire filter like this, you can write the frequency response as a polynomial by using this
change of variables. And so you can use techniques for polynomial approximation to design a fire
filters. Okay, so open source implementations of these things. Here I should first say that
the first implementation of these comes from parts, McLean and maybe Robiner. I don't remember
they are like a couple papers from the 70s. And in one of them it's just parts and McLean,
the other one, there's also Robiner. In one of those papers, they give you the fourth one code
for these. It's like the paper is description of the algorithm and then annex fourth one code.
And I don't know about the license in the 70s, people were not so worried about licenses of these things.
And in any, any how it's a mathematical algorithm, read and down in fourth time.
So what happens to these days, if you are using sci-fi, the ramus algorithm which is in there
is the C translation of these original four ton code. And that means it's horrible. You can look at it
like the names of the variables are the names which were reasonable to put on punch cards for code
and on staff. So okay, horrible and readable code. Then you have, so when you want to choose
your filter design, how my response function is going to look like, how my width function is
going to look like, these are piecewise constant. So for example, for the low pass filter
I gave, you have one in the pass bank, you have zero in the stop bank. And that is fine for basic
filters, but we will see some example use cases where we want to have custom functions for those.
So can you do that with sci-fi? No. And what is the main reason for that? Well,
in the 70s paper, the authors were like, yes, we are shipping you some, I think it's called routines
in 490, it doesn't even have like functions. So we are giving you some very simple routines to
define the constant response and weights. In case you need and they were aware that for some
use cases you need like to put in here a function which is not constant, just modify the code,
which is a reasonable thing to do in 490 and 70s. But these days, you know, we have like
closures, lambda functions and stuff in programming languages. And so we should be able to put
in there any function we want and we will see examples. Then license in BSD license,
and all the things I'm going to explain in more detail in a couple slides,
sometimes the algorithm doesn't converge. So you might be putting a program for the algorithm
and it just fails in whatever surprise and way. Then we have going to radio our octave.
And there is an ore here because basically the going to radio implementation is a copy and paste
from octave and it hasn't been changed in per year decade. More two decades. Okay, so
it is seen written from scratch in the 90s. They are coming saying there are things which are wrong.
It is like author of the code says this routine doesn't make any sense, I should fix it and
the comment is from the 90s and no one has fixed it. Then the same piecewise constant response
on weights, it is GPL 2.0 later. That's not only because of the radio, it's because that's
the octave license and it has lots of convergence issues not to wonder because there are things
which are wrong and the author knew. And there are also nicer new things. So for example,
there is this implementation by S Philip which is basically he wrote a paper about how to implement
these things better and then he has the open source C++ implementation demonstrating how these
works. So this is much nicer. It is modern C++. It uses the eigenlibrary so that means you can
put in there any numerical types you want and it will work. It uses GMP or MAPMPFR for arbitrary
numerical position which is important to solve some of the three case filtering sign problems.
Again, unfortunately, piecewise constant response and weights. I think this would be easy to change.
It's just like an API thing. You could put in there a C++ lambda function. It should work.
This used to be GPL 3.0 later but it switched to BSD in April 2024.
Because this is the guy who knows how to write these things properly. This has good convergence.
It also has good runtime. It runs pretty fast and all the numerical methods,
all the tricks they are in a paper. You can go to the GitHub repository and the paper is linked in there.
And finally, I'm presenting here, there is this PM Remes. This is a vast implementation of
the Remes algorithm. It's based on ideas from S Philip. So basically based on ideas
from these good C++ implementation. It has a Python package. So you can just
pp install PM Remes. If you want to do some filter design, I don't know, on a Jupyter notebook,
on the command line or something. You can put in there are beta responses and weights. You can give it
the rast function. You can give it a Python function and anything which is a Python pullable.
You can pass it down. It's permissively licensed. So the main reason I started implementing
this thing was that the S Philip library was GPL and I needed something permissively licensed for
my SDR. It just happened that right when I finished this, he changed the license. I think it was
a pure coincidence because I was working on this and I only published the thing when I
had it completed and he switched, I think, a few days before I published the thing. So just
technical incidents. And it comes with very good documentation. That's also something else which is
the thing we are going to look at next. Well, actually not next, but in a couple slides.
So why are some of these implementations problematic when it comes to convergence issues?
Well, the thing is you have, so if you have many tabs, you have like many of these interpolation
points, you have tiny errors and all of these can be tricky for floating point calculations,
especially depending on how you do these calculations. And if some techniques you can use to improve
things, one of them is something called barricentric Lagrange interpolation, which is a way of doing
polynomial interpolation that basically reducing numerical error. The other one is TabiShep,
proxy route finding. That's actually a smart thing because when I have been floating things in the
previous slides, I have been saying, we now need to search the local extrema for the error,
blah, blah, blah, I'm floating red points on these things. But the question is, how do you find
the maximum of a function? And there are a whole bunch of ways to do it. In those plots, I was
basically cheating. I was pointing like a thousand points on the plot and taking the point where
it has the maximum value. It's not really the maximum, but it's just very close to it. It's the
maximum in my plot. And some of the implementations like site pie and guru radio, that's what they do.
They take a frying grid of points, they compute the error, that grid, and then off you go. The
point where you get the maximum value that's supposed to be the like the local extrema. And
surely when you have many, many points that's going to run into trouble. So these
CherishF proxy route finding is a clever thing. You approximate your function in a small interval
with a polynomial. And then you use like clever math tricks about finding routes of the derivative
of that polynomial, which is also a polynomial. And so the routes are, you know, when the function
is flat and so has the local extrema. And the thing is finding zeros of a polynomial can be done
as an eigenvalue problem. And so you know there are all these numerical methods to solve eigenvalue
problems. And of course you're going to use a higher precision numerical implementation. That
also helps. So that's the reason why there are better and worse implementations of these algorithm.
So now let's do some practical FIR design. All of the code showing is like in more details
in the Python package documentation. So you can go in there and you actually have like the code
for the plots and everything. So very simple decimation. Let's say I want to decimate by two or by
any other number. To do that without aliasing, I need to design a low pass filter whose
cutoff frequency is basically going to be the new Nyquist frequency after decimation. Okay,
fine. It also needs to have some transition bandwidth because it cannot go immediately from
Python to stop on. That is another design problem. And that is usually given as percentage of the
output bandwidth. So in this case, not point two means that 20% of the output bandwidth is going
to be transition bandwidth. So basically part of the spectrum is that we will not work properly
for your application. Your signal is supposed to be like in the center 80% if you do this.
And then the number of types, for example 35. So you basically use this code and so you get
something like this, which is the thing you wanted to have. Depending on the signals, your
supposed to have depending on applications, you might need better adenuation minus 60 dBs,
which is the thing we're getting here is reasonable in many cases. And also a question you might ask
is, well, I don't know what number of types I'm supposed to use to get something like this. So
you can try. There are also formulas like there's an approximate formula which will give you a
hint for the number of types. And then you can try with that and you can even eat rate. You know,
if the thing is doesn't satisfy my criteria, as I put one tap more and try again and so on and so
forth. And you have like the minimal number of types for you think. How can we improve these design?
So the thing is because of this equi-repled property for the error, with the previous design,
we are getting the same error on the passband and on the stopband. But oftentimes this is not the
thing we want to do. Because in the passband, it depends on the modulation and stuff,
but usually a 1% error is fine. In the stopband it depends on interferes signals. So for example,
you might need to have more attenuation in the stopband or you might want to have a lower number of
tops or something like that. So how can we achieve having different errors in the stopband and
the passband? That's where the weight comes in. And basically, I have the rule of them. There,
if you want 1% passband re-pull, then the weight depends on how much stopband attenuation you want.
For example, here, I think I have a weight of 10. And so I get my 1% error. So these
worse than in the previous plot, but I still get my 60 dB stopband rejection. So that's the thing.
You can play with the weights here. And then, and all the things you can do, Fred Harris proposes
everyone to do this. You can have these 1 over F response on the stopband. So you can see it here.
It's no longer flat, it's decreasing. And the intention of doing this is that when all of these
stopband aliases back into the passband because you are decimating, for here it does a mother, but
let's say you are decimating by 100. You have 100 copies of the noise floor which go on top of
your signal. So if you do not do this thing, then you have 100 times 60 decibels and that's
going to increase the noise floor you get on the output. So by doing this, then you have like
these sum of things which are decreasing and basically converges to constant value rather than
increasing as you increase the decimation. So that's like the reason why it's useful and we can
design this thing by saying at the beginning of the stopband the weight is going to be W at the end
of the stopband the weight is going to be W plus something larger. So the weight is supposed to be a linear
function. It has like these two values at the end points and we can define these with a triple.
The API is we'll understand that this is supposed to be a linear function having these two
values on the end points. What else? We can basically use all of these techniques to make a
PFB prototype filter and I'm not going to spend much time in the interest of the other examples explaining
what the PFB is. The only thing I'm going to say is that because a PFB can have many many
channels here's an example for 248 channels and then you also have like taps per channel. So this is
six taps per channel which is the lowest minimum number of taps per channel to get anything
slightly reasonable. We basically have 6,000 blah blah blah points and this is quite a lot. This is
like close to the point where the algorithm breaks. If you want to have like for example 4,000 and 96
points it will break but then what do you do? You design it for lower number of channels and then you
do linear interpolation of the taps. When you have these many taps they are like a nice continuous
function and you can just interpolate the taps to get a larger number of channels. So that's the trick
you need to know about. My favorite example is this one. CI C filters, maybe some of you are familiar
with these. Desimation filters, they are great because they only use summations. They do not need
multiplications. So if you are RFPJ or something which doesn't know how to multiply or it's very expensive
to multiply, these are ideal. The issue though is that the frequency response is crap. It's the great
thing here. So it droops. It only has good attenuation like in here and here in the notches.
Here the attenuation minus 50 is bad. So what do you do? You use a combination of a CFC filter which
is like doing the bulk of the decimation and then you put a decimation FIR. Usually decimate
by 4 or by 3 to make things nicer. And by makes things nicer what you can do is if you design the
CFIR filter properly, you can compensate for the droop that the CI C filter has in the
stop band. Sorry in the pass band. So you can you can see this here. Basically this is drooping
and it's easy to see in here. And this actually has some gain so that when we have the combination
it is actually flat. And let's see how we are doing in time. We are doing pretty well.
So something which is not that well known but it's like once you do the mathy subfuse
is the design that I need for a compensation filter depends on my DC ear which is the decimation
factor I apply on the CI C. Why is this relevant? The reason this is relevant is that when you
build this architecture where you have a CI C filter followed by a compensation filter it is very
easy to change the decimation factor on the CI C filter just by changing when you look at the
output. Basically you can be decimated by 10. And if you look at every other sample then suddenly
you are decimated by 20. For example this is how a USB or many other SDRs which perform
decimation on the FPGA work. That is the reason why your USB is able to do many different
decimation ratios just by changing the this DC the CI C decimation ratio. But formally the
like the transfer function we need for this compensation filter depends on DC there is bad but
luckily it's more or less the same thing as this this is like the small angle approximation for
the same function. So we get something independent of the DC. That is nice for example if you want to
build a USRP because people are going to change the CI C decimation and you do not want to
design any filter for each of them load them in the FPGA blah blah blah. So how can we compute
these filter with PM remaste. In this case we have a lambda function because the transfer function
we need is the inverse of this thing. So that's fine we can have lambda functions in Python and
we can put it in there in the desired response and that's how we get our nice passband shape which
compensates the FPGA filter group. So that's for a particular DC and if you want to have compensation
independent of DC it's the other formula. So just a different lambda function and this is basically
how it looks like because there is a narrower because we have made a small angle approximation in
that formula. Now this is no longer flat but in many cases the error is small enough that it doesn't matter.
What else? There is this thing called a Hilbert filter and different people call different
filters a Hilbert filter. If you take your radio this is not exactly the thing which is inside
the Hilbert filter but this is the main building block. So how ideas we have complex signals and
because they are complex they have both positive and negative frequencies. We want to do something
which looks rather stupid which is a frequency response of i. So that means plus 90 degrees on
positive frequencies and minus i. So that means minus 90 degrees phase shift on negative frequencies.
Why would we do anything like that? The reason is if you have this filter as a building block you
come for example cancel negative frequencies which is the thing that the GUNU radio Hilbert
filter block does basically using these and there's a plus one somewhere. So or plus i.
So this is another example where other implementations of parks, McLean and algorithm can be
awkward to use because they consider this thing as a special example. But in here it's just a
general example the only things you need to remember it's all taps and an all impulse response
symmetry that's how you get your imaginary frequency response and then you need to design an all
pass filter because ideally you want to have these like on the whole pass band of course that's
impossible because you cannot jump immediately from plus i to minus i as of course frequency
equals zero but we shall try. So then what changes not much n needs to be odd I don't know what
I used to hear you need to specify that the symmetry of the filter is odd because by default it's
going to be even and then you have this transition bandwidth which is basically the same as before.
And so we get the nice response we expect it which is basically flat in most of our frequency
ranks and then because you know you need to change to minus i we get a transition region.
Other more complex applications these are in my blog you can search for them to show that this is
quite general and useful. So there's something called a band h filter and this is available in
the radio it's the blog is called fl band h. I think the main reason to use this filter is to find
the frequency of a signal and it's usually used in a closed loop way so you basically have a
frequency look look that's how you tune in frequency to a signal. Interestingly you can also use it
for symbol clock recovery but you know that's like fret Harry's very clever stuff. So how does this
filter work? This is supposed to be for our RC waveforms that's root raised cosine waveform
like the thing I have in gray there and the idea so the ideal thing you wanted to have is this
green like dashed green transaction which is the derivative of the gray thing basically.
And you can notice that because this thing like has a sharp cut of here this derivative jumps
straight from plus one to zero so we already know we are not going to be able to get this with
any kind of filter we can make an approximation and the usual approximation which is what
fret Harry's normally uses and what's also in the radio is this blue thing which is a half
sine window function and this can be explicitly constructed in the time domain and it's nice.
The issue though is what happens if there is an adjacent signal in here. All that's problematic
because suddenly I'm ingesting part of the signal as part of my frequency error estimator.
So by using PM Remains we can tell Remains hey try to find the thing which is closed as close as possible
to this green thing which is the ideal and of course the thing will say yes yes but leave me some
gap where I try to go like very quickly from plus one to zero and so you leave it at
transition bandwidth but you know you can do something like this and depending on the location you
pay for example by having larger stop band worst stop band at innovation I would say nothing is free
right but in some occasions these PM Remains design might be better
and finally I was picking out polyphase filter banks so the thing with polyphase filter banks
they are like an FFT with a window but the window is an FIR filter rather than just multiplication by
a number and because you have all these flexibility of FIR filters you come make nice things
with nice properties. So for example let's say you want to do
spectrum estimation you want to measure power you might want to measure power at different
resolutions on things like that so you can do a design for that and the idea of this is if you look
here at how the frequency responses for the things are designed they basically overlap in such a way
that the thing is flat they added to something which is flat so for example here I am adding
being zero and being two which are these blue and green thing but if I add it's also flat so this is
well adjusted and what this gives me is now I can be power with a flat frequency response on a range
which is twice as coarse as my channel spacing by adding every other bin so that's the key idea
here you you have like these nice tools rather than you know just do an FFT and maybe use the
flat the window I will be nice it is kind of nice but not when you try for example what is the
summed power of a set of the frequency range you can make it way better with this and the design
trick here the thing you have is the freedom of word to put the cat of frequency if you set it
just right and you can make like an optimization problem where you search for this thing you can
achieve this effect where they are perfectly flat and that's it questions
we also hear once to difficult
and you see are new type as you radio
the different is your calculations
what is your-
Yes, so the question is about GR filter, where you can design filters, but not with the
remnants algorithm and face response, whether this is something you can tune with this algorithm.
So the answer would say is not you do not have any direct control over face response with
basically baked into the way that basically you are approximating a real transfer function.
And so the face response you get is going to be real and it's going to be what it is.
Okay, yeah, the question is whether my algorithm starts with equity space points.
And if there is a better choice, what is the reason?
So my implementation begins with equity space points on the x domain, which is not the same thing.
As the frequency domain, there are actually chairmanship points on the frequency domain.
There are smart choices in some cases and the paper bias Philippe talks about it.
But for simplicity, I did the thing which is which works nice in most of the cases.
And the only improvement you can gain by choosing more cleverly the points is runtime rather than convergence.
Thank you very much.