Math Integral operation in FPGA
Simon
I have been several years experience on FPGA networking application
design, but I come up with the a question about math operation in
FPGA. I need to perform Integral operation. I don't know if there is
existing library which offer this function, or DSP? I think the last
option is to run the integral library in the embedded PowerPC in
Xilinx FPGA, but don't know how the performance is. Could anybody
indicate how to achieve this goal?
Many thanks,
Yixuan
rickman
An integral using continuous variables is essentially a multiply
operation with a continuous sum performed by an "integrator" circuit.
When using discrete time variable in digital logic, the multiply is
straight forward and often omitted since that is just a scale factor
and can be done elsewhere. The sum is now discrete and is a simple
matter of a continuous accumulation using an... accumulator.
If you can provide some details of your signal and exactly what you
want to do, maybe we can give some additional advice?
Rick
Sundar S
Integration is not just about summation of the discrete values. It could
also mean finding the area of the region encompassed by the given numbers
in Cartesian coordinates. Yixuan, if that's what you are looking for, you
may have to implement some adaptive quadrature algorithm.
The simplest way is to calculate the area as the curve builds for each and
every clock signal. Assuming that clock goes at the rate of one unit in
x-axis, you would need to evaluate int(n) = int(n-1) + 0.5 * diff(f(n),
f(n-1)) + min(f(n), f(n-1)), where f(n) is your function value at the nth
clock.
http://sunnyeves.blogspot.com/
Simon
Hi，Rick,
Regarding the integral I'm using it goes something like that:
If f(x) is the function to be integrated then
q = quadl(@(x) f(x),a,b,tol),
where quadl approximately integrates the f(x) from a to b using an
error tolerance tol.
Yixuan
Simon
Hi, Sundar,
I am thinking to use Xilinx Multiply Accumulator (MAC) core
http://www.xilinx.com/support/documentation/ip_documentation/mac.pdf
to do the integral. I don't know if it's feasible, or has anybody done
this, or have better option. Thanks,
Yixuan
rickman
That calculation looks complex, but isn't it really just
int(n) = int(n-1) + 0.5 * (f(n) - f(n-1))
The 0.5 times the difference assumes that the function is a straight
line between the two end points and when added to the min value is
just the average of the two points. This is an approximation, but
depending you your needs will be adequate.
I would argue that for an arbitrary function, there is no advantage to
using the average of each two points over just summing the points.
Consider points 0 to N where N is a large number.
N
< f(n) = f(1) + f(2) + ... + f(N)
<
1
N
< avg(f(n),f(n-1) = 0.5 f(0) + f(1) + f(2) + ... + 0.5 * f(N)
<
1
Notice that the only difference is that the average needs an extra
input point to calculate the first average and that the two end points
of the summation are halved. Numerically the difference between the
two calculations is 0.5 * (f(n) - f(0)). It appears to me to be a
very minuscule error to just add all the points without the complexity
of averaging. I would bet that for any value of N, 256 or over, this
error in the integral is much less than the error you get by the
original straight line average approximation.
In fact, whether you the average is correct or not depends on how you
picture the error formation. This is too complex to draw here, but if
you picture the sample as being centered in the region being
integrated by adding that value, then the error is only a function of
the second order components of f(x). To require an average
calculation you are assuming that the area being calculated for a
given point is the area *between* two points.
I could explain this more fully, but it is very hard to do without a
drawing.
Rick
Muzaffer Kal
wrote:
What you need is called numerical integration. Depending on the level
of accuracy needed, there are various ways of accomplishing it but at
the basic level it is an accumulator of the values of the function
needed to be integrated. Suppose you want to integrate a function f(t)
from t=a to t=b. The crudes integration would be to assume function f
has an average value of (f(b)+f(a))/s over this range so you integral
is (b-a) * (f(b)+f(a))/2. A better way is to divide the range [a b]
into small Dt slots (say N of them) and then iterate over this set
which makes the result (a+b)/N*(sum{(i=0;i<N) f(a+i*Dt)} which
basically says that you make small rectangles from the area under f(t)
and every rectange has an area of its width (Dt) times the value of
the function at the left hand x coordinate of the function. For
"smooth" functions this usually gives a pretty good result if your N
is "large enough" compared to "smoothness" of the function.
Here is a reference http://en.wikipedia.org/wiki/Numerical_integration
which has further references.
-
Muzaffer Kal
DSPIA INC.
ASIC/FPGA Design Services
Sundar S
@Rick,
0.5 * (f(n) - f(0)) is the straight line approximation of the entire
function. It gives the area of a triangle with base equal to 1 and height
equals to f(n) - f(0). But that's not the actual area as f(0) is somewhere
hanging above x-axis. So we have to consider the square that is formed
between the zeroth position and nth position so that you can get area below
the curve.
The calculation that I put up look complex, but if I assume that f(n) >
f(n-1), what I get is: for every clock tic: int(n) = int(n-1) + (f(n) -
f(n-1))>> 1 + f(n). Here as the curve grows by unit distance, the area is
calculated as the area of the rising triangle plus the area of the
rectangle that supports the triangle above x-axis.
@Yixuan,
You gave a vital information by telling that you use quadl of Matlab.
quadl is actually an implementation of lobatto quadrature. Implementing
that in FPGA may be really difficult. But what lobatto gives as output is
the area under the curve with certain approximation. In FPGA that
approximation anyway comes in terms of discretisation. Also in FPGA you
would not get a function to integrate; you usually get the values of the
function. In this case, the triangle-rectangle method would do.
rickman
I understand perfectly your original calculations. They are not
complex to understand. However, they are much more complex than
required. That is why I made my post explaining how in the infinite
series, your approach approximates the simple sum of the input
values. The only difference is that your approach can only produce
(N-1) output values for N input values and as a consequence, the final
output will not include the area of one column because you are always
subtracting out the first one.
Your approach errs in the thinking that the area is defined by the
points "surrounding" the area. What it is ignoring that the points
are *included* in the area. I suppose that initial int(0) could be
accounting for this somehow, without specifying how that
initialization is to be done.
I made a typo in my original equation which is equivalent to your more
complex calculation.
int(n) = int(n-1) + 0.5 * (f(n) - f(n-1))
should have been
int(n) = int(n-1) + 0.5 * (f(n) + f(n-1))
If you need me to, I can show you in a step wise manner how this is
equivalent to your equation,
int(n) = int(n-1) + 0.5 * diff(f(n), f(n-1)) + min(f(n), f(n-1))
Other than the end points of a series, your equation adds in half of
each data point at two separate times. So each point is summed to
produce the integral. Your calculation simply omits half of each end
point.
The mistake that is often made when dealing with discrete time samples
is thinking that they are the same as the instantaneous values of a
continuous function. In reality they are already integrals of the
amplitude and the sample period (1/f). That is why you only need to
sum them to obtain the integral over a series.
A simple sum is the correct way to calculate the integral of a
discrete time data series.
Rick
glen herrmannsfeldt
(snip, someone wrote)
<> http://sunnyeves.blogspot.com/
< That calculation looks complex, but isn't it really just
< int(n) = int(n-1) + 0.5 * (f(n) - f(n-1))
< The 0.5 times the difference assumes that the function is a straight
< line between the two end points and when added to the min value is
< just the average of the two points. This is an approximation, but
< depending you your needs will be adequate.
< I would argue that for an arbitrary function, there is no advantage to
< using the average of each two points over just summing the points.
< Consider points 0 to N where N is a large number.
< N
< < f(n) = f(1) + f(2) + ... + f(N)
< <
< 1
< N
< < avg(f(n),f(n-1) = 0.5 f(0) + f(1) + f(2) + ... + 0.5 * f(N)
< <
< 1
< Notice that the only difference is that the average needs an extra
< input point to calculate the first average and that the two end points
< of the summation are halved. Numerically the difference between the
< two calculations is 0.5 * (f(n) - f(0)). It appears to me to be a
< very minuscule error to just add all the points without the complexity
< of averaging. I would bet that for any value of N, 256 or over, this
< error in the integral is much less than the error you get by the
< original straight line average approximation.
Sure, but for small N it might be important.
< In fact, whether you the average is correct or not depends on how you
< picture the error formation. This is too complex to draw here, but if
< you picture the sample as being centered in the region being
< integrated by adding that value, then the error is only a function of
< the second order components of f(x). To require an average
< calculation you are assuming that the area being calculated for a
< given point is the area *between* two points.
There is a similar explanation of the uselessness of Simpson't
rule in Numerical Recipes. Simpson't rule has 1/3's and 2/3's that
look important. Using a similar argument, NR shows that it
is fancy decoration for the effect at the ends.
< I could explain this more fully, but it is very hard to do without a
< drawing.
-- glen
glen herrmannsfeldt
< I understand perfectly your original calculations. They are not
< complex to understand. However, they are much more complex than
< required. That is why I made my post explaining how in the infinite
< series, your approach approximates the simple sum of the input
< values. The only difference is that your approach can only produce
< (N-1) output values for N input values and as a consequence, the final
< output will not include the area of one column because you are always
< subtracting out the first one.
< Your approach errs in the thinking that the area is defined by the
< points "surrounding" the area. What it is ignoring that the points
< are *included* in the area. I suppose that initial int(0) could be
< accounting for this somehow, without specifying how that
< initialization is to be done.
This is confusing. The points have zero width and integrate
to zero. Left out is the actual limit for the integral. If it
includes one half a sample period before the first point, and
one half after the last point, then just add. If it doesn't,
then you need to divide the first and last point by two.
< I made a typo in my original equation which is equivalent to your more
< complex calculation.
< int(n) = int(n-1) + 0.5 * (f(n) - f(n-1))
< should have been
< int(n) = int(n-1) + 0.5 * (f(n) + f(n-1))
< If you need me to, I can show you in a step wise manner how this is
< equivalent to your equation,
< int(n) = int(n-1) + 0.5 * diff(f(n), f(n-1)) + min(f(n), f(n-1))
< Other than the end points of a series, your equation adds in half of
< each data point at two separate times. So each point is summed to
< produce the integral. Your calculation simply omits half of each end
< point.
< The mistake that is often made when dealing with discrete time samples
< is thinking that they are the same as the instantaneous values of a
< continuous function. In reality they are already integrals of the
< amplitude and the sample period (1/f). That is why you only need to
< sum them to obtain the integral over a series.
For the sampling theorem to work they must be samples of the
instantaneous value. Otherwise the result is filtered and
the results will be wrong when used as filter input (or almost
any other use.)
< A simple sum is the correct way to calculate the integral of a
< discrete time data series.
If you can manage the sample points in the center of the bins, yes.
Certainly that is preferred, but not always possible.
-- glen
rickman
Glen,
What you have said is a self-contradiction. If there is no averaging
at all (an impossibility - all measurements are made within a finite
time window) and the point measures a value at an exact point it time,
then there is no bin to be in the center of or at the edge of. They
are just measurements of a point in time. What would define the bin?
In fact, what *IS* a bin in this context?
It is a matter of perspective which I believe is shown clearly in the
calculations. As I showed above (or I tried to show) the two
calculations performed on an infinite time series are exactly
equivalent and when performed on a finite time series the only
difference is that one method includes all the samples with a
multiplier of 1 and the other uses a multiplier of 0.5 only on the end
points.
The distinction of looking at the "bins" as being around the sample or
between the samples is a red-herring and have no relation to the real
math of integrating a time sampled signal.
| | | |--+--|
| | |--+--| * |
| |--+--| * | * |
|--+--| * | * | * |
| * | * | * | * |
| * | * | * | * |
t0 t1 t2 t3
3 4 5 6
In this example, the + indicate the value of the signal f(t) at that
point in time. The integral calculation assumes that the value of f
(t) between the time samples is a straight line. By treating the
period as "1" and summing the values, you are in effect calculating
the areas shown above with the point centered in the rectangles. The
area on the left and right of the point between the top of the
rectangle and the straight line approximation have opposite signs
since one is making the area larger than it should be and other is
making is smaller. But they are equal in value and cancel.
The other method draws rectangles between the points with a height
equal to the average of each two points. This is a much more complex
calculation and gives the same result. There is no reason to
complicate the calculations using this calculation. Neither
calculation has a basis in the physics of the problem. They are just
equivalent numerical approximations.
Rick
rickman
Actually, the difference is that for a finite series, one approximates
the integral of the curve at N points and the other approximates the
integral of N-1 points. This is easy to see if you consider the width
of the area being approximated. For a constant value using two
points, the simple sum is twice the value of a single point and three
points gives a value half again as large as two points. Using the
average method you can't calculate a value for one point and the value
for three points is *twice* the value for two points. That's because
the area being approximated using three points is only two periods
wide, not three.
Time sampled series are not continuous functions and can not be
analyzed the same way.
Rick
glen herrmannsfeldt
(snip, I wrote)
<> This is confusing. The points have zero width and integrate
<> to zero. Left out is the actual limit for the integral. If it
<> includes one half a sample period before the first point, and
<> one half after the last point, then just add. If it doesn't,
<> then you need to divide the first and last point by two.
(after you wrote)
<> < The mistake that is often made when dealing with discrete time samples
<> < is thinking that they are the same as the instantaneous values of a
<> < continuous function. In reality they are already integrals of the
<> < amplitude and the sample period (1/f). That is why you only need to
<> < sum them to obtain the integral over a series.
<> For the sampling theorem to work they must be samples of the
<> instantaneous value. Otherwise the result is filtered and
<> the results will be wrong when used as filter input (or almost
<> any other use.)
(snip)
< What you have said is a self-contradiction. If there is no averaging
< at all (an impossibility - all measurements are made within a finite
< time window) and the point measures a value at an exact point it time,
< then there is no bin to be in the center of or at the edge of. They
< are just measurements of a point in time. What would define the bin?
< In fact, what *IS* a bin in this context?
I read both comp.dsp and comp.arch.fpga, and sometimes there is
overlap between them. An important part of digital signal processing
is the sampling theorem which, briefly, (and ignoring non-baseband
signals) states that a signal sampled at more than twice its highest
frequency component can be reconstructed from those samples.
For any problem involving sampled signals, the assumption, then,
is that it was properly filtered before sampling. As integration
is a averaging, or low pass filtering, process it tends to reduce
any error in the input data, including sampling error.
at twice its
< It is a matter of perspective which I believe is shown clearly in the
< calculations. As I showed above (or I tried to show) the two
< calculations performed on an infinite time series are exactly
< equivalent and when performed on a finite time series the only
< difference is that one method includes all the samples with a
< multiplier of 1 and the other uses a multiplier of 0.5 only on the end
< points.
Yes. If you add up (n+1) points and want the average value you
should divide by (n+1) not n. Consider the case where all sample
points have the same value. If the length of the continuous integral
is n, adding up (n+1) points will give a result (n+1)/n too large.
On the other hand, the (n+1) points could be over an interval of
length (n+2). If the length is n the sum of the weights of the
points should be n. (That is, for at least a 0th order approximation.)
If the length is n, one solution is to multiply the first and last
points by 0.5, another is to multiply all points by n/(n+1).
< The distinction of looking at the "bins" as being around the sample or
< between the samples is a red-herring and have no relation to the real
< math of integrating a time sampled signal.
< | | | |--+--|
< | | |--+--| * |
< | |--+--| * | * |
< |--+--| * | * | * |
< | * | * | * | * |
< | * | * | * | * |
< t0 t1 t2 t3
< 3 4 5 6
< In this example, the + indicate the value of the signal f(t) at that
< point in time. The integral calculation assumes that the value of f
< (t) between the time samples is a straight line. By treating the
< period as "1" and summing the values, you are in effect calculating
< the areas shown above with the point centered in the rectangles. The
< area on the left and right of the point between the top of the
< rectangle and the straight line approximation have opposite signs
< since one is making the area larger than it should be and other is
< making is smaller. But they are equal in value and cancel.
I agree. What is important is the total length of the region
being integrated. That could be n, n+1, n+2, or almost anything
else. (For uniform sampling, hopefully not too far from one of those.)
< The other method draws rectangles between the points with a height
< equal to the average of each two points. This is a much more complex
< calculation and gives the same result. There is no reason to
< complicate the calculations using this calculation. Neither
< calculation has a basis in the physics of the problem. They are just
< equivalent numerical approximations.
If the actual length of the integral isn't too important, I 100% agree.
If it is somewhat important, as it often is for not so large n, then
I don't agree. Numerical integration routines are often ranked on
the highest order polynomial they give an exact result for. (Even
though it is rare to actually numerically integrate polynomials.)
If you don't even get zero order, that isn't good. Still, in the
large n limit it isn't far off.
-- glen
glen herrmannsfeldt
(snip, someone wrote)
<> < int(n) = int(n-1) + 0.5 * (f(n) - f(n-1))
(snip)
< Actually, the difference is that for a finite series, one approximates
< the integral of the curve at N points and the other approximates the
< integral of N-1 points.
I would call it an integral of length N or N-1, but yes.
< This is easy to see if you consider the width
< of the area being approximated. For a constant value using two
< points, the simple sum is twice the value of a single point and three
< points gives a value half again as large as two points. Using the
< average method you can't calculate a value for one point and the value
< for three points is *twice* the value for two points. That's because
< the area being approximated using three points is only two periods
< wide, not three.
< Time sampled series are not continuous functions and can not be
< analyzed the same way.
Sampled series are often approximations of continuous functions.
For band limited continuous functions the sampled series can exactly
represent the continuous function. (That is, sampled by not
quantized.) As I said previously (before reading your post), it
is usual to indicate the highest order polynomial that a numerical
integration routine gives an exact answer for. (Ignoring any
rounding or quantization errors.) Within certain constraints,
time sampled series can be analyzed in a way comparable to
coninuous functions.
-- glen