⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 新建 文本文档.txt

📁 自己设计的新型FFT算法
💻 TXT
字号:
🔍 
//////////谱减SS//////////////

void CWaveDlg::Ondenoise() 
{
	// TODO: Add your control notification handler code here
	int i,j,FrameNum=nDataSize/256;
	long nDataSize_temp=nDataSize;
	short *SpeechTmp1,*SpeechTmp2,*start,*temp;
	SpeechTmp1=new short[nDataSize];
	temp=SpeechTmp1;
	start=pSpeech;
	SpeechTmp2=pSpeech;
    int  length =DEFAULT_LENGTH;
	int shift=DEFAULT_SHIFT;
	double smul=DEFAULT_MULTIPLE; 
    double lambda=DEFAULT_SMOOTHING; 
    
    
    double *temp1,sum1=0,sum2=0,sum3=0,h=0;
    double win[256], frame[256], ps_noise[256],frame_temp[256], ps_signal[256];
	double pxn[129],os[129],pn[129],px[129];
	double alfac=1,md=0.73;
	double pmean[129],fsnr[129],psnr[129],gain[129];
	double pvar[129];
	//设置
	double fftlen=0.032;int overlap=4;double tm=1.5;int nob=4;
	int ko=4;int fo=800;
	int ni=64;int fs=8000;
	double ti=(double)ni/fs;double zg=exp(-ti/0.04);
	double nw=ni*overlap;
	int nm=1+(int)(fs*tm/(ni*nob));
    for(i=0;i<length;i++)
		win[i]=0.54-0.46*cos(2*PI*i/length);
	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<length;i++)
		frame[i]=(double)(*pSpeech++)/pow(2,16)/nw;//
	rfft(length,frame,frame_temp,1);
	pxn[0]=pn[0]=frame_temp[0]*frame_temp[0];
	pxn[length/2]=ps_noise[length/2]=frame_temp[1]*frame_temp[1];
	for(i=1;i<=length/2-1;i++)
	{
		pxn[i]=pn[i]=frame_temp[2*i]*frame_temp[2*i]+frame_temp[2*i+1]*frame_temp[2*i+1];
		
	}

	for(i=0;i<129;i++)
	{		
		pmean[i]=pn[i];
		pvar[i]=pn[i]*pn[i];
	}

	int im=0;
	for(i=0;i<129;i++)
	{
		fsnr[i]=0;px[i]=0;
		psnr[i]=1;
		gain[i]=0.5;
	}
	
	
	//
	//for(i=0;i<shift;i++)
	//	start++;
	pSpeech=start;

	do
	{
		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++)/nw;
		frame[i]=(double)(*pSpeech++)/pow(2,16)/nw;//
	for(i=0;i<shift;i++)
		start++;   
    pSpeech=start; 
	temp1=frame_temp;
	rfft(length,frame,frame_temp,1);

	//
   
	ps_signal[0]=frame_temp[0]*frame_temp[0];
	
	ps_signal[length/2]=frame_temp[1]*frame_temp[1];                
	for(i=1;i<=length/2-1;i++)
	{
		ps_signal[i]=frame_temp[2*i]*frame_temp[2*i]+frame_temp[2*i+1]*frame_temp[2*i+1];
		ps_signal[length-i]=ps_signal[i];
	}
	
	

///////////////////////////////////////////////////////////////

	int i,j,length=DEFAULT_LENGTH;double a=0,w=0;
	double b[22],sf[22][22],c[22],O[22],T[256],thv1[256];//0~22
	double temp=0,uj=0,ua=0,sfm=0,gm=0,am=0,u=0,f[256],mm;
	for (i=0;i<=21;i++)
	{
		b[i]=c[i]=O[i]=0.0;
		for (j=0;j<=21;j++)
			sf[i][j]=0.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<=21;i++)
		for(j=0;j<=21;j++)
		{
			sf[i][j]=15.81+7.5*((i-j)+0.474)-17.5*sqrt(1+((i-j)+0.474)*((i-j)+0.474));//db
			sf[i][j]=pow(10,sf[i][j]/(20));
		}
	
	for(j=0;j<=21;j++)
	{
		c[j]=0;
		for(i=0;i<=21;i++)
			c[j]=c[j]+b[i]*sf[i][j];
	}
	temp=0.0;
	for(i=0;i<=21;i++)
		temp=temp+b[i];
	ua=temp/256.0;

	temp=0;
	for(i=0;i<length;i++)   
		temp=temp+log10(ps_signal[i]);
	temp=temp/length;
	uj=pow(10,temp);
	sfm=-10*log10(uj/ua);

	u=min(sfm/(-60),1);
	for (i=0;i<=21;i++)
	{
      O[i]=u*(14.5+i)+(1-u)*5.5;
      T[i]=pow(10,log10(c[i])-O[i]/10);
	  T[i]=10*log10(T[i])+90.302;
	  c[i]=T[i];//c[i]暂存T[i]
	}
	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[6]=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);
	}
	double sum1=0,sum2=0;
	for(i=0;i<129;i++)
	{
		sum1=sum1+pxn[i];
		sum2=sum2+ps_signal[i];   
	}
	double alfac1;double gama[129],za[129];
	alfac1=1/pow((1+sum1/sum2-1),2);
	alfac=0.7*alfac+0.3*max(alfac1,0.7);
	for(i=0;i<129;i++)
	{
		gama[i]=(pxn[i]+1.0e-30)/(pn[i]+1.0e-30);
		za[i]=max(0.96*alfac/(1+(gama[i]-1)*(gama[i]-1)),0.3);
	}
	for(i=0;i<129;i++)
	{
		pxn[i]=za[i]*pxn[i]+(1-za[i]*ps_signal[i]);
	}
	im=fmod(im+1,nm);
	if (im)
	{
		for(i=0;i<129;i++)
			mb[i][0]=min(mb[i][0],pxn[i]);
	}
	else
	{
		for(j=3;j>=1;j--)
		{
			for(i=0;i<129;i++)
				mb[i][j]=mb[i][j-1];
		}
		for(i=0;i<129;i++)
			mb[i][0]=pxn[i];
	}
	double beta[129],fvar[129],qeqinv[129],qeq[129],bc[129];
	double qinv=0;
	for(i=0;i<129;i++)
	{
		beta[i]=min(za[i]*za[i],0.8);
		pmean[i]=beta[i]*pmean[i]+(1-beta[i])*pxn[i];
		pvar[i]=beta[i]*pvar[i]+(1-beta[i])*pxn[i]*pxn[i];
		fvar[i]=pvar[i]-pmean[i]*pmean[i];
		qeqinv[i]=min(fvar[i]/(2*pn[i]*pn[i]),0.5);
		qeq[i]=(1/qeqinv[i]-2*md)/(1-md);
		bc[i]=1+(nm-1)*2/qeq[i];
		qinv=qinv+qeqinv[i];
	}
	qinv=qinv/129;
	double binmax,bcc;
	if (qinv<0.03)
		binmax=8;
	else if (qinv<0.05)
		binmax=4;
	else if (qinv<0.06)
		binmax=2;
	else 
		binmax=1.2;
	bcc=1+2.12*sqrt(qinv);
	double bin[129];
	for(i=0;i<129;i++)
	{
		bc[i]=min(bc[i],binmax);
		bin[i]=bcc*bc[i];
	}
	temp=1e10;
	for(i=0;i<129;i++)
	{
		for(j=0;j<4;j++)
		{
			if(mb[i][j]<temp)
				temp=mb[i][j];
		}
		pn[i]=bin[i]*temp;
	}
	double galfa[129],gbeta[129],gain[129];
	for(i=0;i<129;i++)
	{
		os[i]=za[i]*os[i]+(1-za[i])*(1+osf[i]*pn[i]/(pn[i]+pxn[i]));
		px[i]=zg*px[i]+(1-zg)*ps_signal[i];
		psnr[i]=ps_signal[i]/(pn[i]+1.0e-30);
		galfa[i]=15*thv1[i]/(fabs(fsnr[i]-psnr[i]+1)*pn[i]);
		gbeta[i]=1+((galfa[i]-1)/galfa[i])*fsnr[i];
		gain[i]=sqrt(max((galfa[i]-galfa[i]*gbeta[i]/psnr[i]),0));
	}
	frame_temp[0]=frame_temp[0]*gain[0];
	frame_temp[length/2]=frame_temp[length/2]*gain[length/2];
	for(i=1;i<=length/2-1;i++)
	{
		frame_temp[2*i]=frame_temp[2*i]*gain[i]*nw;
		frame_temp[2*i+1]=frame_temp[2*i+1]*gain[i]*nw;
	}
	
	rfft(length,frame_temp,frame,-1);
	
	for(i=length/2;i<length;i++)
		//*SpeechTmp1++=(short)(frame[i]);
		*SpeechTmp1++=(short)(frame[i]*pow(2,16));//

	}while(FrameNum--);
	pSpeech=temp;
}

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -