firdesign.cpp

可以实现FIR滤波器任意阶数和窗函数下的设计。并且有可以在VC下显示幅频响应曲线和滤波效果。

// FIRDesign.cpp: implementation of the CFIRDesign class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Oscillograph.h"
#include "FIRDesign.h"
#include <stdlib.h>
#include <math.h>

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CFIRDesign::CFIRDesign()
{
	m_pWindow = NULL;		//窗函数数组指针
	m_pFIR = NULL;			//FIR滤波器冲击响应函数数组指针

	m_nWindowLength = 0;	//窗函数长度（与滤波器配合使用时，窗函数长度＝滤波器长度＝滤波器阶数＋1）
	m_nFilterOrder = 0;		//滤波器阶数（滤波器长度＝滤波器阶数＋1）
}

CFIRDesign::~CFIRDesign()
{
	if(m_pWindow)
	{
		delete [] m_pWindow;
		m_pWindow = NULL;		//窗函数数组指针
	}
	if(m_pFIR)
	{
		delete [] m_pFIR;
		m_pFIR = NULL;			//FIR滤波器冲击响应函数数组指针
	}
}

//窗函数设计
/*
#define	BOXCAR		1		//矩形窗
#define	BARTLETT	2		//巴特利特窗
#define	HANNING		3		//汉宁窗
#define	HAMMING		4		//哈明窗
#define	BLACKMAN	5		//布莱克曼窗*/
void CFIRDesign::Window(int M, int Type)
{
	if(m_pWindow)
		delete [] m_pWindow;

	m_nWindowLength = M;	//窗函数长度
	m_pWindow = new double[M];
	int i;
	switch(Type)
	{
		//矩形窗
	case BOXCAR:
		for(i=0; i<M; i++)
			*(m_pWindow+i) = 1;
		break;

		//巴特利特窗（窗函数长度为奇数）
	case BARTLETT:

		if(M%2 == 1)	//M为奇数
		{
			for(i=0; i<(M-1)/2; i++)
				*(m_pWindow+i) = 2.0*i/(M-1);

			*(m_pWindow+(M-1)/2) = 1;

			for(i=(M+1)/2; i<M; i++)
				*(m_pWindow+i) = 2 - 2.0*i/(M-1);
		}
		else	//M为偶数
		{
			for(i=0; i<M/2; i++)
				*(m_pWindow+i) = 2.0*i/(M-1);

			for(i=M/2; i<M; i++)
				*(m_pWindow+i) = 2 - 2.0*i/(M-1);
		}
		break;

		//汉宁窗
	case HANNING:
		for(i=0; i<M; i++)
			*(m_pWindow+i) = 0.5 - 0.5*cos(2*PI*(i+1)/(M+1));
//		double temp[10];
//		::memcpy(temp,m_pWindow,8*10);
		break;

		//哈明窗
	case HAMMING:
		for(i=0; i<M; i++)
			*(m_pWindow+i) = 0.54 - 0.46*cos(2*PI*i/(M-1));
		break;

		//布莱克曼窗
	case BLACKMAN:
		for(i=0; i<M; i++)
			*(m_pWindow+i) = 0.42 - 0.5*cos(2*PI*i/(M-1)) + 0.08*cos(4*PI*i/(M-1));
		break;
	default:
		break;
	}
}

//低通、高通滤波器设计
/*
#define	LP			1		//LowPass,低通滤波器
#define	HP			2		//HighPass,高通滤波器
#define	BP			3		//BandPass,带通滤波器
#define	BS			4		//BandStop,带阻滤波器*/
//滤波器阶数（filter_order) = filter_length - 1;
//根据传递函数得到：b(z)=b(1)+b(2)*z(-1)+b(3)*z(-2)+……+b(N+1)*z(-N)
//N：滤波器阶数
//Wc：归一化截止频率（0－1）对应（0－PI），已经归一化
//Type：滤波器类型：LP、HP
void CFIRDesign::FIR(int N, double Wc, int Type, int Window)
{
	//得到滤波器阶数N并创建数据存储数组
	ASSERT(N % 2 == 0); //保证N值为偶数（目前的滤波器设计算法只能在阶数为偶数情况下有效）
	if(m_pFIR)
		delete [] m_pFIR;

	this->m_nFilterOrder = N;
	this->m_pFIR = new double [N+1];
	//窗函数设计
	this->Window(N+1,Window);
	//得到滤波器冲击响应函数
	int i;
	switch(Type)
	{
	//LowPass,低通滤波器
	case LP:
		for(i=0; i<=N; i++)
		{
			m_pFIR[i] = sin((i-N/2)*Wc*PI) / (PI*(i-N/2)) * m_pWindow[i];
			if(i == N/2)
				m_pFIR[i] = Wc;
		}
//		double temp[41];
//		::memcpy(temp,m_pFIR,8*41);
		break;
	//HighPass,高通滤波器
	case HP:
		for(i=0; i<=N; i++)
		{
			m_pFIR[i] = - sin((i-N/2)*Wc*PI) / (PI*(i-N/2)) * m_pWindow[i];
			if(i == N/2)
				m_pFIR[i] = 1 - Wc;
		}
		break;
	default:
		break;
	}
}

//带通、带阻滤波器设计
//滤波器阶数（filter_order) = filter_length - 1;
//根据传递函数得到：b(z)=b(1)+b(2)*z(-1)+b(3)*z(-2)+……+b(N+1)*z(-N)
//N：滤波器阶数
//Wc：归一化截止频率（0－1）对应（0－PI），已经归一化
//Type：滤波器类型：BP、BS
void CFIRDesign::FIR(int N, double Wl, double Wh, int Type, int Window)
{
	//得到滤波器阶数N并创建数据存储数组
	ASSERT(N % 2 == 0); //保证N值为偶数（目前的滤波器设计算法只能在阶数为偶数情况下有效）
	if(m_pFIR)
		delete [] m_pFIR;

	this->m_nFilterOrder = N;
	this->m_pFIR = new double [N+1];
	//窗函数设计
	this->Window(N+1,Window);
	//得到滤波器冲击响应函数
	int i;
	switch(Type)
	{
	//BandPass,带通滤波器
	case BP:
		for(i=0; i<=N; i++)
		{
			m_pFIR[i] = ( sin((i-N/2)*Wh*PI) - sin((i-N/2)*Wl*PI) ) / (PI*(i-N/2)) * m_pWindow[i];
			if(i == N/2)
				m_pFIR[i] = Wh - Wl;
		}
		break;
	//BandStop,带阻滤波器
	case BS:
		for(i=0; i<=N; i++)
		{
			m_pFIR[i] = ( sin((i-N/2)*Wl*PI) - sin((i-N/2)*Wh*PI) ) / (PI*(i-N/2)) * m_pWindow[i];
			if(i == N/2)
				m_pFIR[i] = 1 + Wl - Wh;
		}
		break;
	default:
		break;
	}
}

//滤波器频率相应曲线图（Frequency Respose）
//幅频响应：Amplititude Response; 相频响应：Phase Response
void CFIRDesign::Freqz(double* pAR, double* pPR, int nLen)
{
	// 傅立叶变换的实部(Real)和虚部(Image)数组，长度等于x(n)的长度
	double* pXR,*pXI;
	pXR = new double [nLen];
	pXI = new double [nLen];
	double W = 2*PI/nLen; // W = 2*PI/N
	
	double Real,Image;

	for(int k=0; k<nLen; k++)
	{
		Real = 0;
		Image = 0;
		for(int n=0; n<nLen; n++)
		{
//			pXR[k] += m_pFIR[n]*cos(W*n*k);
//			pXI[k] += m_pFIR[n]*sin(W*n*k);
			Real += m_pFIR[n]*cos(W*n*k);
			Image += m_pFIR[n]*sin(W*n*k);
		}
//		pAR[k] = 10*log10(sqrt(pXR[k]*pXR[k]+pXI[k]*pXI[k]));
		pAR[k] = 10*log10(sqrt(Real*Real+Image*Image));
	}

	delete [] pXR;
	delete [] pXI;
}

//用设计的滤波器进行滤波：pSingal－待滤波信号指针；Length－待滤波信号长度
// y(n) = x(n)*h(n)  （ * 代表卷积）：length(y) = lengyh(x)+length(h)-1
//    jw		jw     jw
// Y(e   ) = X(e   )H(e  )
// 注意：此处滤波后将已经处理的信号放回原始数据的数组中，故其长度为原始数组长度
void CFIRDesign::Filter(double *pSignal, int Length)
{
	ASSERT(Length>0);	//保证待处理信号长度>0
	double* pTemp = new double[Length];
	::memcpy(pTemp,pSignal,sizeof(double)*Length);//待处理信号暂存于pTemp数组中

//		double temp[41];
//		::memcpy(temp,m_pFIR,8*41);
	int	n,k,xMin,xMax;
	double convValue;
	for(n=0; n<Length; n++)
	{
		xMin = __max(0,n-m_nFilterOrder);
		xMax = __min(n,Length-1);
		convValue = 0;	//清零
		for(k=xMin; k<=xMax; k++)
		{
			convValue += pTemp[k] * m_pFIR[n-k];
		}
		pSignal[n] = convValue;
	}

	delete [] pTemp;	//释放内存

}