123view.cpp

滤波程序的vc++实现

// 123View.cpp : implementation of the CMy123View class
//

#include "stdafx.h"
#include "123.h"

#include "123Doc.h"
#include "123View.h"


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

/////////////////////////////////////////////////////////////////////////////
// CMy123View

IMPLEMENT_DYNCREATE(CMy123View, CView)

BEGIN_MESSAGE_MAP(CMy123View, CView)
	//{{AFX_MSG_MAP(CMy123View)
		// NOTE - the ClassWizard will add and remove mapping macros here.
		//    DO NOT EDIT what you see in these blocks of generated code!
	//}}AFX_MSG_MAP
	// Standard printing commands
	ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)
	ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CMy123View construction/destruction

CMy123View::CMy123View()
{
	// TODO: add construction code here

}

CMy123View::~CMy123View()
{
}

BOOL CMy123View::PreCreateWindow(CREATESTRUCT& cs)
{
	// TODO: Modify the Window class or styles here by modifying
	//  the CREATESTRUCT cs

	return CView::PreCreateWindow(cs);
}

/////////////////////////////////////////////////////////////////////////////
// CMy123View drawing

void CMy123View::OnDraw(CDC* pDC)
{
	CMy123Doc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);

	pDC->SetMapMode(MM_TWIPS);  //一个逻辑单位映射成一个点的1/20大小,x轴正向朝右,y轴正向朝下
	CRect rect;
    GetClientRect(&rect);
	//pDC->LPtoDP(rect);
	int xx=rect.right/2;
	int yy=rect.bottom/2;

	double re[N],img[N];
	for (int i=0; i<N; i++)
       re[i]=sin(PI*100*(1.0/fs)*i);
	pDC->SetViewportOrg(0,yy-200); //设置原点坐标
	pDC->MoveTo(0,0);
	pDC->LineTo(xx*100,0);
	pDC->MoveTo(0,yy);      
	pDC->LineTo(0,-yy);    
	pDC->MoveTo(0,0);
    for (i=0; i<N; i++)
       pDC->LineTo(i*100,re[i]*100);
       pDC->TextOut(13000,0,"sin(PI*100*(1.0/fs)*i)",strlen("sin(PI*100*(1.0/fs)*i)"));

	for (i=0; i<N; i++)
       re[i]=5*sin(PI*100*(1.0/fs)*i);
	pDC->SetViewportOrg(0,yy-150); //设置原点坐标
	pDC->MoveTo(0,0);
	pDC->LineTo(xx*100,0);
	pDC->MoveTo(0,yy);      
	pDC->LineTo(0,-yy);    
	pDC->MoveTo(0,0);
    for (i=0; i<N; i++)
       pDC->LineTo(i*100,re[i]*100);
       pDC->TextOut(13000,0,"5*sin(PI*200*(1.0/fs)*i)",strlen("5*sin(PI*200*(1.0/fs)*i)"));

	for (i=0; i<N; i++)
    {
		re[i]=sin(PI*100*(1.0/fs)*i)+5*sin(PI*200*(1.0/fs)*i);
		img[i]=0;
	}


	pDC->SetViewportOrg(0,yy-100); //设置原点坐标
	pDC->MoveTo(0,0);
	pDC->LineTo(xx*100,0);
	pDC->MoveTo(0,yy);      
	pDC->LineTo(0,-yy);    
	pDC->MoveTo(0,0);
    for (i=0; i<N; i++)
       pDC->LineTo(i*100,re[i]*100);

    //生成滤波器
	double H[N];
	double f;
	for(i=0;i<N/2;i++)
	{
        f=(double)i/N*fs;             // 公式:i/N=f/fs
	    if(f<=50) H[i]=1.0;
		else H[i]=0.0;
	}
	for(i=N/2;i<N;i++)  //滤波器是对称的
	  H[i]=H[N-i-1];


	//正变换
    double fre[N],fimg[N];
	FFT(re,img,N,Bit(N),fre,fimg,0,0);
	pDC->SetViewportOrg(0,yy); //设置原点坐标
	pDC->MoveTo(0,0);
	pDC->LineTo(xx*100,0);
	pDC->MoveTo(0,yy);      
	pDC->LineTo(0,-yy);
	char ch[10];
    for (i=0; i<N; i++)
	{   
		itoa(i,ch,10);
		pDC->TextOut(i*100,-200,ch,strlen(ch));
		pDC->MoveTo(i*100,0);      
	    pDC->LineTo(i*100,-100);
	}

	pDC->MoveTo(0,0);
    for (i=0; i<N; i++)
       pDC->LineTo(i*100,-fre[i]*100);
//       pDC->LineTo(i*100,100*sqrt(fre[i]*fre[i]+fimg[i]*fimg[i]));

	//滤波
   for (i=0; i<N; i++)
	{
		fre[i]=fre[i]*H[i];
		fimg[i]=fimg[i]*H[i];
	}



	//反变换
	FFT(fre,fimg,N,Bit(N),re,img,1,0);
    
	//显示变换后的图形
	pDC->SetViewportOrg(0,yy+200); //设置原点坐标
	pDC->MoveTo(0,0);
	pDC->LineTo(xx*100,0);
	pDC->MoveTo(0,yy);      
	pDC->LineTo(0,-yy);
    
	pDC->MoveTo(0,0);
    for (i=0; i<N; i++)
       pDC->LineTo(i*100,re[i]*100);
       pDC->TextOut(11000,0,"滤波后图像",strlen("滤波后图像"));

	// TODO: add draw code for native data here
}

/////////////////////////////////////////////////////////////////////////////
// CMy123View printing

BOOL CMy123View::OnPreparePrinting(CPrintInfo* pInfo)
{
	// default preparation
	return DoPreparePrinting(pInfo);
}

void CMy123View::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add extra initialization before printing
}

void CMy123View::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
	// TODO: add cleanup after printing
}

/////////////////////////////////////////////////////////////////////////////
// CMy123View diagnostics

#ifdef _DEBUG
void CMy123View::AssertValid() const
{
	CView::AssertValid();
}

void CMy123View::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}

CMy123Doc* CMy123View::GetDocument() // non-debug version is inline
{
	ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CMy123Doc)));
	return (CMy123Doc*)m_pDocument;
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CMy123View message handlers

/*
FFT是变换用的子函数
形参说明：
pr,pi:分别是传 送来 的实部和虚部数组
n：数组长度
k：n是2的k次方
fr，fi：分别是变换后的实部和虚部 
l：为0 做正变换，不为0做反变换
il:il不等于0时上述fr,fi数组分别返回的是变换后的模与幅角。
*/
void CMy123View::FFT(double *pr,double *pi,int n,int k,
					 double *fr,double *fi,int l,int il)
{
   int it,m,is,i,j,nv,l0;
   double p,q,s,vr,vi,poddr,poddi;

   for (it=0; it<=n-1; it++)
      { m=it; is=0;
        for (i=0; i<=k-1; i++)
		{ 
			j=m/2; 
			is=2*is+(m-2*j);
			m=j;
		}
        fr[it]=pr[is]; fi[it]=pi[is];
      }

    pr[0]=1.0; pi[0]=0.0;
    p=6.283185306/(1.0*n);
    pr[1]=cos(p); pi[1]=-sin(p);
    if (l!=0) pi[1]=-pi[1];
    for (i=2; i<=n-1; i++)
      { p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
        s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
        pr[i]=p-q; pi[i]=s-p-q;
      }
    for (it=0; it<=n-2; it=it+2)
      { vr=fr[it]; vi=fi[it];
        fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];
        fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];
      }
    m=n/2; nv=2;
    for (l0=k-2; l0>=0; l0--)
      { m=m/2; nv=2*nv;
        for (it=0; it<=(m-1)*nv; it=it+nv)
          for (j=0; j<=(nv/2)-1; j++)
            { p=pr[m*j]*fr[it+j+nv/2];
              q=pi[m*j]*fi[it+j+nv/2];
              s=pr[m*j]+pi[m*j];
              s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
              poddr=p-q; poddi=s-p-q;
              fr[it+j+nv/2]=fr[it+j]-poddr;
              fi[it+j+nv/2]=fi[it+j]-poddi;
              fr[it+j]=fr[it+j]+poddr;
              fi[it+j]=fi[it+j]+poddi;
            }
      }
    if (l!=0)
      for (i=0; i<=n-1; i++)
        { fr[i]=fr[i]/(1.0*n);
          fi[i]=fi[i]/(1.0*n);
        }
    if (il!=0)
      for (i=0; i<=n-1; i++)
        { pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
          if (fabs(fr[i])<0.000001*fabs(fi[i]))
            { if ((fi[i]*fr[i])>0) pi[i]=90.0;
              else pi[i]=-90.0;
            }
          else
            pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
        }
    return;
}

int CMy123View::Bit(int n)
{
  int k=0;
  while(n)
  {
	  n=n>>1;
	  k++;
  }
  return k-1;
}