dfb1doc.cpp

FF快速傅立叶变换算法构建的DFB

		for (int i=0;i<(subdatalength+maxdelay);i++)
		{
			m_ninco_output[i]=0;
		}
	
		for(i=(DD-1);i<subdatalength;i++)      //第一个循环是求出每个时刻从各个channel中输出的值.subdatalength是最大的数值，放在最外一个循环
		{
		//卷积运算
			//考虑DFB中64个subfilter的第一个filter的卷积，因为这个subfilter的参数和其它的subfilter不同
			m_nincoherent[0]=0;
			for(int j=0;j<DD;j++)
			{
				m_nincoherent[0]=m_nincoherent[0]+m_nfir[fftlength*j]*m_ndocinputdata[(i*fftlength-j*fftlength)];
			}
			//再计算余下的63个subfilter的卷积
			for (int sub=1;sub<fftlength;sub++)     //计算每一个channel
			{
                int databias=i*fftlength+sub;     //为了减少计算量。要不下面那个循环的数组index就会写得更复杂
				m_nincoherent[sub]=0;
				for(int jj=0;jj<DD;jj++)
				{
					m_nincoherent[sub]=m_nincoherent[sub]+m_nfir[fftlength*(jj+1)-sub]*m_ndocinputdata[databias-jj*fftlength];
			    	 //原来考虑到subfilters的参数已经按照顺序排好了，所以用的是下面这个循环
					 //m_nincoherent[sub]=m_nincoherent[sub]+m_nfir[fftlength*j+sub]*m_ndocinputdata[databias-j*fftlength];
					 //这种情况下，就不需要单独考虑64个subfilter中那个第一个sudfilter
				}
			}

			//然后将m_nincoherent进行FFT。
			fftw_complex *m_nsuboutput;                                            
			m_nsuboutput= (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (m_ndocchannel+1));   
			fftw_plan p;
			p = fftw_plan_dft_r2c_1d(fftlength, m_nincoherent, m_nsuboutput, FFTW_ESTIMATE);   
   			fftw_execute(p); 
			fftw_destroy_plan(p);
			//做完一个FFT就将channel×2个FFT的输出，做位移，累加到最终的incoherent输出数组m_ninco_output中。这个循环还在主循环内
	    	complex<double> *suboutput;       //定义一个complex类指针，接着把m_nsuboutput的指针赋值给这个指针。这样就可以用complex类的重载运算符计算由fftw计算得来的复数。整个过程中不用对他分配空间。
                                              //他和上面那个指针其实指的是一个地方。之所以又定义了一个指针，只是为了这种指针类型所指的量可以进行复数运算。这种方法就不需要将上面那个指针所指的数据空间重新排列或者重新定义一个空间存放，大大减少了计算量和存储空间
			suboutput=(complex<double> *)m_nsuboutput;    //complex类指针disp_fft_complex指向fftw_complex类型的二维数组。这样才可以使用complex.h函数中的运算符重载，也就是说才可以进行复数运算
			for(sub=0;sub<(m_ndocchannel+1);sub++)
			{
			  int intdelay=int(delaynum*sub);
			  m_ninco_output[i+intdelay]=m_ninco_output[i+intdelay]+abs(suboutput[sub]);    //shifting
			}
			fftw_free(m_nsuboutput);   //释放FFT输出值的空间
            numspp=B/m_ndocchannel*dlg.m_period;   //用于unfolded画图，这个代表显示一个周期需要的点数
		}

		delete []m_nincoherent;    //释放FFT输入值的空间
		delete []m_nfir; 
	    //开始通过积分器
        double temp=double(((dlg.m_period)/dlg.m_bins)/(m_ndocchannel/B));
		int intnum=((dlg.m_period)/dlg.m_bins)/(m_ndocchannel/B);
		if ((temp-intnum)>0.5)
			intnum=intnum+1;      //四舍五入
		if (dlg.m_autobins==TRUE)
			bins=dlg.m_period/(m_ndocchannel/B)/intnum;
		else
			bins=dlg.m_bins;
        double period1=bins*m_ndocchannel/B*intnum;     //实际那么多bins数能表达的周期时间
		double temp1=(dlg.m_period-period1)/dlg.m_period;     //每次连续积分continousN，然后重新设定时间。
		if (temp1<0)
			temp1=-1*temp1;
		double coeff=0.01;         //这个值代表了可以忍受的由于bins的分辨率而带来的time smear
        int continousN=coeff/temp1;
		if (continousN==0)
		{
			AfxMessageBox("too less bins! reset the bins number!",MB_ICONWARNING);
			bins=dlg.m_period/(m_ndocchannel/B)/intnum;
			period1=bins*m_ndocchannel/B*intnum;
			temp1=(dlg.m_period-period1)/dlg.m_period; 
			if (temp1<0)
				temp1=-1*temp1;
			continousN=0.01/temp1;
		}
		m_folded_output=new double[bins];
        for (i=0;i<bins;i++)
			m_folded_output[i]=0;       //清零
		int togglenum=subdatalength/(continousN*intnum*bins);       //togglenum是重新定标的次数
		int step=0;
		if (togglenum==0)
		{
			togglenum=1;
			continousN=subdatalength/(intnum*bins);
		}
		for (int p=0;p<togglenum;p++)              //重新定标
			for (int k=0;k<continousN;k++)         //连续积分的次数。积分完那么多个周期之后，需要重新定标
				for (int i=0;i<bins;i++)    //积分到bins数以后，代表积分了一个周期
					for (int j=0;j<intnum;j++)    //每intnum个点积分一次
					{
						step=i*intnum+j+k*intnum*bins+int(p*continousN*(dlg.m_period*B/m_ndocchannel));
						m_folded_output[i]=m_folded_output[i]+m_ninco_output[step];      //dlg.m_period*B/m_ndocchannel代表每个实际周期的点数。
					}
						int aaa=1;
	}
		
}

void CDFB1Doc::OnMenuitem32772()    //coherent消色散处理函数
{
	// TODO: Add your command handler code here
	Cdialog2 dlg;
	dlg.m_nDM=0.1;      //DM
	dlg.m_nstartfreq=6100000;
	dlg.m_nendfreq=6200000;    //	带宽是K/2=4096/2=2048，k＝dlg.m_nendfreq-dlg.m_nstartfreq
	if(dlg.DoModal()==IDOK)
	{
		//dialog向doc类传值
		m_ndocDM=dlg.m_nDM;
        m_ndocstartfreq=dlg.m_nstartfreq;
		m_ndocendfreq=dlg.m_nendfreq;	

		//定义一些指针和数据
		fftw_complex *disp_fft;    //为了节约内存空间，只定义了1个复数fftw的指针，他开始致输入实数的fft变换之后一半数据（实数fft变换是对称的），之后存放乘上消色散传输函数之后的频谱。
		complex<double> *disp_fft_copy;       //定义一个complex类指针，接着把上面的指针赋值给这个指针。这样就可以用complex类的重载运算符计算由fftw计算得来的复数。整个过程中不用对他分配空间。
                                              //他和上面那个指针其实指的是一个地方。之所以又定义了一个指针，只是为了这种指针类型所指的量可以进行复数运算。这种方法就不需要将上面那个指针所指的数据空间重新排列或者重新定义一个空间存放，大大减少了计算量和存储空间
		double f0=(m_ndocstartfreq+m_ndocendfreq)/2;   //中心频率
		complex<double> constant;     //一个复数常量，用以消色散计算的简化
		constant.real(0);
		constant.imag(-2*3.1415926*m_ndocDM/(2.41e-16)/(f0*f0*f0));	//因为消色散的传输函数是：exp{-i*2pi*(f-f0)^2*D/f0^3},所以定义constant为-i*2pi*D/f0^3
		complex<double> freq;     //常量声明。这个存放频率
		freq.real(m_ndocstartfreq-f0);          //因为消色散的传输函数是：exp{-i*2pi*(f-f0)^2*D/f0^3}，所以计算freq的时候，实际上算的是f-f0
		freq.imag(0);

		//开始计算fft
		fftw_plan p;
		disp_fft = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (m_nreadnum/2+1));    //一组实数的fft计算的结果由这个指针指出。并事先分配好空间。
		p = fftw_plan_dft_r2c_1d(m_nreadnum, m_ndocinputdata, disp_fft, FFTW_ESTIMATE);
		fftw_execute(p);     
		fftw_destroy_plan(p);
		//FFT计算结束

		//频域coherent消色散计算开始
		disp_fft_copy=(complex<double> *)disp_fft;    //complex类指针disp_fft_complex指向fftw_complex类型的二维数组。这样才可以使用complex.h函数中的运算符重载，也就是说才可以进行复数运算
        disp_fft_copy[0]=disp_fft_copy[0]*exp(constant*freq*freq);
		complex<double> step;     //因为要求复数运算中所有的项目都是复数，所以只有把step设为复数
		step.real((m_ndocendfreq-m_ndocstartfreq)*2/m_nreadnum);
		step.imag(0);
		for (int j=1;j<(m_nreadnum/2+1);j++)
		{
			freq=freq+step;     //得到freq[0].....freq[m_nreadnum/2],总共m_nreadnum/2+1个数据。而实数fftw的结果也是N/2+1个有效数据。符合。
			disp_fft_copy[j]=disp_fft_copy[j]*exp(constant*freq*freq);
		}
        
		
        //做逆FFT变换
		fftw_plan pp;
		pp = fftw_plan_dft_c2r_1d(m_nreadnum, disp_fft, m_ndocinputdata, FFTW_ESTIMATE);//默认的c2r的fftw将把输入的复数destory。可以通过将FFTW_ESTIMATE替换成FFTW_PRESERVE_INPUT参数来预防这种destory，但是fftw的效率将会降低。
		fftw_execute(pp);          //fftw中调用IFFT其实就是调用fft，只不过e的符号是负号还是正好。负号代表正变换，正好代表逆变换。所以调用IFFT的时候，需要将结果除于N＝m_nreadnum
		fftw_destroy_plan(pp);
//      fftw_free(disp_fft);    //如果上面line175中的参数改成保留输入值，则必须要添加这个代码释放空间

			}
}


/*  --------------------------------------------------------------------------//复数fft	
		fftw_complex *in;
	    in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * m_nreadnum);	
		for (int i=0;i<m_nreadnum;i++)
		{
			in[i][0]=m_ndocinputdata[i];       //real data FFT
			in[i][1]=0;
		}
		p = fftw_plan_dft_1d(m_nreadnum, in, disp_fft, FFTW_FORWARD, FFTW_ESTIMATE);
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－*/

/*－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
		for (int j=0;j<m_nreadnum;j++)
		{
			m_ndocinputdata[j]=sqrt(disp_fft[j][0]*disp_fft[j][0]+disp_fft[j][1]*disp_fft[j][1]);   用于计算和显示fft之后的输入数据的频谱幅度值
		}
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－*/


void CDFB1Doc::OnFileSaveAs() 
{
	// TODO: Add your command handler code here
	CFileDialog myfiledlg(FALSE,NULL,NULL,0,"文件(*.txt)|*.txt||");
	if(myfiledlg.DoModal()==IDOK)
	{
		CString strFileName=myfiledlg.GetFileName();
		CFile myfile(strFileName, CFile::modeCreate|CFile::modeWrite);
//	    myfile.Write(&m_nc,30);
	}
	
}

void CDFB1Doc::OnFileOpen()    //打开文件
{
	// TODO: Add your command handler code here
	CFileDialog myfiledlg(TRUE,NULL,NULL,0,"文件(*.txt)|*.txt||");
	if(myfiledlg.DoModal()==IDOK)
	{
		CString strFileName=myfiledlg.GetFileName();
		FILE *stream;    		
		m_nnum=10000000;   //一次读的数据大小,最终必须通过别的途径输入，现在只是先假设一下			
		m_ndocinputdata=new double[m_nnum];      //将读到的数据存入m_ndocinputdata[m_ndocN]这个数组中。
		stream = fopen( strFileName, "r+b" );
		m_nreadnum = fread( m_ndocinputdata, sizeof( double ), m_nnum, stream );
		if (m_nreadnum!=m_nnum)
		{
			AfxMessageBox("The data capacity actually read out is less than what you desire!",MB_ICONWARNING);
		}
		fclose( stream );	

		//open the binary file using the fread and fopen function
		//the next step is to display the data using GDI
		
	}


}




void CDFB1Doc::OnUpdateMenuitem32771(CCmdUI* pCmdUI) 
{
	// TODO: Add your command update UI handler code here
	if(m_nreadnum==0)  pCmdUI->Enable(FALSE);    //如果没有输入数据，即实际读入的数据m_nreadnum==0,则将该按钮变灰
	else pCmdUI->Enable(TRUE);
}

void CDFB1Doc::OnUpdateMenuitem32772(CCmdUI* pCmdUI) 
{
	// TODO: Add your command update UI handler code here
	if(m_nreadnum==0)  pCmdUI->Enable(FALSE);    //如果没有输入数据，即实际读入的数据m_nreadnum==0,则将该按钮变灰
	else pCmdUI->Enable(TRUE);
}


void CDFB1Doc::OnUpdateMenuitem32773(CCmdUI* pCmdUI) 
{
	// TODO: Add your command update UI handler code here
	if(m_nreadnum==0)  pCmdUI->Enable(FALSE);    //如果没有输入数据，即实际读入的数据m_nreadnum==0,则将该按钮变灰
	else pCmdUI->Enable(TRUE);
}