API tools faq
paste
Advertisement
r4j

PLCA converted

r4j
Jul 15th, 2015
469
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 19.14 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.  * File: coupled_plca.c
  3.  *
  4.  * MATLAB Coder version            : 2.8
  5.  * C/C++ source code generated on  : 15-Jul-2015 06:42:26
  6.  */
  7.  
  8. /* Include Files */
  9. #include "rt_nonfinite.h"
  10. #include "coupled_plca.h"
  11. #include "coupled_plca_emxutil.h"
  12. #include "sum.h"
  13. #include "bsxfun.h"
  14. #include "diag.h"
  15. #include "conv2.h"
  16. #include "mrdivide.h"
  17. #include "rand.h"
  18. #include "randperm.h"
  19.  
  20. /* Function Declarations */
  21. static void eml_xgemm(int m, const emxArray_real_T *A, int lda, const double B
  22.                       [257275], emxArray_real_T *C, int ldc);
  23.  
  24. /* Function Definitions */
  25.  
  26. /*
  27.  * Arguments    : int m
  28.  *                const emxArray_real_T *A
  29.  *                int lda
  30.  *                const double B[257275]
  31.  *                emxArray_real_T *C
  32.  *                int ldc
  33.  * Return Type  : void
  34.  */
  35. static void eml_xgemm(int m, const emxArray_real_T *A, int lda, const double B
  36.                       [257275], emxArray_real_T *C, int ldc)
  37. {
  38.   int c;
  39.   int cr;
  40.   int i6;
  41.   int ic;
  42.   int br;
  43.   int ar;
  44.   int ib;
  45.   int ia;
  46.   if (m == 0) {
  47.   } else {
  48.     c = ldc * 250;
  49.     cr = 0;
  50.     while ((ldc > 0) && (cr <= c)) {
  51.       i6 = cr + m;
  52.       for (ic = cr; ic + 1 <= i6; ic++) {
  53.         C->data[ic] = 0.0;
  54.       }
  55.  
  56.       cr += ldc;
  57.     }
  58.  
  59.     br = 0;
  60.     cr = 0;
  61.     while ((ldc > 0) && (cr <= c)) {
  62.       ar = -1;
  63.       for (ib = br; ib + 1 <= br + 1025; ib++) {
  64.         if (B[ib] != 0.0) {
  65.           ia = ar;
  66.           i6 = cr + m;
  67.           for (ic = cr; ic + 1 <= i6; ic++) {
  68.             ia++;
  69.             C->data[ic] += B[ib] * A->data[ia];
  70.           }
  71.         }
  72.  
  73.         ar += lda;
  74.       }
  75.  
  76.       br += 1025;
  77.       cr += ldc;
  78.     }
  79.   }
  80. }
  81.  
  82. /*
  83.  * Perform coupled PLCA
  84.  *
  85.  *  [w zf zt h] = coupled_plca(spec_f, spec_t, Bs, Bt, R, iter)
  86.  *
  87.  *  spec_f: spectrogram with high spectral resolution
  88.  *  spec_t: spectrogram with high time resolution
  89.  *  Bs = blurring function for spectral components
  90.  *  Bt = blurring function for temporal components
  91.  *  R = number of components
  92.  *  iter = number of EM iterations
  93.  *
  94.  *  Output parameters
  95.  *  w = spectral basis vectors (hi-res)
  96.  *  zf = component weights for spectral vectors
  97.  *  zt = component weights for temporal vectors
  98.  *  h = temporal vectors (hi-res too!)
  99.  *
  100.  *  Juhan Nam
  101.  *    Apr-11-2010: initial version
  102.  *    Mar-11-2013: normalizing blur_h and blur_w right after blurring
  103.  * Arguments    : const double spec_f[257275]
  104.  *                const double spec_t[257275]
  105.  *                double Bs[257]
  106.  *                const double Bt[7]
  107.  *                double R
  108.  *                double iter
  109.  *                emxArray_real_T *w
  110.  *                emxArray_real_T *zf
  111.  *                emxArray_real_T *zt
  112.  *                emxArray_real_T *h
  113.  * Return Type  : void
  114.  */
  115. void coupled_plca(const double spec_f[257275], const double spec_t[257275],
  116.                   double Bs[257], const double Bt[7], double R, double iter,
  117.                   emxArray_real_T *w, emxArray_real_T *zf, emxArray_real_T *zt,
  118.                   emxArray_real_T *h)
  119. {
  120.   double ind[251];
  121.   int br;
  122.   int i;
  123.   int firstRowA;
  124.   emxArray_real_T *b_w;
  125.   int tmp_size[2];
  126.   double tmp_data[251];
  127.   emxArray_real_T *a;
  128.   emxArray_real_T *b_h;
  129.   emxArray_real_T *r0;
  130.   emxArray_real_T *b_zf;
  131.   emxArray_real_T *b_zt;
  132.   int it;
  133.   emxArray_real_T *blur_h;
  134.   emxArray_real_T *zfh;
  135.   emxArray_real_T *blur_w;
  136.   emxArray_real_T *nw;
  137.   emxArray_real_T *nh;
  138.   emxArray_real_T *C;
  139.   emxArray_real_T *b;
  140.   emxArray_real_T *b_a;
  141.   emxArray_real_T *b_blur_w;
  142.   emxArray_real_T *c_a;
  143.   emxArray_real_T *b_blur_h;
  144.   static double Rf[257275];
  145.   int ar;
  146.   int k;
  147.   int iB;
  148.   int ib;
  149.   int ia;
  150.   short unnamed_idx_1;
  151.   boolean_T b0;
  152.   int na;
  153.   int lastColB;
  154.   int lastColA;
  155.   int c;
  156.   int b_i;
  157.   int a_length;
  158.   double b_tmp_data[251];
  159.   int b_tmp_size[2];
  160.   static double y[257275];
  161.   unsigned int b_unnamed_idx_1;
  162.   static double b_spec_t[257275];
  163.   int nzf_size[2];
  164.   double nzf_data[251];
  165.   int nzt_size[2];
  166.   double nzt_data[251];
  167.   int c_tmp_size[2];
  168.   double a_data[251];
  169.   int a_size[1];
  170.   emxArray_real_T b_a_data;
  171.   double b_nzf_data[251];
  172.   int b_nzf_size[2];
  173.   double B;
  174.   double b_nzt_data[251];
  175.   int b_nzt_size[2];
  176.  
  177.   /*  size check */
  178.   /*  initializations */
  179.   /*  spectral basis vectors, */
  180.   /* w = rand(m, R)+1;  */
  181.   randperm(ind);
  182.   if (1.0 > R) {
  183.     br = 0;
  184.   } else {
  185.     br = (int)R;
  186.   }
  187.  
  188.   i = w->size[0] * w->size[1];
  189.   w->size[0] = 1025;
  190.   w->size[1] = br;
  191.   emxEnsureCapacity((emxArray__common *)w, i, (int)sizeof(double));
  192.   for (i = 0; i < br; i++) {
  193.     for (firstRowA = 0; firstRowA < 1025; firstRowA++) {
  194.       w->data[firstRowA + w->size[0] * i] = spec_f[firstRowA + 1025 * ((int)
  195.         ind[i] - 1)] + 2.2204460492503131E-16;
  196.     }
  197.   }
  198.  
  199.   emxInit_real_T(&b_w, 2);
  200.   sum(w, tmp_data, tmp_size);
  201.   i = b_w->size[0] * b_w->size[1];
  202.   b_w->size[0] = 1025;
  203.   b_w->size[1] = w->size[1];
  204.   emxEnsureCapacity((emxArray__common *)b_w, i, (int)sizeof(double));
  205.   br = w->size[0] * w->size[1];
  206.   for (i = 0; i < br; i++) {
  207.     b_w->data[i] = w->data[i];
  208.   }
  209.  
  210.   bsxfun(b_w, tmp_data, tmp_size, w);
  211.  
  212.   /*  temporal basis vectors */
  213.   c_rand(R, h);
  214.   i = h->size[0] * h->size[1];
  215.   h->size[1] = 251;
  216.   emxEnsureCapacity((emxArray__common *)h, i, (int)sizeof(double));
  217.   br = h->size[0];
  218.   firstRowA = h->size[1];
  219.   br *= firstRowA;
  220.   emxFree_real_T(&b_w);
  221.   for (i = 0; i < br; i++) {
  222.     h->data[i]++;
  223.   }
  224.  
  225.   b_emxInit_real_T(&a, 1);
  226.   emxInit_real_T(&b_h, 2);
  227.   b_sum(h, a);
  228.   i = b_h->size[0] * b_h->size[1];
  229.   b_h->size[0] = h->size[0];
  230.   b_h->size[1] = 251;
  231.   emxEnsureCapacity((emxArray__common *)b_h, i, (int)sizeof(double));
  232.   br = h->size[0] * h->size[1];
  233.   for (i = 0; i < br; i++) {
  234.     b_h->data[i] = h->data[i];
  235.   }
  236.  
  237.   b_bsxfun(b_h, a, h);
  238.   d_rand(R, zf);
  239.   i = zf->size[0] * zf->size[1];
  240.   emxEnsureCapacity((emxArray__common *)zf, i, (int)sizeof(double));
  241.   br = zf->size[0];
  242.   firstRowA = zf->size[1];
  243.   br *= firstRowA;
  244.   emxFree_real_T(&b_h);
  245.   for (i = 0; i < br; i++) {
  246.     zf->data[i]++;
  247.   }
  248.  
  249.   emxInit_real_T(&r0, 2);
  250.   emxInit_real_T(&b_zf, 2);
  251.  
  252.   /*  weights vector for good freq res plca */
  253.   c_sum(zf, r0);
  254.   i = b_zf->size[0] * b_zf->size[1];
  255.   b_zf->size[0] = zf->size[0];
  256.   b_zf->size[1] = zf->size[1];
  257.   emxEnsureCapacity((emxArray__common *)b_zf, i, (int)sizeof(double));
  258.   br = zf->size[0] * zf->size[1];
  259.   for (i = 0; i < br; i++) {
  260.     b_zf->data[i] = zf->data[i];
  261.   }
  262.  
  263.   mrdivide(b_zf, r0, zf);
  264.  
  265.   /*  normalize */
  266.   d_rand(R, zt);
  267.   i = zt->size[0] * zt->size[1];
  268.   emxEnsureCapacity((emxArray__common *)zt, i, (int)sizeof(double));
  269.   br = zt->size[0];
  270.   firstRowA = zt->size[1];
  271.   br *= firstRowA;
  272.   emxFree_real_T(&b_zf);
  273.   for (i = 0; i < br; i++) {
  274.     zt->data[i]++;
  275.   }
  276.  
  277.   emxInit_real_T(&b_zt, 2);
  278.  
  279.   /*  weights vector for good time res plca */
  280.   c_sum(zt, r0);
  281.   i = b_zt->size[0] * b_zt->size[1];
  282.   b_zt->size[0] = zt->size[0];
  283.   b_zt->size[1] = zt->size[1];
  284.   emxEnsureCapacity((emxArray__common *)b_zt, i, (int)sizeof(double));
  285.   br = zt->size[0] * zt->size[1];
  286.   for (i = 0; i < br; i++) {
  287.     b_zt->data[i] = zt->data[i];
  288.   }
  289.  
  290.   mrdivide(b_zt, r0, zt);
  291.  
  292.   /*  normalize */
  293.   emxFree_real_T(&b_zt);
  294.   emxFree_real_T(&r0);
  295.  
  296.   /*  vertical */
  297.   /*  horizontal */
  298.   /*  EM iterations */
  299.   it = 0;
  300.   emxInit_real_T(&blur_h, 2);
  301.   emxInit_real_T(&zfh, 2);
  302.   emxInit_real_T(&blur_w, 2);
  303.   emxInit_real_T(&nw, 2);
  304.   emxInit_real_T(&nh, 2);
  305.   emxInit_real_T(&C, 2);
  306.   emxInit_real_T(&b, 2);
  307.   emxInit_real_T(&b_a, 2);
  308.   emxInit_real_T(&b_blur_w, 2);
  309.   b_emxInit_real_T(&c_a, 1);
  310.   emxInit_real_T(&b_blur_h, 2);
  311.   while (it <= (int)iter - 1) {
  312.     /*  E-step for high frequency resolution spectrogram */
  313.     conv2(h, Bt, blur_h);
  314.     b_sum(blur_h, a);
  315.     i = c_a->size[0];
  316.     c_a->size[0] = a->size[0];
  317.     emxEnsureCapacity((emxArray__common *)c_a, i, (int)sizeof(double));
  318.     br = a->size[0];
  319.     for (i = 0; i < br; i++) {
  320.       c_a->data[i] = a->data[i] + 2.2204460492503131E-16;
  321.     }
  322.  
  323.     i = b_blur_h->size[0] * b_blur_h->size[1];
  324.     b_blur_h->size[0] = blur_h->size[0];
  325.     b_blur_h->size[1] = 251;
  326.     emxEnsureCapacity((emxArray__common *)b_blur_h, i, (int)sizeof(double));
  327.     br = blur_h->size[0] * blur_h->size[1];
  328.     for (i = 0; i < br; i++) {
  329.       b_blur_h->data[i] = blur_h->data[i];
  330.     }
  331.  
  332.     b_bsxfun(b_blur_h, c_a, blur_h);
  333.     diag(zf, a);
  334.     i = zfh->size[0] * zfh->size[1];
  335.     zfh->size[0] = a->size[0];
  336.     zfh->size[1] = 251;
  337.     emxEnsureCapacity((emxArray__common *)zfh, i, (int)sizeof(double));
  338.     br = a->size[0];
  339.     for (i = 0; i < br; i++) {
  340.       for (firstRowA = 0; firstRowA < 251; firstRowA++) {
  341.         zfh->data[i + zfh->size[0] * firstRowA] = a->data[i] * blur_h->
  342.           data[firstRowA];
  343.       }
  344.     }
  345.  
  346.     if ((w->size[1] == 1) || (zfh->size[0] == 1)) {
  347.       for (i = 0; i < 1025; i++) {
  348.         for (firstRowA = 0; firstRowA < 251; firstRowA++) {
  349.           Rf[i + 1025 * firstRowA] = 0.0;
  350.           br = w->size[1];
  351.           for (ar = 0; ar < br; ar++) {
  352.             Rf[i + 1025 * firstRowA] += w->data[i + w->size[0] * ar] * zfh->
  353.               data[ar + zfh->size[0] * firstRowA];
  354.           }
  355.         }
  356.       }
  357.     } else {
  358.       k = w->size[1];
  359.       memset(&Rf[0], 0, 257275U * sizeof(double));
  360.       for (firstRowA = 0; firstRowA < 256252; firstRowA += 1025) {
  361.         for (iB = firstRowA; iB + 1 <= firstRowA + 1025; iB++) {
  362.           Rf[iB] = 0.0;
  363.         }
  364.       }
  365.  
  366.       br = 0;
  367.       for (firstRowA = 0; firstRowA < 256252; firstRowA += 1025) {
  368.         ar = 0;
  369.         i = br + k;
  370.         for (ib = br; ib + 1 <= i; ib++) {
  371.           if (zfh->data[ib] != 0.0) {
  372.             ia = ar;
  373.             for (iB = firstRowA; iB + 1 <= firstRowA + 1025; iB++) {
  374.               ia++;
  375.               Rf[iB] += zfh->data[ib] * w->data[ia - 1];
  376.             }
  377.           }
  378.  
  379.           ar += 1025;
  380.         }
  381.  
  382.         br += k;
  383.       }
  384.     }
  385.  
  386.     for (i = 0; i < 257275; i++) {
  387.       Rf[i] = spec_f[i] / Rf[i];
  388.     }
  389.  
  390.     /*  E for high time resolution spectrogram */
  391.     diag(zt, a);
  392.     i = blur_h->size[0] * blur_h->size[1];
  393.     blur_h->size[0] = a->size[0];
  394.     blur_h->size[1] = 251;
  395.     emxEnsureCapacity((emxArray__common *)blur_h, i, (int)sizeof(double));
  396.     br = a->size[0];
  397.     for (i = 0; i < br; i++) {
  398.       for (firstRowA = 0; firstRowA < 251; firstRowA++) {
  399.         blur_h->data[i + blur_h->size[0] * firstRowA] = a->data[i] * h->
  400.           data[firstRowA];
  401.       }
  402.     }
  403.  
  404.     unnamed_idx_1 = (short)w->size[1];
  405.     i = blur_w->size[0] * blur_w->size[1];
  406.     blur_w->size[0] = 1025;
  407.     emxEnsureCapacity((emxArray__common *)blur_w, i, (int)sizeof(double));
  408.     i = blur_w->size[0] * blur_w->size[1];
  409.     blur_w->size[1] = unnamed_idx_1;
  410.     emxEnsureCapacity((emxArray__common *)blur_w, i, (int)sizeof(double));
  411.     br = 1025 * unnamed_idx_1;
  412.     for (i = 0; i < br; i++) {
  413.       blur_w->data[i] = 0.0;
  414.     }
  415.  
  416.     if ((w->size[1] == 0) || (unnamed_idx_1 == 0)) {
  417.       b0 = true;
  418.     } else {
  419.       b0 = false;
  420.     }
  421.  
  422.     if (!b0) {
  423.       na = w->size[1];
  424.       if (1 <= unnamed_idx_1 - 1) {
  425.         lastColB = 1;
  426.       } else {
  427.         lastColB = unnamed_idx_1;
  428.       }
  429.  
  430.       br = 0;
  431.       while (br <= lastColB - 1) {
  432.         if (na - 1 < unnamed_idx_1 - 1) {
  433.           lastColA = na;
  434.         } else {
  435.           lastColA = unnamed_idx_1;
  436.         }
  437.  
  438.         for (k = 0; k < lastColA; k++) {
  439.           if (k >= 0) {
  440.             firstRowA = k;
  441.           } else {
  442.             firstRowA = 0;
  443.           }
  444.  
  445.           ia = firstRowA * 1025;
  446.           c = k * 1025;
  447.           iB = 0;
  448.           for (i = 0; i < 257; i++) {
  449.             if (i < 128) {
  450.               firstRowA = 128 - i;
  451.             } else {
  452.               firstRowA = 0;
  453.             }
  454.  
  455.             if (i + 1025 <= 1152) {
  456.               b_i = 1025;
  457.             } else {
  458.               b_i = 1153 - i;
  459.             }
  460.  
  461.             a_length = b_i - firstRowA;
  462.             br = c + firstRowA;
  463.             ar = ia;
  464.             for (ib = 1; ib <= a_length; ib++) {
  465.               blur_w->data[ar] += Bs[iB] * w->data[br];
  466.               br++;
  467.               ar++;
  468.             }
  469.  
  470.             iB++;
  471.             if (i >= 128) {
  472.               ia++;
  473.             }
  474.           }
  475.         }
  476.  
  477.         br = 1;
  478.       }
  479.     }
  480.  
  481.     sum(blur_w, tmp_data, tmp_size);
  482.     b_tmp_size[0] = 1;
  483.     b_tmp_size[1] = tmp_size[1];
  484.     br = tmp_size[0] * tmp_size[1];
  485.     for (i = 0; i < br; i++) {
  486.       b_tmp_data[i] = tmp_data[i] + 2.2204460492503131E-16;
  487.     }
  488.  
  489.     i = b_blur_w->size[0] * b_blur_w->size[1];
  490.     b_blur_w->size[0] = 1025;
  491.     b_blur_w->size[1] = blur_w->size[1];
  492.     emxEnsureCapacity((emxArray__common *)b_blur_w, i, (int)sizeof(double));
  493.     br = blur_w->size[0] * blur_w->size[1];
  494.     for (i = 0; i < br; i++) {
  495.       b_blur_w->data[i] = blur_w->data[i];
  496.     }
  497.  
  498.     bsxfun(b_blur_w, b_tmp_data, b_tmp_size, blur_w);
  499.     if ((blur_w->size[1] == 1) || (blur_h->size[0] == 1)) {
  500.       for (i = 0; i < 1025; i++) {
  501.         for (firstRowA = 0; firstRowA < 251; firstRowA++) {
  502.           y[i + 1025 * firstRowA] = 0.0;
  503.           br = blur_w->size[1];
  504.           for (ar = 0; ar < br; ar++) {
  505.             y[i + 1025 * firstRowA] += blur_w->data[i + blur_w->size[0] * ar] *
  506.               blur_h->data[ar + blur_h->size[0] * firstRowA];
  507.           }
  508.         }
  509.       }
  510.     } else {
  511.       k = blur_w->size[1];
  512.       memset(&y[0], 0, 257275U * sizeof(double));
  513.       for (firstRowA = 0; firstRowA < 256252; firstRowA += 1025) {
  514.         for (iB = firstRowA; iB + 1 <= firstRowA + 1025; iB++) {
  515.           y[iB] = 0.0;
  516.         }
  517.       }
  518.  
  519.       br = 0;
  520.       for (firstRowA = 0; firstRowA < 256252; firstRowA += 1025) {
  521.         ar = 0;
  522.         i = br + k;
  523.         for (ib = br; ib + 1 <= i; ib++) {
  524.           if (blur_h->data[ib] != 0.0) {
  525.             ia = ar;
  526.             for (iB = firstRowA; iB + 1 <= firstRowA + 1025; iB++) {
  527.               ia++;
  528.               y[iB] += blur_h->data[ib] * blur_w->data[ia - 1];
  529.             }
  530.           }
  531.  
  532.           ar += 1025;
  533.         }
  534.  
  535.         br += k;
  536.       }
  537.     }
  538.  
  539.     /*  M-step */
  540.     i = b->size[0] * b->size[1];
  541.     b->size[0] = 251;
  542.     b->size[1] = zfh->size[0];
  543.     emxEnsureCapacity((emxArray__common *)b, i, (int)sizeof(double));
  544.     br = zfh->size[0];
  545.     for (i = 0; i < br; i++) {
  546.       for (firstRowA = 0; firstRowA < 251; firstRowA++) {
  547.         b->data[firstRowA + b->size[0] * i] = zfh->data[i + zfh->size[0] *
  548.           firstRowA];
  549.       }
  550.     }
  551.  
  552.     b_unnamed_idx_1 = (unsigned int)b->size[1];
  553.     i = C->size[0] * C->size[1];
  554.     C->size[0] = 1025;
  555.     emxEnsureCapacity((emxArray__common *)C, i, (int)sizeof(double));
  556.     i = C->size[0] * C->size[1];
  557.     C->size[1] = (int)b_unnamed_idx_1;
  558.     emxEnsureCapacity((emxArray__common *)C, i, (int)sizeof(double));
  559.     br = 1025 * (int)b_unnamed_idx_1;
  560.     for (i = 0; i < br; i++) {
  561.       C->data[i] = 0.0;
  562.     }
  563.  
  564.     if (b->size[1] == 0) {
  565.     } else {
  566.       c = 1025 * (b->size[1] - 1);
  567.       for (firstRowA = 0; firstRowA <= c; firstRowA += 1025) {
  568.         for (iB = firstRowA; iB + 1 <= firstRowA + 1025; iB++) {
  569.           C->data[iB] = 0.0;
  570.         }
  571.       }
  572.  
  573.       br = 0;
  574.       for (firstRowA = 0; firstRowA <= c; firstRowA += 1025) {
  575.         ar = 0;
  576.         for (ib = br; ib + 1 <= br + 251; ib++) {
  577.           if (b->data[ib] != 0.0) {
  578.             ia = ar;
  579.             for (iB = firstRowA; iB + 1 <= firstRowA + 1025; iB++) {
  580.               ia++;
  581.               C->data[iB] += b->data[ib] * Rf[ia - 1];
  582.             }
  583.           }
  584.  
  585.           ar += 1025;
  586.         }
  587.  
  588.         br += 251;
  589.       }
  590.     }
  591.  
  592.     i = nw->size[0] * nw->size[1];
  593.     nw->size[0] = 1025;
  594.     nw->size[1] = w->size[1];
  595.     emxEnsureCapacity((emxArray__common *)nw, i, (int)sizeof(double));
  596.     br = w->size[0] * w->size[1];
  597.     for (i = 0; i < br; i++) {
  598.       nw->data[i] = w->data[i] * C->data[i];
  599.     }
  600.  
  601.     i = b_a->size[0] * b_a->size[1];
  602.     b_a->size[0] = w->size[1];
  603.     b_a->size[1] = 1025;
  604.     emxEnsureCapacity((emxArray__common *)b_a, i, (int)sizeof(double));
  605.     for (i = 0; i < 1025; i++) {
  606.       br = w->size[1];
  607.       for (firstRowA = 0; firstRowA < br; firstRowA++) {
  608.         b_a->data[firstRowA + b_a->size[0] * i] = w->data[i + w->size[0] *
  609.           firstRowA];
  610.       }
  611.     }
  612.  
  613.     b_unnamed_idx_1 = (unsigned int)b_a->size[0];
  614.     i = nh->size[0] * nh->size[1];
  615.     nh->size[0] = (int)b_unnamed_idx_1;
  616.     nh->size[1] = 251;
  617.     emxEnsureCapacity((emxArray__common *)nh, i, (int)sizeof(double));
  618.     br = (int)b_unnamed_idx_1 * 251;
  619.     for (i = 0; i < br; i++) {
  620.       nh->data[i] = 0.0;
  621.     }
  622.  
  623.     for (i = 0; i < 257275; i++) {
  624.       b_spec_t[i] = spec_t[i] / y[i];
  625.     }
  626.  
  627.     eml_xgemm(b_a->size[0], b_a, b_a->size[0], b_spec_t, nh, b_a->size[0]);
  628.     i = nh->size[0] * nh->size[1];
  629.     nh->size[0] = blur_h->size[0];
  630.     nh->size[1] = 251;
  631.     emxEnsureCapacity((emxArray__common *)nh, i, (int)sizeof(double));
  632.     br = blur_h->size[0] * blur_h->size[1];
  633.     for (i = 0; i < br; i++) {
  634.       nh->data[i] *= blur_h->data[i];
  635.     }
  636.  
  637.     sum(nw, nzf_data, nzf_size);
  638.     sum(blur_w, nzt_data, nzt_size);
  639.  
  640.     /*  normalize */
  641.     sum(nw, tmp_data, tmp_size);
  642.     c_tmp_size[0] = 1;
  643.     c_tmp_size[1] = tmp_size[1];
  644.     br = tmp_size[0] * tmp_size[1];
  645.     for (i = 0; i < br; i++) {
  646.       b_tmp_data[i] = tmp_data[i] + 2.2204460492503131E-16;
  647.     }
  648.  
  649.     bsxfun(nw, b_tmp_data, c_tmp_size, w);
  650.     b_sum(nh, a);
  651.     a_size[0] = a->size[0];
  652.     br = a->size[0];
  653.     for (i = 0; i < br; i++) {
  654.       a_data[i] = a->data[i] + 2.2204460492503131E-16;
  655.     }
  656.  
  657.     b_a_data.data = (double *)&a_data;
  658.     b_a_data.size = (int *)&a_size;
  659.     b_a_data.allocatedSize = 251;
  660.     b_a_data.numDimensions = 1;
  661.     b_a_data.canFreeData = false;
  662.     b_bsxfun(nh, &b_a_data, h);
  663.     b_nzf_size[0] = 1;
  664.     b_nzf_size[1] = nzf_size[1];
  665.     br = nzf_size[0] * nzf_size[1];
  666.     for (i = 0; i < br; i++) {
  667.       b_nzf_data[i] = nzf_data[i] + 2.2204460492503131E-16;
  668.     }
  669.  
  670.     B = d_sum(b_nzf_data, b_nzf_size);
  671.     i = zf->size[0] * zf->size[1];
  672.     zf->size[0] = 1;
  673.     zf->size[1] = nzf_size[1];
  674.     emxEnsureCapacity((emxArray__common *)zf, i, (int)sizeof(double));
  675.     br = nzf_size[0] * nzf_size[1];
  676.     for (i = 0; i < br; i++) {
  677.       zf->data[i] = nzf_data[i] / B;
  678.     }
  679.  
  680.     b_nzt_size[0] = 1;
  681.     b_nzt_size[1] = nzt_size[1];
  682.     br = nzt_size[0] * nzt_size[1];
  683.     for (i = 0; i < br; i++) {
  684.       b_nzt_data[i] = nzt_data[i] + 2.2204460492503131E-16;
  685.     }
  686.  
  687.     B = d_sum(b_nzt_data, b_nzt_size);
  688.     i = zt->size[0] * zt->size[1];
  689.     zt->size[0] = 1;
  690.     zt->size[1] = nzt_size[1];
  691.     emxEnsureCapacity((emxArray__common *)zt, i, (int)sizeof(double));
  692.     br = nzt_size[0] * nzt_size[1];
  693.     for (i = 0; i < br; i++) {
  694.       zt->data[i] = nzt_data[i] / B;
  695.     }
  696.  
  697.     it++;
  698.   }
  699.  
  700.   emxFree_real_T(&b_blur_h);
  701.   emxFree_real_T(&c_a);
  702.   emxFree_real_T(&b_blur_w);
  703.   emxFree_real_T(&b_a);
  704.   emxFree_real_T(&b);
  705.   emxFree_real_T(&a);
  706.   emxFree_real_T(&C);
  707.   emxFree_real_T(&nh);
  708.   emxFree_real_T(&nw);
  709.   emxFree_real_T(&blur_w);
  710.   emxFree_real_T(&zfh);
  711.   emxFree_real_T(&blur_h);
  712. }
  713.  
  714. /*
  715.  * File trailer for coupled_plca.c
  716.  *
  717.  * [EOF]
  718.  */
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!