dgt_fb.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 dgt_fb(ltfat_complex *f, ltfat_complex *g,
	    const int L, const int gl,
	    const int W, const int R, const int a, const int M, 
	    ltfat_complex *cout)
{
   /*  --------- initial declarations -------------- */

   int b, N;
   int glh, glh_d_a;

   int k, l, m, n, r, w;

   ltfat_complex *gw, *fw;

   fftw_plan p_veryend;

   ltfat_complex *coefsum;

   div_t domod;
   
   /*  ----------- calculation of parameters and plans -------- */
   
   b=L/M;
   N=L/a;

   /* This is a floor operation. */
   glh=gl/2;

   /* This is a ceil operation. */
   glh_d_a=(int)ceil((glh*1.0)/(a));
   
   gw = ltfat_malloc(gl*R*sizeof(ltfat_complex));
   fw = ltfat_malloc(gl*W*sizeof(ltfat_complex));

   /* Create plan. In-place. */
   p_veryend = fftw_plan_many_dft(1, &M, N*R*W,
				  cout, NULL,
				  1, M,
				  cout, NULL,
				  1, M,
				  FFTW_FORWARD, FFTW_OPTITYPE);
      

   /*  ---------- main body ----------- */
  
   /* Do the fftshift of g to place the center in the middle and
    * conjugate it.
    */
   
   for (r=0;r<R;r++)
   {
      for (l=0;l<glh;l++)
      {
	 gw[l+gl*r]=conj(g[l+(gl-glh)+gl*r]);
      }
      for (l=glh;l<gl;l++)
      {
	 gw[l+gl*r]=conj(g[l-glh+gl*r]);
      }
   }

   /*----- Handle the first boundary using periodic boundary conditions.*/
   for (n=0; n<glh_d_a; n++)
   {
      
      /* Matlab code for this:
       * fpart=[f(L-(glh-n*a)+1:L,:);...
       * f(1:gl-(glh-n*a),:)];
       * 
       * fg=fpart.*gw;
       */
      for (r=0;r<R;r++)
      {
	 for (w=0;w<W;w++)
	 {
	    for (l=0;l<glh-n*a;l++)
	    {
	       fw[l+gl*w]=f[L-(glh-n*a)+l+L*w]*gw[l+r*R];
	    }
	    for (l=glh-n*a;l<gl;l++)
	    {
	       fw[l+gl*w]=f[l-(glh-n*a)+L*w]*gw[l+r*R];
	    }
	 }	 
      }

      /* Do the sum (decimation in frequency, Poisson summation) */      
      /*      
       *coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);
       */
      for (m=0;m<M;m++)
      {
	 /* Place coefficient correctly such that we get a
	  * frequency invariant system. Add 2*L to make sure it is always
	  * positive.
	  */
	 domod = div(2*L+m+(n*a-glh), M);
	 coefsum=cout+n*M+domod.rem;
	 for (w=0;w<W;w++)
	 {	 	  
	    coefsum[w*M*N]=0.0;
	    for (k=0;k<gl/M;k++)
	    {
	       coefsum[w*M*N]+= fw[m+k*M+gl*w];
	    }
	 }
      }	     
   }
   
   /* ----- Handle the middle case. --------------------- */
   for (n=glh_d_a; n<(L-(gl+1)/2)/a+1; n++)
   {

      /*
	fg=f(n*a-glh+1:n*a-glh+gl,:).*gw;
	
	% Do the sum (decimation in frequency, Poisson summation)
	coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);
	end;
      */
      for (r=0;r<R;r++)
      {
	 for (w=0;w<W;w++)
	 {
	    for (l=0;l<gl;l++)
	    {
	       fw[l+gl*w]=f[l+n*a-glh+L*w]*gw[l+r*R];
	    }
	 }	 
      }

      /* Do the sum (decimation in frequency, Poisson summation) */      
      /*      
       *coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);
       */
      for (m=0;m<M;m++)
      {
	 /* Place coefficient correctly such that we get a
	  * frequency invariant system. Add 2*L to make sure it is always
	  * positive.
	  */
	 domod = div(2*L+m+(n*a-glh), M);
	 coefsum=cout+n*M+domod.rem;
	 for (w=0;w<W;w++)
	 {	 	  
	    coefsum[w*M*N]=0.0;
	    for (k=0;k<gl/M;k++)
	    {
	       coefsum[w*M*N]+= fw[m+k*M+gl*w];
	    }
	 }
      }
   }
   
   /* Handle the last boundary using periodic boundary conditions. */   
   for (n=(L-(gl+1)/2)/a+1; n<N; n++)
   {
      /*
	fpart=[f((n*a-glh)+1:L,:);... %   L-n*a+glh) elements
	f(1:n*a-glh+gl-L,:)];  %  gl-L+n*a-glh elements
	fg=fpart.*gw;
	
	% Do the sum (decimation in frequency, Poisson summation)
	coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);      
      */
      for (r=0;r<R;r++)
      {
	 for (w=0;w<W;w++)
	 {
	    for (l=0;l<L-n*a+glh;l++)
	    {
	       fw[l+gl*w]=f[n*a-glh+l+L*w]*gw[l+r*R];
	    }
	    for (l=L-n*a+glh;l<gl;l++)
	    {
	       fw[l+gl*w]=f[l-(L-n*a+glh)+L*w]*gw[l+r*R];
	    }
	 }	 
      }

      /* Do the sum (decimation in frequency, Poisson summation) */      
      /*      
       *coef(:,n+1,:)=sum(reshape(fg,M,gl/M,W),2);
       */
      for (m=0;m<M;m++)
      {
	 /* Place coefficient correctly such that we get a
	  * frequency invariant system. Add 2*L to make sure it is always
	  * positive.
	  */
	 domod = div(2*L+m+(n*a-glh), M);
	 coefsum=cout+n*M+domod.rem;
	 for (w=0;w<W;w++)
	 {	 	  
	    coefsum[w*M*N]=0.0;
	    for (k=0;k<gl/M;k++)
	    {
	       coefsum[w*M*N]+= fw[m+k*M+gl*w];
	    }
	 }
      }

   }

   /* FFT to modulate the coefficients. */
   fftw_execute(p_veryend);   

    /* -----------  Clean up ----------------- */   
   ltfat_free(gw);
   ltfat_free(fw);    
}