File: /Users/peter/svn/azumio/trunk/iOSCommon/heartrate/lib/hr_analyzer.c

    1   /*
    2    * hr_analyzer.c
    3    *
    4    * Code generation for function 'hr_analyzer'
    5    *
    6    * C source code generated on: Sun Jan 27 20:19:03 2013
    7    *
    8    */
    9   
   10   /* Include files */
   11   #include "rt_nonfinite.h"
   12   #include "hr_analyzer.h"
   13   
   14   /* Type Definitions */
   15   
   16   /* Named Constants */
   17   
   18   /* Variable Declarations */
   19   
   20   /* Variable Definitions */
   21   
   22   /* Function Declarations */
   23   static void b_eml_xaxpy(real_T a, const real_T x[41], real_T y[41]);
   24   static void b_filter(const real_T x_data[10000], const int32_T x_sizes[1],
   25                        real_T y_data[50000], int32_T y_sizes[1]);
   26   static void b_max(const real_T varargin_1_data[50000], const int32_T
   27                     varargin_1_sizes[1], const real_T varargin_2_data[50000],
   28                     const int32_T varargin_2_sizes[1], real_T maxval_data[50000],
   29                     int32_T maxval_sizes[1]);
   30   static void c_eml_xaxpy(real_T a, const real_T x[41], real_T y[41]);
   31   static void eml_scalexp_alloc(const real_T varargin_1_data[50000], const int32_T
   32     varargin_1_sizes[1], const real_T varargin_2_data[50000], const int32_T
   33     varargin_2_sizes[1], real_T z_data[50000], int32_T z_sizes[1]);
   34   static void eml_xaxpy(int32_T n, real_T a, const real_T x_data[50000], const
   35                         int32_T x_sizes[1], int32_T ix0, real_T y_data[50000],
   36                         int32_T y_sizes[1], int32_T iy0);
   37   static void filter(const real_T x_data[50000], const int32_T x_sizes[1], real_T
   38                      y_data[50000], int32_T y_sizes[1]);
   39   static void flipud(real_T x_data[10000], int32_T x_sizes[1]);
   40   
   41   /* Function Definitions */
   42   
   43   /*
   44    *
   45    */
   46   static void b_eml_xaxpy(real_T a, const real_T x[41], real_T y[41])
   47   {
   48     int32_T ix;
   49     int32_T iy;
   50     int32_T k;
   51     if (a == 0.0) {
   52     } else {
   53       ix = 0;
   54       iy = 0;
   55       for (k = 0; k < 41; k++) {
   56         y[iy] += a * x[ix];
   57         ix++;
   58         iy++;
   59       }
   60     }
   61   }
   62   
   63   /*
   64    *
   65    */
   66   static void b_filter(const real_T x_data[10000], const int32_T x_sizes[1],
   67                        real_T y_data[50000], int32_T y_sizes[1])
   68   {
   69     int32_T j;
   70     int32_T k;
   71     static const real_T dv0[41] = { 0.0023611273282788886, 0.0026754160612104489,
   72       0.0033736336905303569, 0.0044948250919669488, 0.0060621018954461151,
   73       0.00808080523183687, 0.010537404054552698, 0.013399198470849851,
   74       0.016614861289538603, 0.020115812614132048, 0.0238183836076616,
   75       0.027626688387000657, 0.031436089173432075, 0.035137110989079567,
   76       0.038619639773400868, 0.041777222926969913, 0.044511284739703484,
   77       0.046735071297835565, 0.048377150240273885, 0.049384309684601189,
   78       0.049723726903396541, 0.049384309684601189, 0.048377150240273885,
   79       0.046735071297835565, 0.044511284739703484, 0.041777222926969913,
   80       0.038619639773400868, 0.035137110989079567, 0.031436089173432075,
   81       0.027626688387000657, 0.0238183836076616, 0.020115812614132048,
   82       0.016614861289538603, 0.013399198470849851, 0.010537404054552698,
   83       0.00808080523183687, 0.0060621018954461151, 0.0044948250919669488,
   84       0.0033736336905303569, 0.0026754160612104489, 0.0023611273282788886 };
   85   
   86     real_T dbuffer[41];
   87     static const real_T b[41] = { 0.0023611273282788886, 0.0026754160612104489,
   88       0.0033736336905303569, 0.0044948250919669488, 0.0060621018954461151,
   89       0.00808080523183687, 0.010537404054552698, 0.013399198470849851,
   90       0.016614861289538603, 0.020115812614132048, 0.0238183836076616,
   91       0.027626688387000657, 0.031436089173432075, 0.035137110989079567,
   92       0.038619639773400868, 0.041777222926969913, 0.044511284739703484,
   93       0.046735071297835565, 0.048377150240273885, 0.049384309684601189,
   94       0.049723726903396541, 0.049384309684601189, 0.048377150240273885,
   95       0.046735071297835565, 0.044511284739703484, 0.041777222926969913,
   96       0.038619639773400868, 0.035137110989079567, 0.031436089173432075,
   97       0.027626688387000657, 0.0238183836076616, 0.020115812614132048,
   98       0.016614861289538603, 0.013399198470849851, 0.010537404054552698,
   99       0.00808080523183687, 0.0060621018954461151, 0.0044948250919669488,
  100       0.0033736336905303569, 0.0026754160612104489, 0.0023611273282788886 };
  101   
  102     y_sizes[0] = (int16_T)x_sizes[0];
  103     if (x_sizes[0] >= 82) {
  104       j = y_sizes[0];
  105       j--;
  106       for (k = 0; k <= j; k++) {
  107         y_data[k] = 0.0;
  108       }
  109   
  110       for (k = 0; k < 41; k++) {
  111         eml_xaxpy(x_sizes[0] - k, dv0[k], x_data, x_sizes, 1, y_data, y_sizes, k +
  112                   1);
  113       }
  114     } else {
  115       memset((void *)&dbuffer[1], 0, 40U * sizeof(real_T));
  116       for (j = 0; j + 1 <= x_sizes[0]; j++) {
  117         for (k = 0; k < 40; k++) {
  118           dbuffer[k] = dbuffer[k + 1];
  119         }
  120   
  121         dbuffer[40] = 0.0;
  122         if (x_data[j] == 0.0) {
  123           for (k = 0; k < 41; k++) {
  124             dbuffer[k] += x_data[j] * dv0[k];
  125           }
  126         } else {
  127           c_eml_xaxpy(x_data[j], b, dbuffer);
  128         }
  129   
  130         y_data[j] = dbuffer[0];
  131       }
  132     }
  133   }
  134   
  135   /*
  136    *
  137    */
  138   static void b_max(const real_T varargin_1_data[50000], const int32_T
  139                     varargin_1_sizes[1], const real_T varargin_2_data[50000],
  140                     const int32_T varargin_2_sizes[1], real_T maxval_data[50000],
  141                     int32_T maxval_sizes[1])
  142   {
  143     int32_T k;
  144     real_T u0;
  145     real_T u1;
  146     eml_scalexp_alloc(varargin_1_data, varargin_1_sizes, varargin_2_data,
  147                       varargin_2_sizes, maxval_data, maxval_sizes);
  148     for (k = 0; k + 1 <= maxval_sizes[0]; k++) {
  149       u0 = varargin_1_data[k];
  150       u1 = varargin_2_data[k];
  151       u0 = (u0 >= u1) || rtIsNaN(u1) ? u0 : u1;
  152       maxval_data[k] = u0;
  153     }
  154   }
  155   
  156   /*
  157    *
  158    */
  159   static void c_eml_xaxpy(real_T a, const real_T x[41], real_T y[41])
  160   {
  161     int32_T ix;
  162     int32_T iy;
  163     int32_T k;
  164     if (a == 0.0) {
  165     } else {
  166       ix = 0;
  167       iy = 0;
  168       for (k = 0; k < 41; k++) {
  169         y[iy] += a * x[ix];
  170         ix++;
  171         iy++;
  172       }
  173     }
  174   }
  175   
  176   /*
  177    *
  178    */
  179   static void eml_scalexp_alloc(const real_T varargin_1_data[50000], const int32_T
  180     varargin_1_sizes[1], const real_T varargin_2_data[50000], const int32_T
  181     varargin_2_sizes[1], real_T z_data[50000], int32_T z_sizes[1])
  182   {
  183     z_sizes[0] = (uint16_T)varargin_1_sizes[0];
  184   }
  185   
  186   /*
  187    *
  188    */
  189   static void eml_xaxpy(int32_T n, real_T a, const real_T x_data[50000], const
  190                         int32_T x_sizes[1], int32_T ix0, real_T y_data[50000],
  191                         int32_T y_sizes[1], int32_T iy0)
  192   {
  193     int32_T ix;
  194     int32_T iy;
  195     int32_T loop_ub;
  196     int32_T k;
  197     if (a == 0.0) {
  198     } else {
  199       ix = ix0 - 1;
  200       iy = iy0 - 1;
  201       loop_ub = n - 1;
  202       for (k = 0; k <= loop_ub; k++) {
  203         y_data[iy] += a * x_data[ix];
  204         ix++;
  205         iy++;
  206       }
  207     }
  208   }
  209   
  210   /*
  211    *
  212    */
  213   static void filter(const real_T x_data[50000], const int32_T x_sizes[1], real_T
  214                      y_data[50000], int32_T y_sizes[1])
  215   {
  216     int32_T j;
  217     int32_T k;
  218     static const real_T b[41] = { -0.0069687082408333356, -0.0098234327352945336,
  219       -0.015428282490856905, -0.022600435931048331, -0.028736176949146989,
  220       -0.030787403541533767, -0.026899242824176162, -0.017934742285145765,
  221       -0.0079305046999676911, -0.0028724904614616849, -0.00797790963137552,
  222       -0.024556986233664251, -0.048022079814209341, -0.068358785745496475,
  223       -0.073375407143698765, -0.053710110915723232, -0.0075439371231891131,
  224       0.057191162272427772, 0.12442662372882912, 0.17510107665826793,
  225       0.19393482504938614, 0.17510107665826793, 0.12442662372882912,
  226       0.057191162272427772, -0.0075439371231891131, -0.053710110915723232,
  227       -0.073375407143698765, -0.068358785745496475, -0.048022079814209341,
  228       -0.024556986233664251, -0.00797790963137552, -0.0028724904614616849,
  229       -0.0079305046999676911, -0.017934742285145765, -0.026899242824176162,
  230       -0.030787403541533767, -0.028736176949146989, -0.022600435931048331,
  231       -0.015428282490856905, -0.0098234327352945336, -0.0069687082408333356 };
  232   
  233     int32_T jend;
  234     real_T dbuffer[41];
  235     y_sizes[0] = (uint16_T)x_sizes[0];
  236     if (x_sizes[0] >= 82) {
  237       j = y_sizes[0];
  238       j--;
  239       for (k = 0; k <= j; k++) {
  240         y_data[k] = 0.0;
  241       }
  242   
  243       for (k = 0; k < 41; k++) {
  244         if (b[k] == 0.0) {
  245           jend = (k + x_sizes[0]) - k;
  246           for (j = k; j + 1 <= jend; j++) {
  247             y_data[j] += b[k] * x_data[j - k];
  248           }
  249         } else {
  250           eml_xaxpy(x_sizes[0] - k, b[k], x_data, x_sizes, 1, y_data, y_sizes, k +
  251                     1);
  252         }
  253       }
  254     } else {
  255       memset((void *)&dbuffer[1], 0, 40U * sizeof(real_T));
  256       for (j = 0; j + 1 <= x_sizes[0]; j++) {
  257         for (k = 0; k < 40; k++) {
  258           dbuffer[k] = dbuffer[k + 1];
  259         }
  260   
  261         dbuffer[40] = 0.0;
  262         if (x_data[j] == 0.0) {
  263           for (k = 0; k < 41; k++) {
  264             dbuffer[k] += x_data[j] * b[k];
  265           }
  266         } else {
  267           b_eml_xaxpy(x_data[j], b, dbuffer);
  268         }
  269   
  270         y_data[j] = dbuffer[0];
  271       }
  272     }
  273   }
  274   
  275   /*
  276    *
  277    */
  278   static void flipud(real_T x_data[10000], int32_T x_sizes[1])
  279   {
  280     int32_T m;
  281     int32_T md2;
  282     int32_T i;
  283     real_T xtmp;
  284     m = x_sizes[0];
  285     md2 = m / 2;
  286     for (i = 1; i <= md2; i++) {
  287       xtmp = x_data[i - 1];
  288       x_data[i - 1] = x_data[m - i];
  289       x_data[m - i] = xtmp;
  290     }
  291   }
  292   
  293   /*
  294    * function [ value ] = filterForDispaly( value )
  295    */
  296   void filterForDisplay(real_T value_data[50000], int32_T value_sizes[1])
  297   {
  298     int32_T ac_sizes;
  299     static real_T ac_data[50000];
  300     int32_T loop_ub;
  301     int32_T i0;
  302     static real_T b_ac_data[50000];
  303     int32_T tmp_sizes;
  304     static real_T tmp_data[50000];
  305     real_T b_tmp_data[10000];
  306     int32_T acnorm_div_sizes;
  307     real_T acnorm_div_data[10000];
  308   
  309     /* 'filterForDisplay:4' if length(value)>50000 */
  310     /* 'filterForDisplay:12' fs = 30; */
  311     /* 'filterForDisplay:13' hrmin = 30/60; */
  312     /* 'filterForDisplay:16' p1rows = 5; */
  313     /* p = dcblock(hrmin/2,30);              % get filter coefficient */
  314     /* b = [1 -1];                         % set up differentiator */
  315     /* a = [1 -p];                         % set up integrator */
  316     /*  ac = filter(b,a, value(1:cnt,1) );     */
  317     /* 'filterForDisplay:26' fmin = 0.6; */
  318     /* 'filterForDisplay:27' fmax = 3; */
  319     /* onpassfir = fir1(60,[fmin./fs.*2 fmax./fs.*2]); */
  320     /* onpassfir = fir1(10,fmax./fs.*2,'low'); */
  321     /*  ac = filter(onpassfir,1,ac); */
  322     /* 'filterForDisplay:32' filt = pulse_bandpass_const(); */
  323     /* 'filterForDisplay:34' ac = filter(filt,1,value); */
  324     filter(value_data, value_sizes, ac_data, &ac_sizes);
  325   
  326     /* remove more DC */
  327     /*       p = dcblock(hrmin,30);              % get filter coefficient */
  328     /*      b = [1 -1];                         % set up differentiator */
  329     /*      a = [1 -p];                         % set up integrator */
  330     /*      acenv = filter(b,a, value );  */
  331     /* 'filterForDisplay:44' acenv = ac; */
  332     /* 'filterForDisplay:47' envFilterWidth = 40; */
  333     /* 'filterForDisplay:49' coder.varsize('envF',[1,41],[false,false]) */
  334     /* 'filterForDisplay:50' coder.varsize('envT',[10000,1],[true,false]) */
  335     /* 'filterForDisplay:52' envF = fir1(envFilterWidth,0.4/(fs/2),'low'); */
  336     /* do a top enveloper */
  337     /* 'filterForDisplay:57' envT = max(acenv,-acenv); */
  338     loop_ub = ac_sizes - 1;
  339     for (i0 = 0; i0 <= loop_ub; i0++) {
  340       b_ac_data[i0] = -ac_data[i0];
  341     }
  342   
  343     b_max(ac_data, *(int32_T (*)[1])&ac_sizes, b_ac_data, *(int32_T (*)[1])&
  344           ac_sizes, tmp_data, &tmp_sizes);
  345   
  346     /* envT = filtfilt(envF,1,envT.*2); */
  347     /* fake filtfilt */
  348     /* 'filterForDisplay:62' envT = filter(envF,1, envT); */
  349     loop_ub = tmp_sizes - 1;
  350     for (i0 = 0; i0 <= loop_ub; i0++) {
  351       b_tmp_data[i0] = tmp_data[i0];
  352     }
  353   
  354     b_filter(b_tmp_data, *(int32_T (*)[1])&tmp_sizes, tmp_data, &tmp_sizes);
  355   
  356     /* 'filterForDisplay:63' envT = filter(envF,1, flipud(envT)); */
  357     acnorm_div_sizes = tmp_sizes;
  358     loop_ub = tmp_sizes - 1;
  359     for (i0 = 0; i0 <= loop_ub; i0++) {
  360       acnorm_div_data[i0] = tmp_data[i0];
  361     }
  362   
  363     flipud(acnorm_div_data, &acnorm_div_sizes);
  364     b_filter(acnorm_div_data, *(int32_T (*)[1])&acnorm_div_sizes, tmp_data,
  365              &tmp_sizes);
  366   
  367     /* 'filterForDisplay:64' envT = flipud(envT); */
  368     acnorm_div_sizes = tmp_sizes;
  369     loop_ub = tmp_sizes - 1;
  370     for (i0 = 0; i0 <= loop_ub; i0++) {
  371       acnorm_div_data[i0] = tmp_data[i0];
  372     }
  373   
  374     flipud(acnorm_div_data, &acnorm_div_sizes);
  375   
  376     /* 'filterForDisplay:66' acnorm_div = envT.*2; */
  377     acnorm_div_sizes--;
  378     for (i0 = 0; i0 <= acnorm_div_sizes; i0++) {
  379       acnorm_div_data[i0] *= 2.0;
  380     }
  381   
  382     /* 'filterForDisplay:68' value = ac./acnorm_div; */
  383     acnorm_div_sizes = ac_sizes - 40;
  384   
  385     /* 'filterForDisplay:69' value = value(envFilterWidth:end-envFilterWidth); */
  386     if (40 > acnorm_div_sizes) {
  387       i0 = 0;
  388       acnorm_div_sizes = 0;
  389     } else {
  390       i0 = 39;
  391     }
  392   
  393     loop_ub = ac_sizes - 1;
  394     for (tmp_sizes = 0; tmp_sizes <= loop_ub; tmp_sizes++) {
  395       b_ac_data[tmp_sizes] = ac_data[tmp_sizes] / acnorm_div_data[tmp_sizes];
  396     }
  397   
  398     value_sizes[0] = acnorm_div_sizes - i0;
  399     loop_ub = (acnorm_div_sizes - i0) - 1;
  400     for (tmp_sizes = 0; tmp_sizes <= loop_ub; tmp_sizes++) {
  401       value_data[tmp_sizes] = b_ac_data[i0 + tmp_sizes];
  402     }
  403   }
  404   
  405   void hr_analyzer_initialize(void)
  406   {
  407     rt_InitInfAndNaN(8U);
  408   }
  409   
  410   void hr_analyzer_terminate(void)
  411   {
  412     /* (no terminate code required) */
  413   }
  414   
  415   /*
  416    * function [ out ] = hr_dummy( in )
  417    */
  418   real_T hr_dummy(real_T in)
  419   {
  420     /* 'hr_dummy:3' out = in; */
  421     return in;
  422   }
  423   
  424   /* End of code generation (hr_analyzer.c) */
  425