wavedlg.cpp

自己设计的新型FFT算法

		
				}
	     
		   }

		   break;
		   default: break;
			   // return CDialog::WindowProc(message,wParam,lParam);
	}
	
	return CDialog::WindowProc(message, wParam, lParam);
}

//////////谱减SS//////////////

//////////谱减SS//////////////

void CWaveDlg::Ondenoise() 
{
	// TODO: Add your control notification handler code here
	int i,j,k,FrameNum=nDataSize/128;//FrameNum=nDataSize/256;
//	long nDataSize_temp=nDataSize;
	short *SpeechTmp1,*start,*temp;
	SpeechTmp1=new short[nDataSize];
	temp=SpeechTmp1;
	start=pSpeech;
	//SpeechTmp2=pSpeech;
    int  length =DEFAULT_LENGTH;
	int  shift=DEFAULT_SHIFT;
//	int shift=DEFAULT_SHIFT;
//	double smul=DEFAULT_MULTIPLE; 
//  double lambda=DEFAULT_SMOOTHING; 
    

    double sum1=0,sum2=0,sum3=0,h=0;
    double win[256], frame[256], ps_noise[256],
		   ps_signal[256],ps_temp[256],frame1[256],abs_signal[256],
	       pmix_signal[256],ps_last[256],frame_last[256];
	double pn[129];//px[129];
	double alfac=1,md=0.73;
    //double pmean[129],fsnr[129],psnr[129];//,gain[129];
	//double pvar[129];
    double pn0[129],pn1[129],pn2[129],pn3[129],pn4[129];
	double frame_ps,frame_pn,ps_final;
    //double pn5[5][129];
	double pn5[645];
	int hh;
	double da=0,xiao=0,equ=0;

struct compx {float real,imag;};   //定义复数结构体
struct compx EE(struct compx,struct compx);
struct compx compx1[257];

struct compx EE(struct compx b1,struct compx b2)
{struct compx b3;
b3.real =b1.real*b2.real-b1.imag*b2.imag;
b3.imag =b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}
//汉明窗
    for(i=0;i<length;i++)
		win[i]=0.54-0.46*cos(2*PI*i/(length-1));
//	double mb[129][4];
//	for(i=0;i<129;i++)
//	{
//		for(j=0;j<4;j++)
//		{
//			mb[i][j]=nw/2;
//		}
//	}
//	double osf[129];
//	for(i=0;i<129;i++)
//	    osf[i]=4/(1+i*fs/(nw*fo));
for(i=0;i<5;i++)
{	hh=129*i;
	for(j=0;j<length;j++)
		frame[j]=(double)(*pSpeech++)/pow(2,15);///nw;

	for(j=1;j<=length;j++)
	{
		compx1[j].real=frame[j-1];
		compx1[j].imag=0;
	}
	
       fft(compx1,length);
	
//	rfft(8,frame);
	
	pn5[hh]=(compx1[1].real*compx1[1].real)/256.0;

	for(j=1;j<=length/2;j++)
	{
	pn5[hh+j]=(compx1[j+1].real*compx1[j+1].real+compx1[j+1].imag*compx1[j+1].imag)/256.0;

}
for(j=0;j<129;j++)
//   for(i=0;i<5;i++)
//   pn5[i][j]=0.01;
//for(j=0;j<129;j++)
//   for(i=0;i<5;i++)
//	   pn[j]+=pn5[i][j];
pn[j]=(pn5[j]+pn5[j+129]+pn5[j+129*2]+pn5[j+129*3]+pn5[j+129*4])/5.0;
//for(i=0;i<129;i++)
//pn[i]=0.05;
//for(i=0;i<129;i++)
//pn[i]=pn[i]/5.0;

ps_noise[0]=pn[0];
ps_noise[length/2]=pn[length/2];
for(i=2;i<129;i++)
ps_noise[i-1]=ps_noise[length+1-i]=pn[i];

for(i=0;i<256;i++)
frame_pn+=ps_noise[i];

	pSpeech=start;

	for(k=0;k<FrameNum-1;k++)
	{
//		for(i=0;i<129;i++)
//			fsnr[i]=0.98*(gain[i]*px[i]/(pn[i]+1.0e-30))+(1-0.98)*max(psnr[i]-1,0);

	for(i=0;i<length;i++)					
		frame[i]=(double)(*pSpeech++)/pow(2,15);///nw;
//		frame[i]=(double)(*pSpeech++)/pow(2,16)/nw;
	
	for(i=0;i<shift;i++)
		start++;   
    pSpeech=start; 

  //    rfft(8,frame);
  	for(j=1;j<=length;j++)
	{
		compx1[j].real=frame[j];
		compx1[j].imag=0;
	}
  	
     fft(compx1,length);
   
    pmix_signal[0]=(pow(compx1[1].real,2)+pow(compx1[1].imag,2))/256.0;
     for(i=1;i<=length/2-1;i++)
     	{
     	pmix_signal[i]=(pow(compx1[i+1].real,2)+pow(compx1[i+1].imag,2))/256.0;
       pmix_signal[length-i]=pmix_signal[i];
     	}
/*mix_signal[0]=(frame[0]*frame[0])/256.0;
 
	
	pmix_signal[length/2]=(frame[1]*frame[1])/256.0;  
 
	for(i=1;i<=length/2-1;i++)
	{
 
		pmix_signal[i]=(pow(frame[2*i],2)+pow(frame[2*i+1],2))/256.0;
 
		pmix_signal[length-i]=pmix_signal[i];
 
	}
	*/
	for(i=0;i<=256;i++)
	{
    ps_temp[i]=pmix_signal[i];
    abs_signal[i]=pow(pmix_signal[i]*256,0.5);
	}

	double c_n=0.995,c_s=0.999,c_beta=0.1;
	for(i=0;i<=256;i++)
	{
		if(k==0)
		{
			ps_noise[i]=c_n*ps_noise[i]+(1-c_n)*pmix_signal[i];
		    ps_signal[i]=pmix_signal[i];//-c_beta*ps_noise[i];(1-c_s)*
			ps_last[i]=ps_signal[i];
		}
        else 
		{
			ps_noise[i]=ps_noise[i];
		    ps_signal[i]=c_s*ps_last[i]+(1-c_s)*pmix_signal[i]-c_beta*ps_noise[i];
		    ps_last[i]=ps_signal[i];
		}
		ps_signal[i]=fabs(ps_signal[i]);
	}
    //维纳滤波
//	double aa=0.7,bb=2,h;
 //   for(i=0;i<256;i++)
 //   frame_ps+=ps_signal[i];
//	ps_final=frame_ps-0.8*frame_pn;
//	h=pow(ps_final/(ps_final+aa*frame_pn),bb);
 //   for(i=0;i<=256;i++)
//	    ps_signal[i]=ps_signal[i]*pow(h,2);


///////////////////////////////////////////////////////////////////
	int ich=18,index1=1,index2=0;//double lf2;
 //   double w=0;
	double b[22],sf[22],c[22],O[22],T[256],f[256],alfa[256],
		beita[256],a1[256],a2[256],a[256],e_ps_signal[256];//0~22
	double cc_temp[22][22];//double t1[22];
	double temp=0,uj=0,ua=0,sfm=0,gm=0,am=0,u=0,mm,temppow,tmax=0,tmin=0,ww0,ww=0;

	for(i=0;i<22;i++)
	{
		b[i]=c[i]=O[i]=sf[i]=0;
	}
	for(i=0;i<256;i++)
	{
		T[i]=0;
	}


    for(i=0;i<=2;i++)
	    b[0]+=ps_signal[i];
	for(i=3;i<=5;i++)
		b[1]+=ps_signal[i];
	for(i=6;i<=9;i++)
		b[2]+=ps_signal[i];
	for(i=10;i<=12;i++)
		b[3]+=ps_signal[i];
	for(i=13;i<=15;i++)
		b[4]+=ps_signal[i];
	for(i=16;i<=19;i++)
		b[5]+=ps_signal[i];
	for(i=20;i<=24;i++)
		b[6]+=ps_signal[i];
	for(i=25;i<=28;i++)
		b[7]+=ps_signal[i];
	for(i=29;i<=34;i++)
		b[8]+=ps_signal[i];
	for(i=35;i<=40;i++)
		b[9]+=ps_signal[i];
	for(i=41;i<=46;i++)
		b[10]+=ps_signal[i];
	for(i=47;i<=54;i++)
		b[11]+=ps_signal[i];
	for(i=55;i<=63;i++)
		b[12]+=ps_signal[i];
	for(i=64;i<=73;i++)
		b[13]+=ps_signal[i];
	for(i=74;i<=85;i++)
		b[14]+=ps_signal[i];
	for(i=86;i<=100;i++)
		b[15]+=ps_signal[i];
	for(i=101;i<=117;i++)
		b[16]+=ps_signal[i];
	for(i=118;i<=140;i++)
		b[17]+=ps_signal[i];
	for(i=141;i<=169;i++)
		b[18]+=ps_signal[i];
	for(i=170;i<=204;i++)
		b[19]+=ps_signal[i];
	for(i=205;i<=245;i++)
		b[20]+=ps_signal[i];
	for(i=246;i<=255;i++)
		b[21]+=ps_signal[i];
		

	
    for(i=0;i<22;i++)
    sf[i]=15.81+7.5*((i+1)+0.474)-17.5*sqrt(1+((i+1)+0.474)*((i+1)+0.474));
    
	for(i=0;i<22;i++)
       for(j=0;j<22;j++)
        cc_temp[i][j]=b[j]*sf[i];
  
    
    for(i=0;i<22;i++)	
		for(j=0;j<22;j++)
			c[i]+=cc_temp[i][j];
        


	double temp_value=0.0;//temp_value in matlab
	for(i=0;i<22;i++)
		temp_value=temp_value+b[i];
	ua=temp_value/256.0;

	temp_value=0;
	for(i=0;i<length;i++)
	{   ww0=fabs(ps_signal[i]);
		ww=log10(ww0);
       temp_value=temp_value+ww;
	}
//	temp_value=temp_value+log10(abs(ps_signal[i]));
	temp_value=temp_value/length;
	uj=pow(10,temp_value);
	sfm=-10*log10(uj/ua);

	u=min(sfm/(-60),1);
	for (i=0;i<=22;i++)
	{
      O[i]=u*(14.5+i)+(1-u)*5.5;
      T[i]=pow(10,log10(fabs(c[i]))-O[i]/10);
	 // T[i]=(10*log10(T[i]*(b[i]/c[i]))+90.302)*fb[i];
	  c[i]=T[i];
	}

	for (i=0;i<=22;i++)
		T[i]=0;


	for(i=0;i<=2;i++)
	    T[i]=c[0];
	for(i=3;i<=5;i++)
		T[i]=c[1];
	for(i=6;i<=9;i++)
		T[i]=c[2];
	for(i=10;i<=12;i++)
		T[i]=c[3];
	for(i=13;i<=15;i++)
		T[i]=c[4];
	for(i=16;i<=19;i++)
		T[i]=c[5];
	for(i=20;i<=24;i++)
		T[i]=c[6];
	for(i=25;i<=28;i++)
		T[i]=c[7];
	for(i=29;i<=34;i++)
		T[i]=c[8];
	for(i=35;i<=40;i++)
		T[i]=c[9];
	for(i=41;i<=46;i++)
		T[i]=c[10];
	for(i=47;i<=54;i++)
		T[i]=c[11];
	for(i=55;i<=63;i++)
		T[i]=c[12];
	for(i=64;i<=73;i++)
		T[i]=c[13];
	for(i=74;i<=85;i++)
		T[i]=c[14];
	for(i=86;i<=100;i++)
		T[i]=c[15];
      for(i=101;i<=117;i++)
		T[i]=c[16];
	for(i=118;i<=140;i++)
		T[i]=c[17];
	for(i=141;i<=169;i++)
		T[i]=c[18];
	for(i=170;i<=204;i++)
		T[i]=c[19];
       for(i=205;i<=245;i++)
		T[i]=c[20];
	for(i=246;i<=255;i++)
		T[i]=c[21];


	

	
	////计算绝对听阈////
     mm=0.0;
	for (i=0;i<length;i++)
	{
		f[i]=mm;
		mm=mm+(double)8/255;
	}
 
	f[0]=f[1];
	for(i=0;i<256;i++)
	{		
		f[i]=(3.64*pow(f[i],-0.8)-6.5*exp(-0.6*pow((f[i]-3.3),2))+0.001*pow(f[i],4));
	}									   
	
	for(i=0;i<256;i++)
	{
		T[i]=maxNum(T[i],f[i]);
//	
//		thv1[i]=pow(10,(T[i]-90.302)/10);
	}
//	for(i=0;i<256;i++)
//	{
//		if(T[i]>tmax)
//			tmax=T[i];
//		if(T[i]<tmin)
//			tmin=T[i];
//    }  
   

//////////////////////旧算法//////////////////////////
/*	for(i=0;i<256;i++)
	{
			alfa[i]=(tmax*6-T[i]*6+1*T[i]-1*tmin)/(tmax-tmin);
			beita[i]=(tmax*0.02-T[i]*0.02+0*T[i]-0*tmin)/(tmax-tmin);      
            temppow=ps_noise[i]/ps_temp[i];
            
			if(pow(temppow,2.0)<(1.0/(alfa[i]+beita[i])))
				//if(1-alfa[k]*pow(temppow,2)>0)
				frame1[i]=pow((1-alfa[i]*pow(temppow,2)),0.5)*ps_temp[i];
            else
				frame1[i]=pow((beita[i]*pow(temppow,2)),0.5)*ps_temp[i];
    }*/


	////////////////////////新算法//////////////////////////////
	for(i=0;i<256;i++)
	{
	a1[i]=pow((ps_signal[i]+ps_noise[i]),2)/(ps_signal[i])-(ps_signal[i]+ps_noise[i]);
    a2[i]=pow((ps_signal[i]+ps_noise[i]),2)/(T[i])-(ps_signal[i]+ps_noise[i]);
    if(a1[i]>a2[i])
         a[i]=a1[i];
    else a[i]=a2[i];
	}

double landa=5;      //  make up for estimate error of a_l & a_h  
double c_a=0.4;  
double above[256],under[256];          // coefficient of a_l & a_h
for(i=0;i<256;i++)
{
above[i]=(double)pow(pmix_signal[i],2);
under[i]=(double)(pmix_signal[i]+a[i]);
e_ps_signal[i]=above[i]/under[i];
//e_ps_signal[i]=(double)pow(pmix_signal[i],2)/(double)(pmix_signal[i]+a[i]);
}
/////////////////////////////////////////////////////////////////
    for(i=0;i<256;i++)
   frame1[i]=pow(e_ps_signal[i],0.5);

	frame[0]=frame[0]*frame1[0]/abs_signal[0];
	frame[length/2]=frame[length/2]*frame1[length/2]/abs_signal[length/2];
//@@	frame[255]=frame[255]*frame1[length/2]/abs_signal[length/2];
 for(i=1;i<length;i+=2)
   //@@	for(i=1;i<length-1;i+=2)
	{
       frame[i]=frame[i]*frame1[(i+1)/2]/abs_signal[(i+1)/2];
	frame[i+1]=frame[i+1]*frame1[(i+1)/2]/abs_signal[(i+1)/2];
//@@ 		frame[i]=frame[i]*frame1[(i+1)/2]/abs_signal[(i+1)/2];
//@@	    frame[i+1]=frame[i+1]*frame1[(i+1)/2]/abs_signal[(i+1)/2];
	}

//反fft	
	// rfft(length,frame,frame,-1);
	    for(i=1;i<=length/2;i++)
	   	compx1[i].real=frame[2*(i-1)];
	       compx1[i].imag=-frame[2*(i-1)+1];
	    for(i=length/2+1;i<=length;i++)
	    	compx1[i].real=frame[2*(i-1)];
	       compx1[i].imag=frame[2*(i-1)+1];
        //irfft(8,frame);
           fft(compx1,length);
        
            for(i=1;i<=length;i++)
	   	compx1[i].real=compx1[i].real/length;
	       compx1[i].imag=-compx1[i].real/length;

	    for(i=1;i<=length;i++)
	   	frame[i-1]=compx1[i].real;
	      	       
 		if(k==0)
		{
			for(i=0;i<length;i++)
			frame_last[i]=frame[i];
		}
        else 
		{
			for(i=0;i<length/2;i++)
				frame[i]=(frame[i]+frame_last[128+i])/2.0;
			for(i=0;i<length;i++)
		    frame_last[i]=frame[i];
		}

	for(i=0;i<length/2;i++)
		//*SpeechTmp1++=(short)(frame[i]);
		*SpeechTmp1++=(short)(frame[i]*pow(2,15)*5);//*nw);

	}
	pSpeech=temp;
}

void CWaveDlg::FFT(struct compx *xin,int N )
{ 
  int f,m,nv2,nm1,i,k,l,j=1;
  int le,lei,ip;
  int b,n;
  float pi,x,y;