comp_iwfac.c

来自「Matlab时频分析工具箱,希望能对大家有所帮助啊」· C语言 代码 · 共 78 行

#include "mex.h"
#include "config.h"
#include "../src/iwfac.c"
#include "../src/gcd.c"

/* Calling convention:
 *  comp_iwfac(gf,L,a,M);
 */

void mexFunction( int nlhs, mxArray *plhs[], 
		  int nrhs, const mxArray *prhs[] )
     
{ 

  int pqR, L, R, a, M, N, c, d, p, q,h_a,h_m;
  double *g_combined, *gf_combined;
  double *g_r,*g_i,*gf_r,*gf_i;

  int ii;

  /* Get matrix dimensions.*/
  pqR = mxGetM(prhs[0]); 
  
  L=(int)mxGetScalar(prhs[1]);
  a=(int)mxGetScalar(prhs[2]);
  M=(int)mxGetScalar(prhs[3]);
  
  N=L/a;

  c=gcd(a, M,&h_a, &h_m);
  p=a/c;
  q=M/c;
  d=N/q;

  R=pqR/(p*q);
    
  /* Create temporary matrices to convert to correct complex layout. */

  g_combined=mxCalloc(L*R*2,sizeof(double));
  gf_combined=mxCalloc(L*R*2,sizeof(double));

  /* Copy the data. */

  g_r=mxGetPr(prhs[0]);
  for (ii=0;ii<L*R; ii++)
    {
      g_combined[2*ii]=g_r[ii];
    }

  if (mxIsComplex(prhs[0]))
    {
      g_i=mxGetPi(prhs[0]);
      for (ii=0;ii<L*R; ii++)
	{
	  g_combined[2*ii+1]=g_i[ii];
	}
    }

  iwfac((ltfat_complex*)g_combined, L, R, a, M, (ltfat_complex*)gf_combined);

  plhs[0] = mxCreateDoubleMatrix(L, R, mxCOMPLEX);

  gf_r=mxGetPr(plhs[0]);
  gf_i=mxGetPi(plhs[0]);
  
  for (ii=0;ii<L*R; ii++)
    {
      gf_r[ii]=gf_combined[2*ii];
      gf_i[ii]=gf_combined[2*ii+1];
    }

  mxFree(gf_combined);
  return;
  
}