Bug of matlab coder (2013a) for IFFT?

Asked by laoya about 11 hours ago
Latest activity Commented on by laoya about 7 hours ago

Dear all,

I am trying to use the matlab coder to generate C++ code for IFFT. However, results from matlab and C++ code have greate difference. Could anyone help me to take a look at this problem?

The matlab code to do IFFT is as follows:

z=ifft(y,[],1);

The data of y are as follows:

487924.35768878186	
481026.24091831001	
478940.94803639117	
480926.64214804117	
478909.92810891091	
475591.66558611894	
472961.12419697485	
478082.98773687298	
475963.06721005152	
471624.91585664771	
473035.01798060897	
470885.04387918825	
465204.72500935785	
465559.98854017764	
468284.20173525374	
467779.54004903068	
467358.62894582795	
461658.53205369966	
461058.87787927384	
461185.53453021083	
455868.37107306172	
462911.95616388740	
459834.12223788921	
458194.44688116410	
456595.68046306813	
451450.72108931403	
448915.17218837264	
451224.11133871856	
448603.78107835044	
442943.27244258055	
445878.56931863062	
451008.71128154272	
451212.54731590918	
455779.62766682287	
454335.89646260266	
450941.74302688910	
451440.15563114837	
451448.07972477260	
449064.57034377591	
448331.01187180565	
449871.73269351595	
449496.69114039827	
445229.11575907044	
444080.89526586136	
443461.78518676467	
437449.62073759746	
443516.19085396524	
438574.10116169654	
434403.16106846917	
433751.71037962160	
426817.29095258669	
422633.17624841497	
418137.67043422104	
413457.43016094179	
411045.60664053162	
403724.29173170537	
392643.02660994284	
380096.57725000178	
357723.60756075662	
332111.19753204676	
317106.73298398359	
311326.13889686763	
311656.30946454400	
310753.51039134728	
313637.04296473821	
311065.22361899650	
310994.64764692046	
309178.16192473430	
310382.78521443135	
310153.53409929451	
308257.54974676436	
309423.90546237369	
310570.45094392239	
308946.62377444573	
310466.56800921442	
309782.13247536187	
320356.73878161400	
347570.20226357155	
368694.03200706805	
384869.13643517962	
398566.44765287120	
407737.20091712585	
412582.04329780280	
418070.05794431345	
422272.98431790201	
424134.63093139383	
430451.22874984005	
433175.44194491615	
437468.43240701430	
438150.54869395244	
439138.45124420297	
443049.79674182995	
441973.82723235723	
440130.96007651964	
447706.58685186348	
449457.42500166676	
454555.58068421762	
449613.45878018625	
452690.29414154071	
451600.11924471724	
452466.32575548667	
454145.23503918055	
449159.24071435159	
456451.88715437474	
450130.32872447267	
450516.57997100998	
452118.30987127440	
449933.06388969702	
452284.90910795948	
452236.75252272276	
453562.39540488424	
455035.18161905993	
456251.33671980369	
457868.23836029979	
460014.57253253541	
462557.04696245730	
461771.33764667553	
462317.90683617146	
466403.38912295533	
467626.46975267958	
466153.36142087687	
467724.10360542370	
470285.06758713548	
469912.31306916947	
469812.71429890068	
469016.11740732624	
470695.70313803683	
473460.18103648809	
475000.22541118169	
476990.88013428752	
477611.31089565484	
475078.08124162792	
479902.43694122758	
478081.34493697528	
480799.95328534150	
484503.82281941688	

Thanks, Tang Laoya

0 Comments

laoya

Products

No products are associated with this question.

2 Answers

Answer by Sean de Wolski about 10 hours ago
Edited by Sean de Wolski about 10 hours ago

I'm not surprised to see a small difference. Building a mex file from this:

function y = doIFFT(x)
      y = ifft(x,2^nextpow2(numel(x)),1);
      %I had to add the nextpow2 part to get it to build
  end

Then looking at the differences:

norm(doIFFT(x)-doIFFT_mex(x))
ans =
   4.695591604029106e-11

The difference is very small. In fact it's close to:

eps(487924.35768878186)

So this is just a small numerical round-off difference in the approach taken by the libraries MATLAB sits on and the generated c-code.

0 Comments

Sean de Wolski
Answer by laoya about 10 hours ago
Edited by laoya about 10 hours ago

Hi Sean de Wolski,

Thanks for your quick reply. However, from here the difference is not so small. The results from matlab are: 1.0e+005 *

4.313512092938006                     
0.300430732843556 - 0.005783095388433i
-0.147813810099868 + 0.000951173782191i
0.133108569308995 - 0.003373860065190i
-0.061656447655613 + 0.003328454544711i
0.027192760859384 - 0.000603671006620i
-0.005797878506124 - 0.000492192018536i
-0.001061055121786 + 0.001998424724582i
0.021577859002924 - 0.002456894682518i
-0.021811429158957 + 0.003207842695919i
0.017735324252322 - 0.003028416146602i
-0.014305192477317 - 0.000474115939438i
0.013515213478997 - 0.002166798336642i
-0.004684409887853 + 0.003175414468610i
-0.000443185935397 + 0.000308035524811i
0.003853153787452 - 0.002023631130418i
-0.005850698801703 + 0.000216651950420i
0.012856303304599 - 0.001080645604734i
-0.005791634751852 + 0.002526085886894i
0.005578657992319 - 0.000926823744157i
-0.000112906707825 + 0.001805996486923i
0.001274493397202 + 0.000127337940853i
-0.001556944764524 - 0.000427023062860i
-0.002102718086437 - 0.000349736278996i
0.002944493706193 - 0.003170105190536i
-0.001861635715701 - 0.000525536225213i
0.002975662968139 - 0.000719539330124i
-0.002191858687055 + 0.000974062580574i
0.000491769705110 - 0.001428164102981i
-0.000607238680120 - 0.000256915753116i
0.002503276464045 + 0.001201049865073i
0.000502580943466 - 0.000828435382671i
-0.001255667660083 + 0.000254113854892i
0.001911608261996 - 0.002132432317272i
-0.000122452068500 - 0.001590327883673i
0.000678876401114 - 0.002190333520453i
0.000850110098716 + 0.001294448469562i
0.002221809217307 - 0.000687101621501i
0.001690063348324 - 0.000139630162640i
0.001283562042216 + 0.000855341907948i
0.002625543849526 - 0.001037104423102i
0.000443377947227 + 0.002073096558817i
-0.000133656504791 - 0.000054353863241i
-0.001507985869583 - 0.001123830171823i
0.001803796928052 - 0.000319702526871i
0.000802884469623 + 0.001769241259875i
0.000017848403353 + 0.000495263516297i
-0.000878046846166 - 0.002090322061183i
-0.001684884158707 + 0.000640259092086i
0.000349058948051 - 0.002674201272202i
0.001317964739399 - 0.000682460501570i
0.000136079066329 - 0.000196295344922i
-0.000751139021633 - 0.000771625053355i
0.001165570869108 + 0.001735364213629i
-0.000174024655019 + 0.001713751491137i
-0.000142209151696 - 0.001924788489301i
-0.001616076044371 - 0.000170224894442i
0.000779708646816 + 0.001697818033402i
0.002037664934465 - 0.000305846247259i
0.000752313211467 - 0.000118969236097i
-0.001010826774948 + 0.001686222917862i
-0.000086432763696 + 0.000856642401964i
0.001333583110647 + 0.000499768614078i
-0.000941810898597 - 0.000492021314954i
0.000162814303691 + 0.000686033402377i
0.001251368703072 + 0.000230732674254i
0.000737206223204 + 0.000479604300745i
-0.000537515690459 + 0.000473585154103i
0.000927698765764                     
-0.000537515690459 - 0.000473585154103i
0.000737206223204 - 0.000479604300745i
0.001251368703072 - 0.000230732674254i
0.000162814303691 - 0.000686033402377i
-0.000941810898597 + 0.000492021314954i
0.001333583110647 - 0.000499768614078i
-0.000086432763696 - 0.000856642401964i
-0.001010826774948 - 0.001686222917862i
0.000752313211467 + 0.000118969236097i
0.002037664934465 + 0.000305846247259i
0.000779708646816 - 0.001697818033402i
-0.001616076044371 + 0.000170224894442i
-0.000142209151696 + 0.001924788489301i
-0.000174024655019 - 0.001713751491137i
0.001165570869108 - 0.001735364213629i
-0.000751139021633 + 0.000771625053355i
0.000136079066329 + 0.000196295344922i
0.001317964739399 + 0.000682460501571i
0.000349058948051 + 0.002674201272202i
-0.001684884158707 - 0.000640259092086i
-0.000878046846166 + 0.002090322061183i
0.000017848403353 - 0.000495263516297i
0.000802884469623 - 0.001769241259875i
0.001803796928052 + 0.000319702526871i
-0.001507985869583 + 0.001123830171823i
-0.000133656504791 + 0.000054353863241i
0.000443377947227 - 0.002073096558817i
0.002625543849526 + 0.001037104423102i
0.001283562042216 - 0.000855341907948i
0.001690063348324 + 0.000139630162640i
0.002221809217307 + 0.000687101621501i
0.000850110098716 - 0.001294448469563i
0.000678876401114 + 0.002190333520453i
-0.000122452068500 + 0.001590327883673i
0.001911608261996 + 0.002132432317272i
-0.001255667660083 - 0.000254113854892i
0.000502580943466 + 0.000828435382671i
0.002503276464045 - 0.001201049865073i
-0.000607238680120 + 0.000256915753116i
0.000491769705110 + 0.001428164102981i
-0.002191858687055 - 0.000974062580574i
0.002975662968139 + 0.000719539330124i
-0.001861635715701 + 0.000525536225213i
0.002944493706193 + 0.003170105190536i
-0.002102718086437 + 0.000349736278996i
-0.001556944764524 + 0.000427023062860i
0.001274493397202 - 0.000127337940853i
-0.000112906707825 - 0.001805996486923i
0.005578657992319 + 0.000926823744157i
-0.005791634751852 - 0.002526085886894i
0.012856303304599 + 0.001080645604734i
-0.005850698801703 - 0.000216651950420i
0.003853153787452 + 0.002023631130418i
-0.000443185935397 - 0.000308035524811i
-0.004684409887853 - 0.003175414468610i
0.013515213478997 + 0.002166798336642i
-0.014305192477317 + 0.000474115939438i
0.017735324252322 + 0.003028416146602i
-0.021811429158957 - 0.003207842695919i
0.021577859002924 + 0.002456894682518i
-0.001061055121786 - 0.001998424724582i
-0.005797878506124 + 0.000492192018536i
0.027192760859384 + 0.000603671006620i
-0.061656447655613 - 0.003328454544711i
0.133108569308995 + 0.003373860065190i
-0.147813810099868 - 0.000951173782191i
0.300430732843556 + 0.005783095388433i

But results from C++ code are as follows:

   72422.083779590306     -21175.927029026894     	
   11125.992351540337     3106.1219836735977     	
   6843.8805862924828     3554.4005362181247     	
   4356.1650364638035     2570.7688328224426     	
   3420.3547292047629     0.31787260351634089     	
   4700.1741689654591     -2577.2663742471132     	
   7910.0731580255633     -3855.9217257203227     	
   13419.434500006109     -4115.2819077652211     	
   -13282.696522193037     -30704.471141475420     	
   7971.4850412689193     3265.5656871257138     	
    5862.2579664486675     4181.4237102004063     	
    3696.9860231590023     2941.4488545757381     	
    3051.9471877097153     129.42359763279092     	
    4852.0985869896149     -2863.1910672093150     	
    9114.7261030580539     -5404.2257888381619     	
    48675.808343791687     -29709.095883131220     	
    -34398.075003605700     -18366.992612829283     	
    5138.7854450924133     4996.7355815242800     	
    4654.2058745860568     4864.4564801580163     	
    1543.8705522082814     4071.7857881112859     	
    1819.7875797722779     457.90605746160725     	
    2830.8298912636742     -2275.1392367139283     	
    9716.4329594613992     -7783.0865968464159     	
    -3630.7897718976387     9142.7833067804131     	
    -8845.6457668685016     414.34287820314398     	
    4992.1599224602105     4868.9723525938753     	
    3229.5822994476157     4878.3164657138086     	
    736.67044680706999     3430.1969587796079     	
    957.35710922778412     -1031.5653556290436     	
    2480.5289862458612     -5208.5851937693360     	
    25700.781958828022     -39844.385866482604     	
    -394.49267492773350     7591.7737507502407     	
    -20665.267900238770     27776.650484342088     	
    -1473.5210380044171     4785.8222036708621     	
    60.051900129337277     7119.4630475685044     	
    -830.62218704141992     3490.4524725087017     	
    -1074.6118697358522     -875.00645669774576     	
    2877.6149243907121     -9498.2795773337984     	
    -9229.4478008981059     18304.491631238121     	
    308.96239530547393     4249.9275609692568     	
    -2404.7107375548189     8985.3059745066676     	
    2132.7840149464282     4447.5810777565584     	
    -477.38110141003563     4978.7152523136738     	
    -2549.9004910493386     2030.8234215058337     	
    -3073.6612221284904     -4054.7995961807842     	
    5432.9358533986097     -38643.332288489612     	
    -3605.4539005091733     9775.7604208193752     	
    45.552641515049103     4113.9140568298253     	
    221.05280647733633     12207.662006532393     	
    -278.76886434015228     7710.8290184780490     	
    -1890.2508641515733     3916.7204912494567     	
    -4536.8419186779274     2104.3932757967714     	
    -5059.1810510076248     -8576.6255176689756     	
    -6110.5398786129972     29323.824846673553     	
    -2794.5666980626256     4779.6064599904876     	
    -198.13055532286543     1728.5801399849456     	
    1797.5246407960512     3609.7437062598865     	
    -576.50490459068919     3309.6370867571586     	
    -3325.1497702358815     1971.7321719294284     	
    -6523.0645758753290     -1698.9721656553158     	
    -12240.582997565247     -30058.889007986600     	
    -2717.1305815325768     11518.450610493081     	
    -1454.2548795338439     4192.0167255739652     	
    52.786090168379495     1026.7177630023957     	
    -10338.642087783273     13642.323636409885     	
    -6812.9003180914069     2866.9799795768654     	
    -3792.8505170487479     -230.78805726075026     	
    -7213.3223499798705     -7679.6190964047646     	
    321.28615576961795     36259.999009415274     	
    -4383.7395961326183     3415.7340629150949     	
    -2309.0277783778693     -1088.5516338423911     	
    -180.79148605519302     -1836.5173221714188     	
    743.63999390828008     -6101.2152527950202     	
    -1134.2989609607107     -524.23819381373767     	
    -4792.5955569409125     -2348.0148101847226     	
    -14766.177531123119     -21182.107722402230     	
    -439.08920376328257     10344.607622734815     	
    -2542.4461870583073     599.14559236859304     	
    -845.19793586519643     -2599.3865384547471     	
    1187.5292526634566     -2785.6446093624072     	
    -9038.0486912758315     -7065.9521600213020     	
    80.625318829726467     10.714344970158622     	
    -5653.8415272376815     -4732.4052208226767     	
    21190.907021163606     37007.587554493541     	
    -852.94447933068420     3637.7765490408310     	
    -2098.0122985887024     -2686.0869132157095     	
    261.25956450206797     -4322.5060757986057     	
    1970.2539612731198     -4518.7374882269842     	
    -7274.3941270798468     -4518.5287317769698     	
    -330.11620341411918     -2225.0447706025757     	
    -11852.288982194095     -10043.221245743609     	
    6679.5226265230413     7775.3855646021684     	
    1377.9236897985454     1025.0031917079912     	
    113.67914951110266     -3664.2480443867912     	
    3090.0901443598441     -4740.7616094401865     	
    4285.8937306032358     -4849.8787958145267     	
    -24356.495771448343     14480.960943266133     	
    -5295.8195730284415     -3506.5828383575276     	
    40818.185859963669     30465.471866954791     	
    3565.6822875299658     3109.8019153641180     	
	    997.99930535255316     -1037.8513161753970     	
	    1596.4431752630942     -4058.2212596692852     	
	    2371.8530430438618     -6235.0706981591920     	
	    5157.0953715735923     -5258.0134317301263     	
	    -5483.7061459089282     6393.6879372332833     	
	    -2472.6887179699484     -4024.0136908333516     	
	    11191.498476511286     4243.0178602305305     	
	    4326.6288886882658     1068.9429182750932     	
	    2523.3147519425652     -1914.9842759916771     	
	    4750.3542365514413     -4293.9197562996178     	
	    4884.6046593881738     -6045.4671383659079     	
	    7723.2247447638074     -4995.1376654136848     	
	    -8572.1840474269666     38548.652251609477     	
	    57941.834946725779     12129.156253256031     	
	    8352.6956653903162     2187.0454560172848     	
	    3106.6583415751998     386.09826727500115     	
	    3429.7058523298415     -1988.1225111991816     	
	    1628.0516246254488     -4714.1787317279168     	
	    5902.3957729349340     -5738.5628963074414     	
	    8049.8073106827114     -4652.4255896710165     	
	    9145.3967079252434     32125.876528317396     	
	    12010.057611526385     766.20021100067629     	
	    8174.8653008944102     1558.3663819684973     	
	    4622.4202209866007     570.29637943909381     	
	    5678.5596863576793     -2317.1898630674232     	
	    3994.7403360780213     -4371.9898884749864     	
	    8688.5972400898154     -5958.6490558767273     	
	    10928.004975388818     -5130.2471041175613     	
	    2965.8936925455027     8488.5318855891073     	
	    1091.4535707318753     -1704.1399129821730     	
	    -3792.8505170487479     -230.78805726075026     	
	    -7213.3223499798705     -7679.6190964047646     	
	    0.00000000000000000     0.00000000000000000     	
	    0.00000000000000000     0.00000000000000000     	
	    0.00000000000000000     0.00000000000000000     	
	    0.00000000000000000     0.00000000000000000     	

C++ codes generated by matlab coder are:

/// ifft.c

/* * ifft.c * * Code generation for function 'ifft' * * C source code generated on: Thu Aug 01 16:42:35 2013 * */

/* Include files */

#include "rt_nonfinite.h"

#include "Dissociation.h"

#include "ifft.h"

#include "Dissociation_emxutil.h"

#include "Dissociation_rtwutil.h"

/* Function Definitions */

/* * */ void b_ifft(const emxArray_creal_T *x, emxArray_creal_T *y)

{

uint32_T sz[2];
int32_T ju;
int32_T minval;
emxArray_real_T *costab1q;
int32_T nd2;
int32_T ixDelta;
int32_T nRowsD2;
int32_T nRowsD4;
int32_T lastChan;
real_T e;
int32_T k;
emxArray_real_T *costab;
emxArray_real_T *sintab;
int32_T ix;
int32_T chanStart;
int32_T i;
boolean_T tst;
real_T temp_re;
real_T temp_im;
int32_T iDelta;
int32_T iDelta2;
int32_T iheight;
int32_T ihi;
real_T twid_im;
for (ju = 0; ju < 2; ju++) {
  sz[ju] = (uint32_T)x->size[ju];
}
ju = y->size[0] * y->size[1];
y->size[0] = x->size[0];
y->size[1] = (int32_T)sz[1];
emxEnsureCapacity((emxArray__common *)y, ju, (int32_T)sizeof(creal_T));
if ((x->size[0] == 0) || (x->size[1] == 0)) {
} else {
  minval = x->size[0];
  b_emxInit_real_T(&costab1q, 2);
  nd2 = (x->size[0] - minval) + 1;
  if (1 >= nd2) {
    ixDelta = 1;
  } else {
    ixDelta = nd2;
  }
    ju = x->size[0];
    nRowsD2 = ju / 2;
    nRowsD4 = nRowsD2 / 2;
    lastChan = x->size[0] * (div_s32(x->size[0] * x->size[1], x->size[0]) - 1);
    e = 6.2831853071795862 / (real_T)x->size[0];
    ju = costab1q->size[0] * costab1q->size[1];
    costab1q->size[0] = 1;
    costab1q->size[1] = nRowsD4 + 1;
    emxEnsureCapacity((emxArray__common *)costab1q, ju, (int32_T)sizeof(real_T));
    costab1q->data[0] = 1.0;
    nd2 = nRowsD4 / 2;
    for (k = 1; k <= nd2; k++) {
      costab1q->data[k] = cos(e * (real_T)k);
    }
    for (k = nd2 + 1; k < nRowsD4; k++) {
      costab1q->data[k] = sin(e * (real_T)(nRowsD4 - k));
    }
    b_emxInit_real_T(&costab, 2);
    b_emxInit_real_T(&sintab, 2);
    costab1q->data[nRowsD4] = 0.0;
    nd2 = (costab1q->size[1] - 1) << 1;
    ju = costab->size[0] * costab->size[1];
    costab->size[0] = 1;
    costab->size[1] = nd2 + 1;
    emxEnsureCapacity((emxArray__common *)costab, ju, (int32_T)sizeof(real_T));
    ju = sintab->size[0] * sintab->size[1];
    sintab->size[0] = 1;
    sintab->size[1] = nd2 + 1;
    emxEnsureCapacity((emxArray__common *)sintab, ju, (int32_T)sizeof(real_T));
    costab->data[0] = 1.0;
    sintab->data[0] = 0.0;
    for (k = 1; k < costab1q->size[1]; k++) {
      costab->data[k] = costab1q->data[k];
      sintab->data[k] = costab1q->data[(costab1q->size[1] - k) - 1];
    }
    for (k = costab1q->size[1]; k <= nd2; k++) {
      costab->data[k] = -costab1q->data[nd2 - k];
      sintab->data[k] = costab1q->data[(k - costab1q->size[1]) + 1];
    }
    emxFree_real_T(&costab1q);
    ix = 0;
    for (chanStart = 0; chanStart <= lastChan; chanStart += x->size[0]) {
      ju = 0;
      nd2 = chanStart;
      for (i = 1; i < minval; i++) {
        y->data[nd2] = x->data[ix];
        nd2 = x->size[0];
        tst = TRUE;
        while (tst) {
          nd2 >>= 1;
          ju ^= nd2;
          tst = ((ju & nd2) == 0);
        }
        nd2 = chanStart + ju;
        ix++;
      }
      y->data[nd2] = x->data[ix];
      ix += ixDelta;
      ju = (chanStart + x->size[0]) - 2;
      if (x->size[0] > 1) {
        for (i = chanStart; i <= ju; i += 2) {
          temp_re = y->data[i + 1].re;
          temp_im = y->data[i + 1].im;
          y->data[i + 1].re = y->data[i].re - y->data[i + 1].re;
          y->data[i + 1].im = y->data[i].im - y->data[i + 1].im;
          y->data[i].re += temp_re;
          y->data[i].im += temp_im;
        }
      }
      iDelta = 2;
      iDelta2 = 4;
      k = nRowsD4;
      iheight = 1 + ((nRowsD4 - 1) << 2);
      while (k > 0) {
        i = chanStart;
        ihi = chanStart + iheight;
        while (i < ihi) {
          nd2 = i + iDelta;
          temp_re = y->data[nd2].re;
          temp_im = y->data[nd2].im;
          y->data[i + iDelta].re = y->data[i].re - y->data[nd2].re;
          y->data[i + iDelta].im = y->data[i].im - y->data[nd2].im;
          y->data[i].re += temp_re;
          y->data[i].im += temp_im;
          i += iDelta2;
        }
        nd2 = chanStart + 1;
        for (ju = k; ju < nRowsD2; ju += k) {
          e = costab->data[ju];
          twid_im = sintab->data[ju];
          i = nd2;
          ihi = nd2 + iheight;
          while (i < ihi) {
            temp_re = e * y->data[i + iDelta].re - twid_im * y->data[i + iDelta]
              .im;
            temp_im = e * y->data[i + iDelta].im + twid_im * y->data[i + iDelta]
              .re;
            y->data[i + iDelta].re = y->data[i].re - temp_re;
            y->data[i + iDelta].im = y->data[i].im - temp_im;
            y->data[i].re += temp_re;
            y->data[i].im += temp_im;
            i += iDelta2;
          }
          nd2++;
        }
        k /= 2;
        iDelta = iDelta2;
        iDelta2 <<= 1;
        iheight -= iDelta;
      }
    }
    emxFree_real_T(&sintab);
    emxFree_real_T(&costab);
    if (y->size[0] > 1) {
      e = 1.0 / (real_T)y->size[0];
      ju = y->size[0] * y->size[1];
      emxEnsureCapacity((emxArray__common *)y, ju, (int32_T)sizeof(creal_T));
      nd2 = y->size[0];
      ju = y->size[1];
      nd2 *= ju;
      for (ju = 0; ju < nd2; ju++) {
        y->data[ju].re *= e;
        y->data[ju].im *= e;
      }
    }
  }
}

/* * */

void ifft(const emxArray_real_T *x, emxArray_creal_T *y)

{

int32_T ju;
int32_T minval;
emxArray_real_T *costab1q;
int32_T nd2;
int32_T ixDelta;
int32_T nRowsD2;
int32_T nRowsD4;
int32_T lastChan;
real_T e;
int32_T k;
emxArray_real_T *costab;
emxArray_real_T *sintab;
int32_T ix;
int32_T chanStart;
int32_T i;
boolean_T tst;
real_T temp_re;
real_T temp_im;
int32_T iDelta;
int32_T iDelta2;
int32_T iheight;
int32_T ihi;
real_T twid_im;
ju = y->size[0];
y->size[0] = x->size[0];
emxEnsureCapacity((emxArray__common *)y, ju, (int32_T)sizeof(creal_T));
if (x->size[0] == 0) {
} else {
  minval = x->size[0];
  b_emxInit_real_T(&costab1q, 2);
  nd2 = (x->size[0] - minval) + 1;
  if (1 >= nd2) {
    ixDelta = 1;
  } else {
    ixDelta = nd2;
  }
    ju = x->size[0];
    nRowsD2 = ju / 2;
    nRowsD4 = nRowsD2 / 2;
    lastChan = x->size[0] * (div_s32(x->size[0], x->size[0]) - 1);
    e = 6.2831853071795862 / (real_T)x->size[0];
    ju = costab1q->size[0] * costab1q->size[1];
    costab1q->size[0] = 1;
    costab1q->size[1] = nRowsD4 + 1;
    emxEnsureCapacity((emxArray__common *)costab1q, ju, (int32_T)sizeof(real_T));
    costab1q->data[0] = 1.0;
    nd2 = nRowsD4 / 2;
    for (k = 1; k <= nd2; k++) {
      costab1q->data[k] = cos(e * (real_T)k);
    }
    for (k = nd2 + 1; k < nRowsD4; k++) {
      costab1q->data[k] = sin(e * (real_T)(nRowsD4 - k));
    }
    b_emxInit_real_T(&costab, 2);
    b_emxInit_real_T(&sintab, 2);
    costab1q->data[nRowsD4] = 0.0;
    nd2 = (costab1q->size[1] - 1) << 1;
    ju = costab->size[0] * costab->size[1];
    costab->size[0] = 1;
    costab->size[1] = nd2 + 1;
    emxEnsureCapacity((emxArray__common *)costab, ju, (int32_T)sizeof(real_T));
    ju = sintab->size[0] * sintab->size[1];
    sintab->size[0] = 1;
    sintab->size[1] = nd2 + 1;
    emxEnsureCapacity((emxArray__common *)sintab, ju, (int32_T)sizeof(real_T));
    costab->data[0] = 1.0;
    sintab->data[0] = 0.0;
    for (k = 1; k < costab1q->size[1]; k++) {
      costab->data[k] = costab1q->data[k];
      sintab->data[k] = costab1q->data[(costab1q->size[1] - k) - 1];
    }
    for (k = costab1q->size[1]; k <= nd2; k++) {
      costab->data[k] = -costab1q->data[nd2 - k];
      sintab->data[k] = costab1q->data[(k - costab1q->size[1]) + 1];
    }
    emxFree_real_T(&costab1q);
    ix = 0;
    for (chanStart = 0; chanStart <= lastChan; chanStart += x->size[0]) {
      ju = 0;
      nd2 = chanStart;
      for (i = 1; i < minval; i++) {
        y->data[nd2].re = x->data[ix];
        y->data[nd2].im = 0.0;
        nd2 = x->size[0];
        tst = TRUE;
        while (tst) {
          nd2 >>= 1;
          ju ^= nd2;
          tst = ((ju & nd2) == 0);
        }
        nd2 = chanStart + ju;
        ix++;
      }
      y->data[nd2].re = x->data[ix];
      y->data[nd2].im = 0.0;
      ix += ixDelta;
      nd2 = (chanStart + x->size[0]) - 2;
      if (x->size[0] > 1) {
        for (i = chanStart; i <= nd2; i += 2) {
          temp_re = y->data[i + 1].re;
          temp_im = y->data[i + 1].im;
          y->data[i + 1].re = y->data[i].re - y->data[i + 1].re;
          y->data[i + 1].im = y->data[i].im - y->data[i + 1].im;
          y->data[i].re += temp_re;
          y->data[i].im += temp_im;
        }
      }
      iDelta = 2;
      iDelta2 = 4;
      k = nRowsD4;
      iheight = 1 + ((nRowsD4 - 1) << 2);
      while (k > 0) {
        i = chanStart;
        ihi = chanStart + iheight;
        while (i < ihi) {
          nd2 = i + iDelta;
          temp_re = y->data[nd2].re;
          temp_im = y->data[nd2].im;
          y->data[i + iDelta].re = y->data[i].re - y->data[nd2].re;
          y->data[i + iDelta].im = y->data[i].im - y->data[nd2].im;
          y->data[i].re += temp_re;
          y->data[i].im += temp_im;
          i += iDelta2;
        }
        nd2 = chanStart + 1;
        for (ju = k; ju < nRowsD2; ju += k) {
          e = costab->data[ju];
          twid_im = sintab->data[ju];
          i = nd2;
          ihi = nd2 + iheight;
          while (i < ihi) {
            temp_re = e * y->data[i + iDelta].re - twid_im * y->data[i + iDelta]
              .im;
            temp_im = e * y->data[i + iDelta].im + twid_im * y->data[i + iDelta]
              .re;
            y->data[i + iDelta].re = y->data[i].re - temp_re;
            y->data[i + iDelta].im = y->data[i].im - temp_im;
            y->data[i].re += temp_re;
            y->data[i].im += temp_im;
            i += iDelta2;
          }
          nd2++;
        }
        k /= 2;
        iDelta = iDelta2;
        iDelta2 <<= 1;
        iheight -= iDelta;
      }
    }
    emxFree_real_T(&sintab);
    emxFree_real_T(&costab);
    if (y->size[0] > 1) {
      e = 1.0 / (real_T)y->size[0];
      ju = y->size[0];
      emxEnsureCapacity((emxArray__common *)y, ju, (int32_T)sizeof(creal_T));
      nd2 = y->size[0];
      for (ju = 0; ju < nd2; ju++) {
        y->data[ju].re *= e;
        y->data[ju].im *= e;
      }
    }
  }
}

/* End of code generation (ifft.c) */

Thanks, Tang Laoya

4 Comments

laoya about 9 hours ago

Hi Sean de Wolski,

Thank you very much for your kindly reply.

I noticed that you have added 'nextpow2' in calling IFFT, however, the results are different from just placing [] at there. This function is from here:

http://www.biomecardio.com/matlab/dctn.m

so I can't modify the code.

Thanks, Tang Laoya

Sean de Wolski about 7 hours ago

If you download the same file from here:

http://www.mathworks.com/matlabcentral/fileexchange/26385-resize-any-arrays-and-images

Then you can make changes to it since it's covered under BSD.

laoya about 7 hours ago

Dear Sean de Wolski,

Thanks again for your kindly reply. Could yoy please give me more suggestion about how to resize the array to get almost the same result as the matlab code?

y = ifft(x,[],1);

What I mean that I can't modify the code is that that's a classical function and I have got reasonable results by matlab code using that function. However, if I modify the code to calculate the IFFT, the final results are totally different from what I got currently.

In fact, I have found another FFT source code by by Jens Joergen Nielsen (<http://www.jjj.de/fft/mixfft03.zip>) from here:

http://www.jjj.de/fft/fftpage.html

and write my own code as follows:

function y=ifft2(x)

x=conj(x);

% y=fft(x,[],1);

y=myfft(x,length(x)); % from mixfft03

y=conj(y)/length(x);

The results are the same as that from matlab. So it shows that it is a bug of matlab coder.

Thanks,

Tang Laoya

laoya

Contact us