Function got different results with the same inputs when directly called in MATLAB and called in exported .dll c++ library?

Hi everyone,

I'm working recently on some MATLAB-C++ combined programming work. I had successfully got my MATLAB function exported as a C++ dynamic library using the MATLAB Coder and have successfully called it in my C++ project. Everything seemed to be working fine.

However, then I found the results of calculation I got in the C++ program very strange so I checked it both in my C++ program and directly in MATLAB with the SAME function. To my surprise, they are different. For example, as I wanted to get a sum of some impulses with the following formular:

    for i=1:N               % elements of vector
        for j=ps(i)*fs+1:(ps(i)+pw(i))*fs
            curve(j)=curve(j)+lvl(i)*cos(double(2*pi*(j-ps(i)*fs-1)/fs*fo(i)))*exp(double(-alpha(i)*((j-ps(i)*fs-1)/fs)));
        end
    end

With the same input data N, ps, pw, fo, alpha and fs, what I got in my MATLAB function

 [sampled_time,sampled_curve,sampled_freq,sampled_spec] =  PIN_D_computation(rr,pw,ps,lvl,fo,fs,N,time_end,points)

was :

 10.0000   -6.7010    4.1943   -2.4138  1.2275    8.5102    2.7659   -6.7054 ...

while what I got corresponding in the exported C++ function

 extern void PIN_D_computation(emxArray_real_T *rr, emxArray_real_T *pw, emxArray_real_T *ps, emxArray_real_T *lvl, emxArray_real_T *fo, int32_T fs, int32_T N, int32_T time_end, int32_T points, emxArray_real_T *sampled_time, emxArray_real_T *sampled_curve, emxArray_real_T *sampled_freq, emxArray_real_T *sampled_spec);

was:

 10.0000   -9.8999    2.7735   -3.0868    0.5030    8.6276    0.5921   -4.5843...

I guess this is due to a wrong configuration in the MATLAB Coder, but I cannnot make it sure.

I am using a PC with Intel Core2 Duo CPU E8400 @ 3GHz CPU and Win7 64 bit OS, working with a 32 bit MATLAB 2012a and a 32 bit MSVC compiler for my programming. My configuration for the MATLAB coder is:

 Saturate on integer overflow: Yes
Support only purely-integer numbers: No
Support non-finite numbers: Yes
Code replacement library: Intel IPP(ISO)  (I've tried Intel IPP(GNU) and Intel IPP(ANSI) as well but none works)
Test hardware is the same as production hardware: Yes
Device Vendor: Intel
Device Type: x86/Pentium
Largest atomic integer size: long
Largest atomic floating-point size: double
Signed integer division rounds to: floor
Language: C++

Could anyone maybe help me please on this issue? Thanks a lot!

Answer 1

[José-Luis]

https://www.mathworks.com/matlabcentral/answers/57444-faq-why-is-0-3-0-2-0-1-not-equal-to-zero

The order in which operations are performed matters. Add to that issues of numerical stability, and it could explain the difference in you results. These issues tend to be amplified when there are exponents and trigonometry is involved. Try feeding slightly different values to your matlab or C functions to see if numerical stability might be an issue. If you get very different results you might be one step closer to identifying the culprit.

[José-Luis]

It could be, to answer your first question. But such large differences would have me suspect a numerical stability problem.

Is there a way to solve the problem: probably, but I can't think of an easy one. You could go step by step in the code and check where things go awry.

[Bin Han]

Sorry I removed my first comment before I realized that you had already replied it, just because I have found the problem. It was the default data-type of a multiply operation between a double variable and an integer variable. In MATLAB it is int*double=double, but in the exported C++ code it is int*double=int...

Answer 2

[Bin Han]

Problem solved.

It was the default data-type of a multiply operation between a double variable and an integer variable. In MATLAB it is int*double=double, but in the exported C++ code it is int*double=int...

As I changed the code from

 curve(j)=curve(j)+lvl(i)*cos(double(2*pi*(j-ps(i)*fs-1)/fs*fo(i)))*exp(double(-alpha(i)*((j-ps(i)*fs-1)/fs)));

to:

 curve(j)=curve(j)+lvl(i)*cos(2.0*pi*(double(j)-ps(i)*double(fs)-1.0)/double(fs)*fo(i))*exp(-alpha(i)*((double(j)-ps(i)*double(fs)-1.0)/double(fs)));

it began to work right! :)