iwfac.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 iwfac(ltfat_complex *gf, const int L, const int R,
	   const int a, const int M, ltfat_complex *g)
{
   
   int b, N, c, d, p, q, h_a, h_m, ld2, ld3;
   
   int l, k, r, s, w;
   div_t domod;
   
   double scaling;
   
   ltfat_complex *tmp_gf;
   fftw_plan p_before;
   
   int ldf;
   
   b=L/M;
   N=L/a;
   
   c=gcd(a, M,&h_a, &h_m);
   p=a/c;
   q=M/c;
   d=b/p;

   /* division by d is because of the way FFTW normalizes the transform. */
   scaling=1.0/sqrt(M)/d;
   
   /* for testing, set ldf=L. */
   ldf=L;

   /* Create plans. In-place.
    * We must copy data to a temporary variable in order not
    * to destroy the input.
    */
   tmp_gf = (ltfat_complex*)ltfat_malloc(L*R*sizeof(ltfat_complex));

   /* This plan is not in-place. It copies to the work array. */
   p_before = fftw_plan_many_dft(1, &d, c*p*q*R,
				 gf, NULL,
				 c*p*q*R, 1,
				 tmp_gf, NULL,
				 c*p*q*R, 1,
				 FFTW_BACKWARD, FFTW_OPTITYPE);
   
  
  
   fftw_execute(p_before);	  

   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	  
		  g[r+domod.rem+L*w]    = tmp_gf[k+(l+q*w)*p+r*ld2+s*ld3]*scaling;
#else
		  g[r+domod.rem+L*w][0] = tmp_gf[k+(l+q*w)*p+r*ld2+s*ld3][0]*scaling;
		  g[r+domod.rem+L*w][1] = tmp_gf[k+(l+q*w)*p+r*ld2+s*ld3][1]*scaling;
#endif
	       }
	    }
	 }
      }
   }           
   
   /* Clear the work-array. */
   ltfat_free(tmp_gf);
}


/* This routine only returns the real values part of the output. Use this ONLY if
 * you know beforehand that the output is real.
 */
void iwfac_r(ltfat_complex *gf, const int L, const int R,
	   const int a, const int M, double *g)
{
   
   int b, N, c, d, p, q, h_a, h_m, ld2, ld3;
   
   int l, k, r, s, w;
   div_t domod;
   
   double scaling;
   
   ltfat_complex *tmp_gf;
   fftw_plan p_before;
   
   int ldf;
   
   b=L/M;
   N=L/a;
   
   c=gcd(a, M,&h_a, &h_m);
   p=a/c;
   q=M/c;
   d=b/p;

   /* division by d is because of the way FFTW normalizes the transform. */
   scaling=1.0/sqrt(M)/d;
   
   /* for testing, set ldf=L. */
   ldf=L;

   /* Create plans. In-place.
    * We must copy data to a temporary variable in order not
    * to destroy the input.
    */
   tmp_gf = (ltfat_complex*)ltfat_malloc(L*R*sizeof(ltfat_complex));

   /* This plan is not in-place. It copies to the work array. */
   p_before = fftw_plan_many_dft(1, &d, c*p*q*R,
				 gf, NULL,
				 c*p*q*R, 1,
				 tmp_gf, NULL,
				 c*p*q*R, 1,
				 FFTW_BACKWARD, FFTW_OPTITYPE);
   
  
  
   fftw_execute(p_before);	  

   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 the real part */
		 g[r+domod.rem+L*w] = tmp_gf[k+(l+q*w)*p+r*ld2+s*ld3][0]*scaling;
	       }
	    }
	 }
      }
   }           
   
   /* Clear the work-array. */
   ltfat_free(tmp_gf);
}