wfac.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 "gcd.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;
   div_t domod;

   double sqrtM;

   int ldf;
   
   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);
   
   /* for testing, set ldf=L. */
   ldf=L;
     
   /* Create plan. In-place. */
   p_before = fftw_plan_many_dft(1, &d, c*p*q*R,
				 gf, NULL,
				 c*p*q*R, 1,
				 gf, NULL,
				 c*p*q*R, 1,
				 FFTW_FORWARD, FFTW_OPTITYPE);
   


   ld2=  p*q*R;
   ld3=c*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
		   gf[k+(l+q*w)*p+r*ld2+s*ld3]=sqrtM*g[r+domod.rem+L*w];
#else
		   gf[k+(l+q*w)*p+r*ld2+s*ld3][0]=sqrtM*g[r+domod.rem+L*w][0];
		   gf[k+(l+q*w)*p+r*ld2+s*ld3][1]=sqrtM*g[r+domod.rem+L*w][1];
#endif
		}
	     }
	  }
      }
   }

   fftw_execute(p_before);	  

}


/* wfac for real valued input. */
void wfac_r(double *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;
   div_t domod;

   double sqrtM;

   int ldf;
   
   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);
   
   /* for testing, set ldf=L. */
   ldf=L;
     
   /* Create plan. In-place. */
   p_before = fftw_plan_many_dft(1, &d, c*p*q*R,
				 gf, NULL,
				 c*p*q*R, 1,
				 gf, NULL,
				 c*p*q*R, 1,
				 FFTW_FORWARD, FFTW_OPTITYPE);
   


   ld2=  p*q*R;
   ld3=c*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++)
		{
		  /* Only copy into real part, set imaginary part to zero. */
		   gf[k+(l+q*w)*p+r*ld2+s*ld3][0]=sqrtM*g[r+domod.rem+L*w];
		   gf[k+(l+q*w)*p+r*ld2+s*ld3][1]=0.0;
		}
	     }
	  }
      }
   }

   fftw_execute(p_before);	  

}