wfac_nos.c

Matlab时频分析工具箱,希望能对大家有所帮助啊

#include "../config.h"
#ifdef HAVE_COMPLEX_H
#include <complex.h>
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <fftw3.h>
#include "dgt.h"

void wfac(ltfat_complex *g, const int L, const int R, const int a, const int M,
	  ltfat_complex *gf)
{
  
   int b, N, c, d, p, q, h_a, h_m;
   
   int w, k, l, r, s, ld2, ld3, ld2n, ld3n;
   div_t domod;

   double sqrtM;

   int ldf;

   ltfat_complex *tmp_gf;
   
   fftw_plan p_before;         

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

   sqrtM=sqrt(M);

   tmp_gf = ltfat_malloc(L*R*sizeof(ltfat_complex));
   
   /* for testing, set ldf=L. */
   ldf=L;
     
   /* Create plan. In-place. Continuous vectors. */
   p_before = fftw_plan_many_dft(1, &d, c*p*q*R,
				 tmp_gf, NULL,
				 1, d,
				 tmp_gf, NULL,
				 1, d,
				 FFTW_FORWARD, FFTW_OPTITYPE);
   


   /* Leading dimensions of the final array. */
   ld2=  p*q*R;
   ld3=c*p*q*R;

   /* Leading dimensions of the temporary array suited for contiguous
    * FFTs.
    */
   ld2n=d*p;
   ld3n=d*p*q*R;
   for (w=0;w<R;w++)
   {
      for (s=0;s<d;s++)
      {
	for (l=0;l<q;l++)
	{
	   for (k=0;k<p;k++)
	   {
	      /*gf(k+1,l+1+q*w,r+1,s+1)=g(r+mod(k*M-l*a+s*p*M,L)+1,w+1);*/
	      /*Add L to make sure it is positive */
	      domod= div(k*M-l*a+s*p*M+L, L);
	      for (r=0;r<c;r++)
	      {
#ifdef HAVE_COMPLEX_H
		 tmp_gf[s+k*d+(l+q*w)*ld2n+r*ld3n]   =sqrtM*g[r+domod.rem+L*w];
#else
		 tmp_gf[s+k*d+(l+q*w)*ld2n+r*ld3n][0]=sqrtM*g[r+domod.rem+L*w][0];
		 tmp_gf[s+k*d+(l+q*w)*ld2n+r*ld3n][1]=sqrtM*g[r+domod.rem+L*w][1];
#endif
	      }
	   }
	}
      }
   }
   
   fftw_execute(p_before);

   /* Do the dimension permutation. 
    * The k, l, r, and w loop have been condensed into the same loop.
    */
   for (s=0;s<d;s++)
   {
      for (k=0;k<p*q*R*c;k++)
      {
#ifdef HAVE_COMPLEX_H
	 gf[k+s*ld3]   =tmp_gf[s+k*d];
#else
	 gf[k+s*ld3][0]=tmp_gf[s+k*d][0];
	 gf[k+s*ld3][1]=tmp_gf[s+k*d][1];
#endif
      }
   }

   /* Clear the work-array. */
   ltfat_free(tmp_gf);


}