前两步，我已经做好了将心音信号转化为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 ;
}下面是初始化特性：
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;
}最后一步：寻找最佳阈值：对于每个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;
扫一扫
3872
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
