filter.~cpp

来自「脑电信号分析软件」· ~CPP 代码 · 共 488 行

//---------------------------------------------------------------------------


#pragma hdrstop

#include "Filter.h"

//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
#ifndef EEGFILTER_H
#define EEGFILTER_H

#include "Filter.h"

using namespace std;

const int APASS = 1;
const int ASTOP = 60;

#define WIN_OS

#ifdef WIN_OS

// 这是 Windows 下多线程工作的 P 操作
#define P(S) WaitForSingleObject(S, INFINITE)
// 这是 Windows 下多线程工作的 V 操作
#define V(S) ReleaseSemaphore(S, 1, NULL)

#endif


Filter::~Filter()
{
}

/*!----------------------------------------------------------------------------
 * \brief 构造 
 *
 * 限定死构造的时候必须给出下面俩个参数
 *
 * \param fs 采样率
 * \param leadCount 导联数
 *
 * ----------------------------------------------------------------------------*/ 
IIRFilter::IIRFilter(int fs, int leadCount)     //IIRFilter 构造函数
{
	this->fs = fs;
	this->leadCount = leadCount;
	//bActive = true;
    bActive = false;

#ifdef WIN_OS
	semaphore = CreateSemaphore(NULL, 1, 2, NULL);
#endif
	
} //当运行到此时认为空间开辟结束,则需要再次运行一遍(关闭空间?)

IIRFilter::~IIRFilter()
{
    RemoveAll();

#ifdef WIN_OS
	CloseHandle(semaphore);
#endif
}

/*!----------------------------------------------------------------------------
 * \brief 
 *     是否激活滤波器
 *
 * \param bActive 是否激活滤波器
 *
 * ----------------------------------------------------------------------------*/ 
void IIRFilter::SetStatus(bool bActive)
{
	this->bActive = bActive;
}

/*!----------------------------------------------------------------------------
 * \brief  
 *     返回滤波器状态
 *
 * \return 滤波器状态
 *
 * ----------------------------------------------------------------------------*/
bool IIRFilter::GetStatus()
{
	return bActive;
}

void IIRFilter::Reset()
{
	FilterDesign *pFilterDesign;
	FilterMapIterator it;

	try
	{
		for(it = filterList.begin(); it != filterList.end(); ++it)
		{
			pFilterDesign = it->second;
			pFilterDesign->Reset();
		}
	}
	catch(...)
	{
		FILTER_ERROR("复位滤波器出错！");
	}
}

/*!----------------------------------------------------------------------------
 * \brief 得到导联数
 *
 * \return 数据缓冲区的导联数
 *
 * ----------------------------------------------------------------------------*/ 
int		IIRFilter::GetLeadCount()
{
	return leadCount;
}
/*!----------------------------------------------------------------------------
 * \brief 设置导联数
 *
 * \param 缓冲区顺序存储的导联数
 *
 * ----------------------------------------------------------------------------*/ 
void	IIRFilter::SetLeadCount(int leadCount)
{
	this->leadCount = leadCount;
}

/*!----------------------------------------------------------------------------
 * \brief 获得采样率
 *
 * 采样率是实际值。 如100，200……
 *
 * \return 波形数据的采样率
 *
 * ----------------------------------------------------------------------------*/ 
int	IIRFilter::GetFs()
{
	return fs;
}
/*!----------------------------------------------------------------------------
 * \brief 设置采样率
 *
 * ----------------------------------------------------------------------------*/ 
void	IIRFilter::SetFs(int fs)
{
	this->fs = fs;
}

/*!----------------------------------------------------------------------------
 * \brief 滤波
 *
 *     对给定的缓冲区进行滤波。滤波后的数据存在所给的缓冲区
 *
 * \param buf	数据缓冲区
 * \param points 缓冲区点数  就是单导缓冲长度
 * 使用已经设置的导联数
 * 使用已经设置的采样率
 *
 * ----------------------------------------------------------------------------*/
void IIRFilter::Filtering(short int *buf, int points)
{
	Filtering(buf, points, this->leadCount);
}

/*!----------------------------------------------------------------------------
 * \brief 滤波
 *
 *     对给定的缓冲区进行滤波，滤波后的数据存在所给的缓冲区。不改变已经设置的
 * 导联数和采样率
 *
 * \param buf	数据缓冲区
 * \param points	缓冲区点数
 * \param leadCount 导联数
 * \param fs	采样率
 *
 * ----------------------------------------------------------------------------*/
void IIRFilter::Filtering(short int *buf, int points, int leadCount)
{
	//int NewFilterSampLv = fs / 100;
	//for (int  j=0; j < size*NewFilterSampLv; j++)
//    double inData = 0;
//    double temp = 0;

	try
	{
		if(bActive)
		{
			P(semaphore);
			
			for (int  j=0; j < points; j++)
			{
				//        double temp = buf[75*1200];
				for (int i=0;i<leadCount;i++)
				{
//                    inData = (double)buf[j*leadCount+ i];
					buf[j*leadCount+ i] = FilterProc(i, (double)buf[j*leadCount+ i]);
//                    temp = FilterProc(i, inData);
//                    buf[j*leadCount + i] = temp;
				}
			}

			V(semaphore);
			
		}
	}
	catch(...)
	{
		FILTER_ERROR("采样数据送入FilterProc出错！");
	}
}

/* TODO 对缓冲区中特定的导联做滤波其它不做。
   virtual void FilteringNo(short int *buf, int points, int leadNo) = 0;
   virtual void FilteringNo(short int *buf, int points, int leadNo, int leadCount) = 0;*/

/*!----------------------------------------------------------------------------
 * \brief 
 *     计算滤波器幅频曲线
 *
 * \param freq 要求幅频值的频率点
 *
 * \return 返回该点的幅频值
 *
 * ----------------------------------------------------------------------------*/
double IIRFilter::MagResponse(double freq)
{
	double range = 1.0;

	FilterMapIterator it;
	FilterDesign * pFilterDesign;

	for(it = filterList.begin(); it != filterList.end(); ++it)   //说明成功开辟内存空间,即成功设计滤波器
	{
		pFilterDesign = it->second;
		range *= pFilterDesign->MagResponse(freq);
	}

	return range;
}

void IIRFilter::AddLowPass(FilterType filterType, double fPass, double fStop,
                                    double aPass, double aStop)
{
	AddPass(filterType, fPass, fStop, aPass, aStop, FilterDesign::LOW_PASS);
}

void IIRFilter::AddLowPass(FilterType filterType, 
	double fPass, double fStop)
{
	AddPass(filterType, fPass, fStop, APASS, ASTOP, FilterDesign::LOW_PASS);
}

/* 这个需要对照表的接口先去掉。到下一层的类再说吧
void IIRFilter::AddLowPass(FilterType filterType, double fPass,
                                    double aPass, double aStop)
{
}


void IIRFilter::AddLowPass(FilterType filterType, double fPass)
{
}*/
		
void IIRFilter::AddHeightPass(FilterType filterType, double fStop, double fPass,
                            double aPass, double aStop)
{
	AddPass(filterType, fStop, fPass, aPass, aStop, FilterDesign::HEIGHT_PASS);
}


void IIRFilter::AddHeightPass(FilterType filterType,
	double fStop, double fPass)
{
	AddPass(filterType, fStop, fPass, APASS, ASTOP, FilterDesign::HEIGHT_PASS);
}

/* 这个需要对照表的接口先去掉。到下一层的类再说吧
void IIRFilter::AddHeightPass(FilterType filterType, double fPass,
                                    double aPass, double aStop)
{
}

void IIRFilter::AddHeightPass(FilterType filterType, double fPass)
{
}*/

//! \todo 处理aStop aPass
void IIRFilter::AddPass(FilterType filterType, double lowFre, double heightFre, 
	   	double aPass, double aStop, FilterDesign::PassType passType)
{
	//只处理高低通两种情况
	if(passType == FilterDesign::HEIGHT_PASS)
	{
	     	switch(filterType)    //判断寻找想要设计实现的滤波器类型,若找到,则记录(开辟内存)下该滤波器
	     	{
	     		case BUTTER:
	     			filterList["HEIGHT_PASS"] = new Butter(passType, fs, lowFre, heightFre, aPass, aStop);
				break;

			case CHEBY_ONE:
				filterList["HEIGHT_PASS"] = new ChebyOne(passType, fs, lowFre, heightFre, aPass, aStop);
				break;

			case CHEBY_TWO:
				filterList["HEIGHT_PASS"] = new ChebyTwo(passType, fs, lowFre, heightFre, aPass, aStop);
	      	    	break;

			default:
				break;
		}
	}
	else
	{
	 	switch(filterType)
		{
			case BUTTER:
				filterList["LOW_PASS"] = new Butter(passType, fs, lowFre, heightFre, aPass, aStop);
				break;

			case CHEBY_ONE:
				filterList["LOW_PASS"] = new ChebyOne(passType, fs, lowFre, heightFre, aPass, aStop);
				break;

			case CHEBY_TWO:  
				filterList["LOW_PASS"] = new ChebyTwo(passType, fs, lowFre, heightFre, aPass, aStop);
	 			break;

			default:
				break;
		}
	}

}

void IIRFilter::AddNotch(int fre, string notchName)
{
	filterList[notchName] = new Notch(FilterDesign::NOTCH, fs, fre);
}

int IIRFilter::RemoveNotch(string notchName)
{
	return RemoveByName(notchName);
}

int IIRFilter::RemoveLowPass()
{
	return RemoveByName("LOW_PASS");
}

int IIRFilter::RemoveHeightPass()
{
	return RemoveByName("HEIGHT_PASS");
}

int IIRFilter::RemoveByName(std::string name)
{
	/*for(FilterMapIterator it = filterList.begin();
			it != filterList.end();)
	{
		if(it->first == name)
		{
			delete it->second;
			filterList.erase(it);
		}
	}*/
    
	P(semaphore);

	FilterMapIterator it;
	int num = 0;

	try
	{
		if((it = filterList.find(name)) != filterList.end())
		{
			delete it->second;
			filterList.erase(it);
			num++;
		}
	}
	catch(...)
	{
		FILTER_ERROR("移除滤波器错误！");
	}

	V(semaphore);
    
	return num;
}

int IIRFilter::RemoveAll()
{
	P(semaphore);

    int num = 0;
	try
	{
		for(FilterMapIterator it = filterList.begin();
				it != filterList.end();)
		{
			delete it->second;
			filterList.erase(it++);
			num++;
		}
	}
	catch(...)
	{
		FILTER_ERROR("移除滤波器错误！");
	}

	V(semaphore);

    return num;
}


/*!----------------------------------------------------------------------------
 * \brief  对第X导的一个点进行滤波
 *
 * \param leadNo 第X导
 * \param inDatum 数据
 *
 * ----------------------------------------------------------------------------*/
double IIRFilter::FilterProc(int leadNo, double inDatum)
{
	double *pNum;
	double *pDen;
	double *pZ;
	int    order;

	double datum = inDatum;

	FilterDesign *pFilterDesign;
	FilterMapIterator it;

	try
	{
      		for(it = filterList.begin(); it != filterList.end(); ++it)  //
		{
			pFilterDesign = it->second;


			pNum = pFilterDesign->GetNum();
			pDen = pFilterDesign->GetDen();
			pZ   = pFilterDesign->GetZ();
			order = pFilterDesign->GetOrder();

                        if(pNum != NULL && pDen != NULL && pZ != NULL)
                         {
			    //datum = Transform(order, pNum, pDen, datum, pZ, false, leadNo);
       			    datum = Transform(order, pDen, pNum, datum, pZ, false, leadNo);
                         }
                         else
                          {
                                  FILTER_ERROR("传递函数出错！");
                          }
                      //  if (it != filterList.end())
		}

            }
	catch(...)
	{
		FILTER_ERROR("传递函数出错！");
	}

	return datum;
}


double IIRFilter::Transform(int N,double *Fenmu,double *Fenzi,double x,
					double *z,bool flag,int lead)
{
     double output;
    int length = 2*N;
    output = Fenzi[0]*x + z[length*lead + 0];
    for(int j = 0; j < N; ++j)
    {
        z[length*lead + j] = Fenzi[j + 1]*x + z[length*lead + j + 1] - Fenmu[j + 1]*output;
    }

    return (output);
}

#endif