No BSD License  

3.75

3.8 | 4 ratings Rate this file 23 Downloads (last 30 days) File Size: 3.09 KB File ID: #8644

Fast Trapezoidal Integration

by

 

05 Oct 2005 (Updated )

Similar to TRAPZ but (up to) 16x-times faster. Equally spaced data NOT required.

| Watch this File

File Information
Description

F=trapzf(X,Y) computes the integral of Y with respect to X using the trapezoidal method. X and Y must be vectors of the same length.
For trapzf to work properly:
X must be a STRICTLY monotonically increasing vector.

Example:
>> x=[1:1:1000];
>> y=log(sqrt(x+1.001)-1.001);
>> area=trapzf(x,y);

Warning: not much in the way of error checking, since this slows things down, so pay attention to the argument passed to the function!!!

Installation (Windows users):
------------------------------

Simply copy the trapzf.dll file in a directory recognizable by MATLAB. You won't need trapzf.c to perform computations, but it is useful if you want to read the user instructions or to customize your code.
If you modify trapzf.c then you must run the following command

mex trapzf.c

in order to obtain the new dll file.

Installation (other OS users):
-------------------------------

A precompiled c-mex file is not provided.
You have to compile trapzf.c with your favourite C compiler (write "mex -setup" at the MATLAB prompt), then write

mex trapzf.c

in order to obtain the corresponding mex-file (see the table below)

        OS ---> Extension
      ------------------------------
        sol2, SunOS 5.x ---> .mexsol
        hpux ---> .mexhpux
        hp700 ---> .mexhp7
        ibm_rs ---> .mexrs6
        sgi ---> .mexsg
        alpha ---> .mexaxp
        glnx86 ---> .mexglx

PLEASE RATE MY SUBMISSION!

MATLAB release MATLAB 6.5 (R13)
Other requirements Non-Windows users need to compile the c-source file; use the "mex -setup" command utility to search for the available compiler(s).
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (5)
19 Sep 2008 Lance Optican

trapzf does not allow for unevenly sampled data, even though the comment claims it does. Matlab's trapz does allow for uneven data. To correct this, the following change must be made to the C code of trapzf, here called trapzfu. Note that the speedup on a 64-bit WinXP computer of trapfzu over trapz was about 10-fold. Unrolling the loop in trapzfu2.c provides another 20% improvement in speed.

--- trapzfu.c ------
/* F=trapzfu(X,Y) computes the integral of Y with respect to X using
the trapezoidal method. X and Y must be vectors of the same length.
For trapzf to work properly:
X must be a STRICTLY monotonically increasing vector.
X does not need to be evenly spaced.

C. Quaia, National Eye Institute, NIH, DHHS
*/

#include "mex.h"

void trapzfu( double *area, double *xv, double *yv, int N) {
register int i;

*area = 0;
for (i = 0; i < (N-1); i++) {
*area = *area + (xv[i+1] - xv[i]) * (yv[i] + yv[i+1]);
}
*area = *area / 2;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
int N, N1, N2;
double *area;
double *xv, *yv ;

N1 = mxGetNumberOfElements(prhs[0]);
N2 = mxGetNumberOfElements(prhs[1]);
if (N1 != N2)
mexErrMsgTxt("Input x and y must have the same length");
else
N = N1;
if (N < 2)
mexErrMsgTxt("The minimum length is 2");

xv = mxGetPr(prhs[0]);
yv = mxGetPr(prhs[1]);
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
area = mxGetPr(plhs[0]);

trapzfu(area, xv, yv, N);
}

------------ trapzfu2.c --------
/* F=trapzfu2(X,Y) computes the integral of Y with respect to X using
the trapezoidal method. X and Y must be vectors of the same length.
For trapzf to work properly:
X must be a STRICTLY monotonically increasing vector.
X does not need to be evenly spaced.

L. Optican, NEI/NIH/DHHS

*/

#include "mex.h"

void trapzfu2( double *area, double *xv, double *yv, int N) {
register int i, j,r;

*area = 0;
i = 0;
r = (N - 1) % 4;
for(; i < r; i++) {
// *area += (xv[i+1] - xv[i]) * (yv[i] + yv[i+1]);
*area += (xv[i+1] - xv[i]) * (yv[i+1] + yv[i]);
}

for(; i <= (N-4); i += 4) {
*area += (xv[i+1] - xv[i]) * (yv[i+1] + yv[i]);
*area += (xv[i+2] - xv[i+1]) * (yv[i+2] + yv[i+1]);
*area += (xv[i+3] - xv[i+2]) * (yv[i+3] + yv[i+2]);
*area += (xv[i+4] - xv[i+3]) * (yv[i+4] + yv[i+3]);
}
*area = *area / 2;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
int N, N1, N2;
double *area;
double *xv, *yv ;

N1 = mxGetNumberOfElements(prhs[0]);
N2 = mxGetNumberOfElements(prhs[1]);
if (N1 != N2)
mexErrMsgTxt("Input x and y must have the same length");
else
N = N1;
if (N < 2)
mexErrMsgTxt("The minimum length is 2");

xv = mxGetPr(prhs[0]);
yv = mxGetPr(prhs[1]);
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
area = mxGetPr(plhs[0]);

trapzfu2(area, xv, yv, N);
}

06 Mar 2008 lll hhh

not good

14 Feb 2008 mm prabhu2003

thanks

12 Mar 2006 sind baad  
08 Feb 2006 Jay Leno  
Updates
19 Oct 2005

no modifications: I just added my request to give a rate to this file.

14 Nov 2005

Now the inputs can be either row or column arrays.

15 Nov 2005

no significant changes

27 Dec 2005

this version has been created with the LCC compiler (instead of Visual Studio 6); further, the installation steps description have been improved.

28 Dec 2005

Requirements added

22 Feb 2006

no significative changes

20 Apr 2006

No changes

09 Aug 2006

no significative changes

11 Oct 2006

1) Bug corrected: now trapzf(x,y) returns zero if length(x)=length(y)=1.
2) A bit of effort on dimension check of x and y has been added (just a little...to not decrease the algorithm speed)

Contact us