dfb1doc.cpp

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

// DFB1Doc.cpp : implementation of the CDFB1Doc class
//

#include "stdafx.h"
#include "DFB1.h"
#include <stdio.h>
#include <math.h>
#include <complex> 
#include <stdlib.h>
using namespace std;

#include "fftw3.h"
#include "DFB1Doc.h"
#include "dialog1.h"
#include "dialog2.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CDFB1Doc

IMPLEMENT_DYNCREATE(CDFB1Doc, CDocument)

BEGIN_MESSAGE_MAP(CDFB1Doc, CDocument)
	//{{AFX_MSG_MAP(CDFB1Doc)
	ON_COMMAND(ID_MENUITEM32771, OnMenuitem32771)
	ON_COMMAND(ID_MENUITEM32772, OnMenuitem32772)
	ON_COMMAND(ID_FILE_SAVE_AS, OnFileSaveAs)
	ON_COMMAND(ID_FILE_OPEN, OnFileOpen)
	ON_UPDATE_COMMAND_UI(ID_MENUITEM32771, OnUpdateMenuitem32771)
	ON_UPDATE_COMMAND_UI(ID_MENUITEM32772, OnUpdateMenuitem32772)
	ON_UPDATE_COMMAND_UI(ID_MENUITEM32773, OnUpdateMenuitem32773)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CDFB1Doc construction/destruction

CDFB1Doc::CDFB1Doc()
{
	// TODO: add one-time construction code here

}

CDFB1Doc::~CDFB1Doc()
{
}

BOOL CDFB1Doc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// TODO: add reinitialization code here
	// (SDI documents will reuse this document)

	return TRUE;
}



/////////////////////////////////////////////////////////////////////////////
// CDFB1Doc serialization

void CDFB1Doc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add loading code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CDFB1Doc diagnostics

#ifdef _DEBUG
void CDFB1Doc::AssertValid() const
{
	CDocument::AssertValid();
}

void CDFB1Doc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CDFB1Doc commands

void CDFB1Doc::DeleteContents() 
{
	// TODO: Add your specialized code here and/or call the base class
	m_ndocchannel=0;    
	m_ndocDM=0;      
	m_ndocstartfreq=0;
	m_ndocendfreq=0;
	m_nreadnum=0;
	m_nnum=0;
	m_nfirnum=0;

	CDocument::DeleteContents();

}

void CDFB1Doc::OnMenuitem32771() 
{
	// TODO: Add your command handler code here
	Cdialog1 dlg;
	dlg.m_nchannel=32;   //32;    
	dlg.m_nDM=0.1;   //0.1520; 
	dlg.m_nstartfreq=6100000; //1638400;    //单位HZ
	dlg.m_nendfreq=6200000;  //1640448;    //	带宽是K/2=4096/2=2048，k＝dlg.m_nendfreq-dlg.m_nstartfreq 
	dlg.m_fir_num=512; 
	dlg.m_bins=256;  
	dlg.m_period=0.901;     //单位ms
	dlg.m_autofirnum =0;
	if(dlg.DoModal()==IDOK)
	{   	
		m_ndocDM=dlg.m_nDM;
        m_ndocstartfreq=dlg.m_nstartfreq;
		m_ndocendfreq=dlg.m_nendfreq;
		double B=m_ndocendfreq-m_ndocstartfreq;
		double fo=(m_ndocendfreq+m_ndocstartfreq)/2;
		//通过DM，B和最低频率可以计算最佳channel数。但是理论值不一定是最佳值，所以也可以自己输入
		if(dlg.m_calcu_channel==TRUE)
		{
			m_ndocchannel=int(pow(2*m_ndocDM/(2.41*m_ndocstartfreq),0.5)*(B/m_ndocstartfreq)*pow(10,8));
		}
		else
		{
			m_ndocchannel=dlg.m_nchannel;
		}
		//如果是从文件中获取fir参数
		if(dlg.m_nfir==TRUE)
		{
			CFileDialog myfiledlg(TRUE,NULL,NULL,0,"文件(*.txt)|*.txt||");
	        if(myfiledlg.DoModal()==IDOK)
			{
				CString strFileName=myfiledlg.GetFileName();
				FILE *stream;    		
				m_nfirnum=10000;      //先设IFR滤波器的参数个数是1024
				m_nfir=new double[m_nfirnum];      //将读到的数据存入m_nfir这个数组中。
				stream = fopen( strFileName, "r+b" );
				m_nfirnum = fread(m_nfir, sizeof( double ),m_nfirnum, stream );    //m_nfirnum必须是fftlength的倍数
				fclose( stream );	
			}
		}			
		//否则自己利用窗函数法生成fir的参数
		else
		{				
			m_nfirnum=dlg.m_fir_num;
			double a=double((m_nfirnum-1))/2;
			m_nfir=new double[m_nfirnum];
			double *hd;
			hd=new double[m_nfirnum];
			double wc=3.1415926/(2*m_ndocchannel);
			for (int r=0;r<m_nfirnum;r++)
			{
				hd[r]=sin(wc*(r-a))/((r-a)*wc);
			}	
			switch(dlg.m_windowtype )
			{
			case 0:		//矩形窗
				for(r=0;r<m_nfirnum;r++)
				{
					m_nfir[r]=hd[r];
				}
				break;
			case 1:     //Bartlett
				for(r=0;r<m_nfirnum/2;r++)
				{
					m_nfir[r]=hd[r]*r*2/m_nfirnum;
				}
				for(r=m_nfirnum/2;r<m_nfirnum;r++)
				{
					m_nfir[r]=hd[r]*(2-r*2/m_nfirnum);
				}
				break;
			case 2:     //Hanning
				for(r=0;r<m_nfirnum;r++)
				{
					m_nfir[r]=hd[r]*0.5*(1-cos(2*3.1415926*r/m_nfirnum));
				}
				break;
			case 3:     //Hamming
				for(r=0;r<m_nfirnum/2;r++)
				{
					m_nfir[r]=hd[r]*(0.54-0.46*cos(2*3.1415926*r/m_nfirnum));
				}
				break;
			case 4:     //balckman
				for(r=0;r<m_nfirnum;r++)
				{
					m_nfir[r]=hd[r]*(0.42-0.5*cos(2*3.1415926*r/m_nfirnum)+0.08*cos(4*3.1415926*r/m_nfirnum));
				}
				break;
			}

		}
		int fftlength=m_ndocchannel*2;				 //例子中fft为64。其实实际中是将数据分为64份，而不是32份。
		int DD=m_nfirnum/fftlength;					 //every subfilter's coefficients  numbers。 每一个subfilter的参数个数。例子中DD＝8；
        int subdatalength=m_nreadnum/fftlength;      //分配给每一个filter处理的数据是输入总的数据量除以fft长度。例子中subdatalength＝1024

/*********这种DFB第一个循环是计算每一个channel中各个数据经过之后的结果，一次有subdatalength＝1024个结果，而第二种DFB的第一个循环是计算每时刻从各个channel出来的数据，每次结果有fftlength＝64个结果，这个结果可以直接fft。
	for (int sub=0;sub<fftlength;sub++)     //计算每一个channel
		{
			m_ncoherent=new double[subdatalength];    
			int firbias=sub*DD;
			for(int q=0;q<(DD-1);q++)
			{
				m_ncoherent[q]=0;             //为了简化卷积运算，所以前面一些不是8个累加之合的结果都设为0；
			}
			for(int i=(DD-1);i<subdatalength;i++)             //下面两个循环完成的是每一个channel中的卷积运算.例子中subdatalength＝1024
			{
				m_ncoherent[i]=0;
                int databias=i*fftlength+sub;
				for(int j=0;j<DD;j++)
				{
					m_ncoherent[i]=m_ncoherent[i]+m_nfir[firbias+j]*m_ndocinputdata[databias-j*fftlength];
				}
			}
*//////////////第二种DFB的第一个循环是计算每时刻从各个channel出来的数据，每次结果有fftlength＝64个结果，这个结果可以直接fft。
     	double *m_nincoherent;    //每次FFT的fftlength个输入值
		m_nincoherent=new double[fftlength]; 
		double ts=B*2*m_ndocDM*pow(10,16)/(2.41*fo*fo*fo*m_ndocchannel);
		double delaynum=ts*2*B/(fftlength);
		//定义一些FFT的相关变量
        int maxdelay=int(delaynum*(m_ndocchannel+1));
        m_ninco_output=new double[subdatalength+maxdelay];    //incoherent's output,这里为1024+37个,因为延时的作用，最高频率的带宽要延后37个单位，所以需要1024＋37个
		//将最终要输出的数组先清零