(三）：初步寻找区分病人的最佳高频/低频阈值

关注

原创最新推荐文章于 2026-09-26 10:29:00 发布 1.5k 阅读 AI 写同款· GEO 优化 ›

信号处理进阶：小波阈值去噪在医学影像中的应用实践本文深入探讨了小波阈值去噪技术在医学影像处理中的核心应用。针对传统滤波方法易丢失诊断细节的痛点，文章详细解析了小波变换的多尺度分析优势，并重点阐述了阈值选择、噪声估计及不同阈值函数的实战策略。通过结合CT、MRI等具体场景，为如何在抑制噪声与保留关键病理特征之间取得最佳平衡提供了专业指导。阅读详情

前两步，我已经做好了将心音信号转化为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;