cantight_fac.c

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

#include "config.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "ltfat_blas.h"
#include "ltfat_lapack.h"
#include "dgt.h"

void cantight_fac(ltfat_complex *gf, const int L, const int R, const int a, const int M, ltfat_complex *gtightf)
{
  
   int N, c, d, p, q, h_a, h_m;
   
   int w, rs;

   ltfat_complex *Sf, *U, *VT, *gfwork;
   double *S;

   ltfat_complex zzero, alpha;
   
   zzero[0]=0.0; zzero[1]=0.0;
   alpha[0]=1.0; alpha[1]=0.0;

   N=L/a;
   
   c=gcd(a, M,&h_a, &h_m);
   p=a/c;
   q=M/c;
   d=N/q;
      
   S  = (double*)ltfat_malloc(p*sizeof(double));
   Sf = (ltfat_complex*)ltfat_malloc(p*p*sizeof(ltfat_complex));
   U  = (ltfat_complex*)ltfat_malloc(p*p*sizeof(ltfat_complex));
   VT = (ltfat_complex*)ltfat_malloc(p*q*R*sizeof(ltfat_complex));
   gfwork = (ltfat_complex*)ltfat_malloc(L*R*sizeof(ltfat_complex));

   /* Copy the contents of gf to gfwork because LAPACK overwrites
    * the input.
    */
   memcpy(gfwork,gf,sizeof(ltfat_complex)*L*R);
  
   for (w=0;w<R;w++)
   {
      for (rs=0;rs<c*d;rs++)
      {
	/* Compute the thin SVD */
	ltfat_zgesvd(p, q*R, gfwork+rs*p*q*R, p,
		     S, U, p, VT, p);

	/* Combine U and V */
	ltfat_zgemm(CblasNoTrans,CblasNoTrans,p,q*R,p,
		    &alpha,U,p,VT,p,
		    &zzero,gtightf+rs*p*q*R, p);


      }
   }

   ltfat_free(gfwork);
   ltfat_free(Sf);
   ltfat_free(S);
   ltfat_free(U);
   ltfat_free(VT);


}