iwfac_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 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, ld2n, ld3n;
   
   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_malloc(L*R*sizeof(ltfat_complex));

   /* In-place. Contiguous IFFT. */
   p_before = fftw_plan_many_dft(1, &d, c*p*q*R,
				 tmp_gf, NULL,
				 1, d,
				 tmp_gf, NULL,
				 1, d,
				 FFTW_BACKWARD, 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;
     
   /* Do the dimension permutation. 
    * The r, w and l 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
	 tmp_gf[s+k*d]   =gf[k+s*ld3];
#else
	 tmp_gf[s+k*d][0]=gf[k+s*ld3][0];
	 tmp_gf[s+k*d][1]=gf[k+s*ld3][1];
#endif
      }
   }

   fftw_execute(p_before);	  

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