| Author |
 Topic  |
|
|
malgoska
Poland
36 Posts |
 Posted - 04/08/2013 : 12:11:39 PM
|
Origin Ver. OriginPro 9.0.0 (64-bit)SR1
Operating System: Windows 8 (64-bit)
Hi,
I'm trying to make a fitting function that would be a convolution of a piecewise function f(x)=C*(p1(x)+p2(x)), where
p1(x)=1/sqrt(a+x) for x>-a and x<2*a, elsewhere p1=0;
p2(x)=1/sqrt(a-x) for x>-2*a and x<a, elsewhere p2=0;
with gauss function.
Parameter a is one of the parameters to be fitted.
I followed this tutorial
http://www.originlab.com/www/helponline/Origin/en/mergedProjects/Tutorial/Tutorial/Fitting_with_Convolution_of_Two_Functions.html
but still have a problem with the range of p(x).
This is my function:
#pragma numlittype(push, TRUE)
#pragma warning(error : 15618)
#include <origin.h>
#include <ONLSF.H>
#include <fft_utils.h>
//----------------------------------------------------------
//
void _nlsfbrpake(
// Fit Parameter(s):
double a, double w, double C, double y0, double x0,
// Independent Variable(s):
double x,
// Dependent Variable(s):
double& y)
{
// Beginning of editable part
NLFitContext *pCtxt = Project.GetNLFitContext();
if ( pCtxt )
{
// Vector for the output in each iteration.
vector vX, vY;
static int nSize;
BOOL bIsNewParamValues = pCtxt->IsNewParamValues();
// If parameters were updated, we will recalculate the convolution result.
if ( bIsNewParamValues )
{
//Sampling Interval
double dx = 0.1;
vX.Data(-10, 10, dx);
nSize = vX.GetSize();
vector vF, vG, vp1, vp2, vBase;
vp1.SetSize(nSize);
vp2.SetSize(nSize);
//Function f(x)
for(int i=0; i<nSize; i++){
if(vX[i]>-2*(a-x0) && vX[i]<a-x0) vp1[i]=1/sqrt(a-(vX[i]-x0));
else vp1[i]=0;
}
for(int j=0; j<nSize; j++){
if(vX[j]<2*(a-x0) && vX[j]>-(a-x0)) vp2[j]=1/sqrt(a+(vX[j]-x0));
else vp1[j]=0;
}
vF=C*(vp1+vp2);
//Function g(x)
vG = 1/(w*sqrt(pi/2))*exp(-2*vX^2/w^2);
//Pad zeroes at the end of f and g before convolution
vector vA(2*nSize-1), vB(2*nSize-1);
vA.SetSubVector( vF );
vB.SetSubVector( vG );
//Perform circular convolution
int iRet = fft_fft_convolution(2*nSize-1, vA, vB);
//Truncate the beginning and the end
vY.SetSize(nSize);
vA.GetSubVector( vY, floor(nSize/2), nSize + floor(nSize/2)-1 );
//Baseline
vBase = y0;
//Fitted Y
vY = dx*vY + vBase;
}
//Interpolate y from x for the fitting data on the convolution result.
ocmath_interpolate( &x, &y, 1, vX, vY, nSize );
}
// End of editable part
}
It compiles but when I try to fit or simulate curve I get "unknown error" message and it does not iterate at all.
I've also tried to include the f(x) function as a sum of p1 and p2 without specifying any range (with no i/j loops) and suceeded to fit, but the function was calculated only on the common x range (-a;a).
What's wrong with function I've wrote? Any help?
Thanks,
Margaret
|
|
|
Sam Fang
293 Posts |
Posted - 04/10/2013 : 10:29:27 PM
|
The problem is due to the code:
----------------------
vector vX, vY;
----------------------
It should be as the tutorial:
---------------------------
static vector vX, vY;
---------------------------
Here we use the static because we store the convolution result for interpolation when parameters are not changed. In this way it can fast the speed.
Besides the above problem, I need say:
1. Maybe you need check your p1 and p2 definition. As I understand from your description, it should look like as follows?
for(int i=0; i<nSize; i++){
if(vX[i]>x0-a && vX[i]<x0+2*a) vp1[i]=1/sqrt(a+(vX[i]-x0));
else vp1[i]=0;
}
for(int j=0; j<nSize; j++){
if(vX[j]<x0+a && vX[j]>x0-2*a) vp2[j]=1/sqrt(a-(vX[j]-x0));
else vp2[j]=0;
}
2. You need change the code based on your x range.
double dx = 0.1;
vX.Data(-10, 10, dx);
If your x range is not located at [-10,10], you need adjust the code.
Sam
OriginLab Technical Services |
Edited by - Sam Fang on 04/10/2013 10:31:22 PM |
 |
|
| |
Topic |
|
|
|
fitting a convolution of two functions