testfilterdlg.cpp

主要功能是任意波段的滤波计算

	}
	else
	{
		AfxMessageBox(_T("N不是2的整数次幂"),MB_OK);

	}
}

//  快速傅立叶变换函数
VOID CTestFilterDlg::kkfft(float pr[],float pi[],int n,int k,float fr[],float fi[],int l,int il)
{
	int it,m,is,i,j,nv,l0;
    float 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.283185306f/(1.0f*n);
    pr[1]=(FLOAT)cos(p); 
	pi[1]=-(FLOAT)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.0f*n);
			fi[i]=fi[i]/(1.0f*n);
        }
	}
	
	if (il!=0)
	{
		for (i=0; i<=n-1; i++)
		{ 
			pr[i]=(FLOAT)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]=(FLOAT)atan(fi[i]/fr[i])*360.0f/6.283185306f;
			}
		}
	}
}

	/// 获得  （去掉直流分量之后） 电压波形 的 均方根
void		CTestFilterDlg::GetWavePower(float	Wave[],  int ptNum, float  * pfPower)
{
	float	fSum=0;
	int		i;
	for(i=0;i<ptNum;i++)
	{
	fSum = fSum + Wave[i]*Wave[i];
	}
	fSum= fSum/ptNum;
	*pfPower= sqrt(fSum);
}


void CTestFilterDlg::OnBtnFft() 
{
	// TODO: Add your control notification handler code here
		// 计算频谱
//	AppFFT( 1024 ,  fAccWave);

	int		i;
	float	xreal[1024],ximag[1024];
	float	fWindow[1024];   //// 加窗口
	float	fFreqDT;
	UpdateData(TRUE);
	fFreqDT =   m_fFsamp/1024;

	for(i=0;i< m_fFreqLow/fFreqDT ;i++)
	{
	fWindow[i]=0.0;
	}
	for(i=m_fFreqLow/fFreqDT;i<m_fFreqHigh/fFreqDT;i++)
	{
	fWindow[i]=1.0;
	}
	for(i=m_fFreqHigh/fFreqDT ;i< 1024-m_fFreqHigh/fFreqDT;i++)
	{
	fWindow[i]=0.0;
	}
	for(i=1024-m_fFreqHigh/fFreqDT;i<1024-m_fFreqLow/fFreqDT;i++)
	{
	fWindow[i]=1.0;
	}
	for(i=1024-m_fFreqLow/fFreqDT;i<1024;i++)
	{
	fWindow[i]=0.0;
	}
	/////////////////

	for(i=0;i<1024;i++)
	{
	xreal[i]	=fAccWave[i];
	ximag[i]	=0.0;
	}
	FFT(xreal,ximag,1024);
	if  (m_bFilter)
	{
		//////裁减频谱，取其中的一部分
		for(i=0;i<1024;i++)
		{
		xreal[i]= xreal[i] * fWindow[i];
		}
	}
	////
	IFFT(xreal,ximag,1024);

	
	for(i=0;i<1024;i++)
	{
	fFreqData[i]=xreal[i] ;
	}
}
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
inline void swap (float &a, float &b)
{
    float t;
    t = a;
    a = b;
    b = t;
}
void bitrp (float xreal [], float ximag [], int n)
{
    // 位反转置换 Bit-reversal Permutation
    int i, j, a, b, p;

    for (i = 1, p = 0; i < n; i *= 2)
		{
        p ++;
        }
    for (i = 0; i < n; i ++)
		{
        a = i;
        b = 0;
        for (j = 0; j < p; j ++)
			{
            b = (b << 1) + (a & 1);    // b = b * 2 + a % 2;
            a >>= 1;			// a = a / 2;
            }
        if ( b > i)
			{
            swap (xreal [i], xreal [b]);
            swap (ximag [i], ximag [b]);
            }
        }
}

void FFT(float xreal [], float ximag [], int n)
{
    // 快速傅立叶变换，将复数 x 变换后仍保存在 x 中，xreal, ximag 分别是 x 的实部和虚部
    float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg;
    int m, k, j, t, index1, index2;
    
    bitrp (xreal, ximag, n);

    // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
    arg = - 2 * PI / n;
    treal = cos (arg);
    timag = sin (arg);
    wreal [0] = 1.0;
    wimag [0] = 0.0;
    for (j = 1; j < n / 2; j ++)
		{
        wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
        wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
        }

    for (m = 2; m <= n; m *= 2)
	{
        for (k = 0; k < n; k += m)
		{
            for (j = 0; j < m / 2; j ++)
				{
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
                timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];
                ureal = xreal [index1];
                uimag = ximag [index1];
                xreal [index1] = ureal + treal;
                ximag [index1] = uimag + timag;
                xreal [index2] = ureal - treal;
                ximag [index2] = uimag - timag;
                }
            }
        }
}


void  IFFT (float xreal [], float ximag [], int n)
{
    // 快速傅立叶逆变换
    float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg;
    int m, k, j, t, index1, index2;
    
    bitrp (xreal, ximag, n);

    // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
    arg = 2 * PI / n;
    treal = cos (arg);
    timag = sin (arg);
    wreal [0] = 1.0;
    wimag [0] = 0.0;
    for (j = 1; j < n / 2; j ++)
		{
        wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
        wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
        }

    for (m = 2; m <= n; m *= 2)
		{
        for (k = 0; k < n; k += m)
			{
            for (j = 0; j < m / 2; j ++)
				{
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
                timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];
                ureal = xreal [index1];
                uimag = ximag [index1];
                xreal [index1] = ureal + treal;
                ximag [index1] = uimag + timag;
                xreal [index2] = ureal - treal;
                ximag [index2] = uimag - timag;
                }
            }
        }

    for (j=0; j < n; j ++)
		{
        xreal [j] /= n;
        ximag [j] /= n;
        }
}



void FFT_test ()
	{
    char	inputfile [] = "input.txt";    // 从文件 input.txt 中读入原始数据
    char	outputfile [] = "output.txt";    // 将结果输出到文件 output.txt 中
    float	xreal [N] ;
	float	ximag [N];
    int n, i;
    FILE *input, *output;

    if (!(input = fopen (inputfile, "r")))
		{
        printf ("Cannot open file. ");
        exit (1);
        }
    if (!(output = fopen (outputfile, "w")))
		{
        printf ("Cannot open file. ");
        exit (1);
        }
        
    i = 0;
    while ((fscanf (input, "%f%f", xreal + i, ximag + i)) != EOF)
		{
        i ++;
        }
    n = i;    // 要求 n 为 2 的整数幂
    while (i > 1)
		{
        if (i % 2)
            {
            fprintf (output, "%d is not a power of 2! ", n);
            exit (1);
            }
        i /= 2;
        }

    FFT (xreal, ximag, n);
    fprintf (output, "FFT:    i     real imag ");
    for (i = 0; i < n; i ++)
		{
        fprintf (output, "%4d    %8.4f    %8.4f ", i, xreal [i], ximag [i]);
        }
    fprintf (output, "================================= ");

    IFFT (xreal, ximag, n);
    fprintf (output, "IFFT:    i     real imag ");
    for (i = 0; i < n; i ++)
		{
        fprintf (output, "%4d    %8.4f    %8.4f ", i, xreal [i], ximag [i]);
        }

    if ( fclose (input))
		{
        printf ("File close error. ");
        exit (1);
        }
    if ( fclose (output))
		{
        printf ("File close error. ");
        exit (1);
        }
}

void CTestFilterDlg::OnBtnReaddata() 
{
	// TODO: Add your control notification handler code here
	CString strTransmit  = _T("");
	DWORD dwCharToWrite  =0;
	DWORD dwBytesWritten =0;
	TCHAR c0, c1 ;
	TCHAR	c2,c3,c4,c5;
	int		i;
	int nValue;
	float	fSum=0.0;
	float	fAvg=0.0;

	GetDlgItemText(IDC_EDIT_RECEIVE_ACC, strTransmit);
	if (strTransmit.GetLength() < (1024*6+8) )
	{
		return ;
	}

	c0=strTransmit.GetAt( 0);
	c1=strTransmit.GetAt( 1);
	//如果为加速度
	if ((c0=='4')  )
	{
	////////////////
	GetDlgItemText(IDC_EDIT_RECEIVE_ACC, strTransmit);
	for (i = 0; i < 1024; i++)
	{
		c0=  strTransmit.GetAt(8+i*6+0);
		c1=  strTransmit.GetAt(8+i*6+1);
		c2=  strTransmit.GetAt(8+i*6+2);
		c3=  strTransmit.GetAt(8+i*6+3);
		c4=  strTransmit.GetAt(8+i*6+4);
		c5=  strTransmit.GetAt(8+i*6+5);

		nValue =  CharToINT(c0,c1,c2,c3,c4,c5);
		fAccWave[i] = 5.0f* (nValue- 512.0f)/512.0f;
		fSum += fAccWave[i];
	}

	fAvg = fSum / 1024.0f;

	// 扣除直流分量，才不会产生静态电压误差
	for (i = 0; i < 1024; ++i)
	{
		fAccWave[i] = fAccWave[i] - fAvg;
	}
	/////////////////////
	}	
}

	// 字符转换为整型
INT CTestFilterDlg::CharToINT(TCHAR c0, TCHAR c1, TCHAR c2,
							TCHAR c3, TCHAR c4, TCHAR c5)
{
	INT   uValue;
	INT	u1,u2;
	unsigned  char	ch1,ch2;

	if  (((c1>='0')&&(c1<='9'))||((c1>='A')&&(c1<='F')) &&
		 ((c2>='0')&&(c2<='9'))||((c2>='A')&&(c2<='F')))
	{
		if ((c1>='0')&&(c1<='9'))
		{
			ch1=c1-48;  // wei shuzi 
		}
		else
		{
			ch1=c1-55;  // wei  zimu
		}

		if ((c2>='0')&&(c2<='9'))
		{
			ch2=c2-48;	// wei shuzi
		}
		else
		{
			ch2=c2-55 ;;   //  zimu
		}

		u1  =  ch1*16+ch2; 
	}

	if  (((c4>='0')&&(c4<='9'))||((c4>='A')&&(c4<='F')) &&
		 ((c5>='0')&&(c5<='9'))||((c5>='A')&&(c5<='F')))
	{
		if  ((c4>='0')&&(c4<='9'))
		{
			ch1=c4-48;  // wei shuzi 
		}
		else
		{
			ch1=c4-55;  // wei  zimu
		}

		if ((c5>='0')&&(c5<='9'))
		{
			ch2=c5-48;	// wei shuzi
		}
		else
		{
			ch2=c5-55 ;;   //  zimu
		}

		u2  =  ch1*16+ch2; 
	}		
	
	uValue = u1*256+ u2;
	
	return  uValue;
}

void CTestFilterDlg::OnBtnJifen() 
{
	// TODO: Add your control notification handler code here
	UpdateData(TRUE);
	int		i;
	fVocWave[0] = fAccWave[0]* m_fDT;
	for(i= 1 ;i<1024;i++)
	{
	fVocWave[i]= fVocWave[i-1] + fAccWave[i]*m_fDT;
	}
	
}