前两步,我已经做好了将心音信号转化为txt文件,和基二FFT。
对于病人和健康人,其心音频率成分是不同的。病人的高频成分多一些。我们可以选择一个频率阈值X:大于X的频率称为高频HP,小于X的频率称为低频LP;由于病人和健康人的频率成分不同,必然存在一个最佳阈值,使得以此为阈值时病人和健康人的HP/LP差别最大。例如假设最佳频率阈值为X,此时求得病人和健康人的HP/LP,设定一标准值Y,我们说:HP/LP大于Y的都是病人,HP/LP小于Y的都是健康人。这样我们就完成了由心音频率识别病人和健康人的工作。
当然。具体的问题远没有这么简单,X和Y也不是这么容易就确定的,尤其是Y,所以我们这里只能初步确定X,怎么确定X?
我们这里用的判断方法是ROC曲线法。具体可参考:ROC曲线-阈值评判标准和ROC曲线-维基百科。
首先,我们知道了X的大体范围80-800HZ,我们以10为步长,逐一递加比较,选择最优的频率阈值。对于某一特定的频率阈值,怎么判断指标?假设我们有100个病人,100个健康人,我们选定频率阈值为100hz,对于某一个值Y1,我们找出HP/LP值大于Y1的病人个数,找出HP/LP大于Y1的健康人个数,这就确定了ROC曲线中的一个点。我们使用的是用ROC曲线的面积判断100hz的好坏,所以我们还需要好多点来确定面积。所以:我们使用[Y1,Yn],步长为y 个值来确定。
原理就是这样,下面是c++程序了。编程水平是个硬伤,我写了一遍程序,改了一遍,到最后还是发现问题不少,我如果再写一遍,估计能改进不少,不过我实在是不愿意再改了,等水平提高了,把这一整个小项目的代码都改了。
下面是我使用的命名规则:成员函数和变量使用m_X;参数类使用_X,首字母大写;函数的参数使用__x,变量使用__x;(后来我发现,虽然STL这样命名非常漂亮,不过我的编程水平却让它们非常难看。。
首先是矫正我以前抄的一个程序(得出结论:抄的程序永远不是自己的。。):
//遍历文件夹及其子文件夹下所有文件
#ifndef READ
#define READ
#include<string>
#include<io.h>
template<class _Functional>
bool _travel_folder(const std::string& _foldername,_Functional _op)
{
_finddata_t _file_info;
long _handle=_findfirst(_foldername.c_str(),&_file_info);
if(_handle==-1){
return false;
}
do
{
_op(std::string(_file_info.name));
}while(_findnext(_handle,&_file_info)==0);
return true;
}
template<class _Functional>
bool _travel_directory(const std::string& _dirname,_Functional _op)
{
_finddata_t _file_info;
std::string _findname=_dirname+"\\*";
long _handle=_findfirst(_findname.c_str(),&_file_info);
if(_handle==-1){
return false;
}
do
{
//判断是否有子目录
if(_file_info.attrib& _A_SUBDIR){
if( (strcmp(_file_info.name,".")!=0)&&(strcmp(_file_info.name,"..")!=0)){
std::string _newpath=_dirname+"\\"+_file_info.name;
_travel_directory(_newpath,_op);
}
}else{
_op(_dirname+"\\"+_file_info.name);
}
}while(_findnext(_handle,&_file_info)==0);
_findclose(_handle);
return true;
}
#endif
以前的程序,没有普遍性,不符合c++风格。
下面是类的设计:
/*关于这个类,首先我们要解决的问题是,给定两个样本组,我们设定某一评判办法,使得计算出ROC曲线中的一个点*/
#include<vector>
#include<algorithm>
#include<utility>
#include<functional>
#include<iterator>
#include"FFT.h"
#include"read.h"
using std::pair;
using std::vector;
//样本类,主要维持一个样本的特性,提供一个方法初始化这个特性
template<class _Feature,class _Argument>
class Sample
{
public:
_Feature m_feature;
public:
Sample(_Argument __arg){}
virtual bool m_init_feature(_Argument& __arg){return true;}
virtual Sample* m_clone() const{return new Sample(*this);}//支持句柄类
virtual _Feature* m_getFeature(){return &m_feature;}
Sample(){}
};
我们的目的是设计这样一个程序:有样本类,有样本组类:里面装着样本;样本最起码的接口是:构造函数,不过构造函数是什么参数还不确定;还有就是样本类维护的是样本的特性,特性是什么不确定,可能是一个数值,也可能是一个字符串;m_init_feature函数来构造样本的特性;
由于样本要放在样本组里,那么我们默认使用vector来维持一个样本组,容器与多态的冲突,导致我们使用句柄类。
//我们要定义一个样本组,使用vector容器,由于对象不是多态的,多态行为不能通过容器来实现,所以我们只能在vector中
//保存Sample*,于是我们需要维护这个指针(主要是保证解决内存问题),可以使用的一种方式是:auto_ptr,但是auto_ptr
//的复制行为不适合用于容器,所以我们需要一个管理指针的句柄类:
//首先定义一个泛型句柄:
template<class _Base>
class Handle
{
public:
Handle(_Base* __p=0):m_ptr(__p),m_use(new size_t(1)){}
_Base& operator*();
_Base* operator->();
/*
const _Base& operator*()const;
const _Base* operator->()const;
*/
Handle(const _Base& __base):m_ptr(__base.m_clone()),m_use(new std::size_t(1)){}
Handle(const Handle& __h):m_ptr(__h.m_ptr),m_use(__h.m_use){++*m_use;}
Handle& operator=(const Handle& __h);
~Handle(){m_rem_ref();}
private:
_Base* m_ptr;
size_t* m_use;
void m_rem_ref(){if(--*m_use==0){delete m_ptr;delete m_use;}}
};
template<class _Base>
inline Handle<_Base>& Handle<_Base>::operator=(const Handle& rhs)
{
++*rhs.m_use;
m_rem_ref();
ptr=rhs.m_ptr;
use=rhs.m_use;
return *this;
}
template<class _Base>
inline _Base& Handle<_Base>::operator*()
{
if(m_ptr) return *m_ptr;
cerr<<"error";
exit();
}
template<class _Base>
inline _Base* Handle<_Base>::operator->()
{
if(m_ptr) return m_ptr;
cerr<<"error";
exit(EXIT_FAILURE);
}
下面是样本组类:
//样本组类,维持一个样本组的特性
template<class _Feature,class _Argument>
class Sample_Group
{
public:
typedef Sample<_Feature,_Argument> Sample;
typedef Handle<Sample> Handle;
//private:
public:
vector<Handle> m_group;
public:
Sample_Group(){}
template<class _Sample_Iter>
Sample_Group(_Sample_Iter __first,_Sample_Iter __last)
:groups(__first,__last){}
int m_count()const {return m_group.size();}
void m_addItem(Sample* __sample){m_group.push_back(Handle(__sample));}
};
样本组中维持着一个vector,里面装的是样本的句柄,样本组提供添加样本功能,删除样本功能(我没有写。。现在才发现 ),使用样本迭代器区间构造样本组;样本组的大小;
对于ROC曲线,我们假设最常用的算法是返回ROC中的一点:
//Roc算法,输入两个组别,和一个判别的基准,返回TPR和FPR组,即返回ROC曲线中的一点
template<class _Group,class _Function>
pair<float,float> roc_point(_Group& __good,_Group& __bad,_Function __op)
{
int __P=__good.m_count();
int __TP=__op(__good);
int __N=__bad.m_count();
int __FP=__op(__bad);
return pair<float,float>((float)__TP/__P,(float)__FP/__N);
}
现在程序的大体框架完成了。下面是针对我们这个题目的具体解答:
(1)子类化样本类:我们从一开始的wav文件转化成的txt文件中构造样本,这个样本维持的特性是以80HZ:10HZ:800HZ为频率阈值,的HP/LP的值:
//对于我们具体的应用,我们每一个样本特性都是心音信号的以某个值为阈值的高频/低频
//就此问题而言,我们求出以80::10::800hz为阈值的高频/低频
class Wav_People:public Sample<vector<float>, const std::string&>
{
public:
enum{Hz_Beg=80,Hz_Step=10,Hz_End=800};
enum{SampleRate=4000,Totalnum=120000};
Wav_People(){}
Wav_People(const std::string& __file_name){ m_init_feature(__file_name);}
bool m_init_feature(const std::string& __file_name);//虚函数重载
//辅助函数
static int m_index(int __hz){return (__hz-Hz_Beg)/Hz_Step;}
};
其中构造函数的参数是const string&格式,传入的是txt文件名;维持的特性是vector<float>
为了完成这一目的,我们首先需要将文件中的数据进行傅里叶变换,变换到频率域,所以我们需要一些辅助程序:
//辅助程序1:将wav文件转化txt文件后的数据进行fft:
template<class _Tp,class _OutputIter>//_Tp表征的是txt中数据格式
bool fft_txtfile(const std::string& __file_name,_OutputIter& __outIter)
{
std::ifstream in(__file_name.c_str());
vector<_Tp> __tmp;
//由于开始测心音时数据不稳定,所以删除前五百个数据
_Tp __t;
//按理来说,应该先读取头信息,然而在这里,我们忽略头信息,因为我们现在处理的数据采样率为4000hz,数据120000个
for(int i=0;i<500;++i){
in>>__t;
}
//此处的固有数据可以变为读取值
int __newsize=truncation(Wav_People::Totalnum-500);//文件中的数据个数-我们不使用的头部的数据个数
for(int i=0;i<__newsize;++i){
in>>__t;
__tmp.push_back(__t);
}
//if(in.badbit()) return false;
in.close();
//数据读完
complex<float>* h=new complex<float>[__newsize];//好像不能用complex<_Tp>...
RFFT(__tmp.begin(),__tmp.end(),h,h+__newsize);
for(int i=0;i<__newsize;++i){
*__outIter=sqrt(h[i].module());
++__outIter;
}
delete[] h;
return true ;
}
这里输出是输出到迭代器,一开始的时候我为了输出,还写了个Out_Place类,表示输出,现在在心音转化txt中还使用的是那个。。确实不如用迭代器简练。
下面是初始化特性:
bool Wav_People::m_init_feature(const std::string& __file_name)
{
if(__file_name.substr(__file_name.size()-3,3)!="txt")
return false;
vector<float> __fft;
std::back_insert_iterator<vector<float> > __iter(__fft);
if(!fft_txtfile<float>(__file_name,__iter)) return false;
/*此时的频率信息都包含在了_fft中,我们根据频率区间的信息,产生(end-beg)/step个数据,
*其值为频率小于beg到end某一值的频率总和。首先我们需要找到_fft索引与频率的对应关系:
*[0,_fft.size())对应[0,SampleRate]hz,其中4000代表的是采样率
*所以对于_fft[i]其对应于(i/size())*SampleRate的频率,我们只使用0-SampleRate/2的频率
*我们最后产生一个数据,代表总的频率和
*/
//求解频率为x的对应的幅值,就是__fft[index],其中index=(int)(x*__fft.size/SampleRate),
int x=__fft.size();
int __num_beg=0,__num_end=0;
float __tmp=0;
for(int i=0;i<=static_cast<int>((Hz_End-Hz_Beg)/Hz_Step);++i){//m_feature中的元素个数
__num_end=static_cast<int>((Hz_Beg+i*Hz_Step)*__fft.size()/SampleRate);
for(int j=__num_beg;j<__num_end;++j){
__tmp+=__fft[j];
}
__num_beg=__num_end;
m_feature.push_back(__tmp);
}//m_feature完成
//我们需要的是高频/低频(hp/lp),现在m_feature中都是低频(lp)信息,hp/lp=(p-lp)/lp,其中p为总频率幅值和
//我们只需要在m_feature最后保存一个总的频率值即可,这么做的好处是,我们只有在真正应用时,才会计算hp/lp的
//值,所以对于某些只需要计算一部分hp/lp信息的应用来说,这种策略是最好的;
//对于我们的应用,我们会使用所有的hp/lp信息,所以采用直接计算法:
//计算总频率值:
for(int j=__num_beg;j<(int)(__fft.size()/2);++j){
__tmp+=__fft[j];
}
//此时__tmp为p
for(int i=0;i<m_feature.size();++i){m_feature[i]=__tmp/m_feature[i]-1;}
return true;
}
现在我们就能从一个文件中构造一个我们需要的样本了,那么怎么构造样本组呢?
我们用一个目录来构造一个样本组,这里只要构造出遍历目录时用到的函数对象即可:
//对于一个目录中的每一个文件,产生一个样本
class __Create_Group
{
private:
Sample_Group<vector<float>,const std::string&>* m_group;
public:
__Create_Group(Sample_Group<vector<float>,const std::string&>* __group):m_group(__group){}
void operator()(const std::string& __file_name){m_group->m_addItem(new Wav_People(__file_name));}
};
我们现在有样本组了,怎么完成我们的目的呢?首先对于每一个频率阈值对应一条ROC曲线,对于一条特定的ROC曲线,对于特定的HP/LP阈值,又对应特定的点。我们已经有寻找特定点的函数了,这个函数缺一个函数对象:
//我们知道两个组数据(特性),怎么获得ROC曲线中的一点?针对此应用,我们对于某个组,首先需要制定某一频率hz
//即:m_feature的索引值i,从每个样本中得到m_fearture[i],然后和某一标准对比,比这一标准大的为一类,小的为一类
class Roc_Func
{
private:
int m_hz;
float m_thoreshold;
public:
Roc_Func(int __hz,float __t):m_hz(__hz),m_thoreshold(__t){}
int operator()(Sample_Group<vector<float>,const std::string&>& __group);
};
int Roc_Func::operator()(Sample_Group<vector<float>,const std::string&>& __group)
{
int __result=0;
for(int i=0;i<__group.m_count();++i){
if(__group.m_group[i]->m_getFeature()->at(Wav_People::m_index(m_hz))>m_thoreshold)
++__result;
}
return __result;
}
这样我们就可以利用roc_point获得一个点了,只是我们对于这个需要HP/LP阈值:我们首先找出好样本中HP/LP和坏样本中HP/LP中的最大值max和最小值min,假设我们需要ROC曲线的N个点,那么阈值设定:min:(max-min)/N:max,对于每一个HP/LP阈值,调用一次roc_point:
//针对此应用,我们需要ROC曲线的面积作为评判标准,所以我们需要特殊的计算ROC曲线面积的方法.
//假设取曲线上PointNumber个点计算面积
const int PointNumber=10;
typedef Sample_Group<vector<float>,const std::string&> Spe_Group;
float roc_area(Spe_Group& __good,Spe_Group& __bad,int __hz)
{
//寻找最小值和最大值:
float __max,__min;
__max=__min=__good.m_group[0]->m_getFeature()->at(Wav_People::m_index(__hz));
for(int i=1;i<__good.m_count();++i){
float __num=__good.m_group[i]->m_getFeature()->at(Wav_People::m_index(__hz));
if(__num<__min){
__min=__num;
}else if(__num>__max){
__max=__num;
}
}
for(int i=1;i<__bad.m_count();++i){
float __num=__bad.m_group[i]->m_getFeature()->at(Wav_People::m_index(__hz));
if(__num<__min){
__min=__num;
}else if(__num>__max){
__max=__num;
}
}
float __result=0;
float __step=(__max-__min)/(PointNumber-1);
for(;__min<=__max;__min+=__step){
__result+=roc_point(__good,__bad,Roc_Func(__hz,__min)).first;
}
return __result;
}
这里求ROC面积,我使用了比较简单的方法,根据微积分思想,面积等于各点纵坐标之和。
最后一步:寻找最佳阈值:对于每个HZ,都有ROC曲线,找出ROC曲线最大值对应的HZ:
int hz_threshold(Spe_Group& __good,Spe_Group& __bad)
{
float __max_area=0;
int __result;
float __tmp;
for(int __hz=Wav_People::Hz_Beg;__hz<=Wav_People::Hz_End;__hz+=Wav_People::Hz_Step){
__tmp=roc_area(__good,__bad,__hz);
if(__tmp>__max_area){
__max_area=__tmp;
__result=__hz;
}
}
return __result;
}
使用程序:
Sample_Group<vector<float>,const std::string&>* group=new Sample_Group<vector<float>,const std::string&>;
_travel_directory("E:\\XXX",__Create_Group(group));
Sample_Group<vector<float>,const std::string&>* group2=new Sample_Group<vector<float>,const std::string&>;
_travel_directory("E:\\XXX2",__Create_Group(group2));
int hz=hz_threshold(*group,*group2);
std::cout<<hz;