signal.cpp

DSP图像处理程序源码

// Dsignal.cpp: implementation of the Dsignal class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "signal.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
Dsignal::Dsignal(void){
	length=1;
	sgnt=(double *)malloc(sizeof(double));
	sgnw=(complex *)malloc(sizeof(complex));
	*sgnt=0;
	FFT();
}
Dsignal::Dsignal(int N){
	length=N;
	sgnt=(double *)malloc(sizeof(double)*N);
	sgnw=(complex *)malloc(sizeof(complex)*N);
	int i;
	for (i=0;i<N;i++){
		sgnt[i]=0.0;
	}
	FFT();
}
Dsignal::Dsignal(double *xn,int N){
	length=N;
	sgnt=(double *)malloc(sizeof(double)*N);
	sgnw=(complex *)malloc(sizeof(complex)*N);
	int i;
	for (i=0;i<N;i++){
		sgnt[i]=xn[i];
	}
	FFT();
}
void Dsignal::FFT(){
	CFFT::FFT(sgnt,sgnw,length);
}
void Dsignal::IFFT(){
	CFFT::IFFT(sgnw,sgnt,length);
}
Dsignal& Dsignal::operator = (const Dsignal Orign){
	int i;
	length=Orign.length;
	(*this).fitForLen(length);
	for (i=0;i<Orign.length;i++){
		sgnt[i]=Orign.sgnt[i];
		sgnw[i]=Orign.sgnw[i];
	}
	return *this;
}
Dsignal& Dsignal::operator = (const double *Orign){
	int i;
	for (i=0;i<length;i++)
		sgnt[i]=Orign[i];
	FFT();
	return *this;
}
Dsignal& Dsignal::operator = (const complex *Orign){
	int i;
	for (i=0;i<length;i++)
		sgnw[i]=Orign[i];
	IFFT();
	return *this;
}
void Dsignal::fitForLen(int N){
	double *sgntemp;
	int i,j;
	sgntemp=(double *)malloc(sizeof(double)*length);
	for (j=0;j<length;j++)
		sgntemp[j]=sgnt[j];
	i=length;
	length=N;
	sgnt=(double *)malloc(sizeof(double)*length);
	sgnw=(complex *)malloc(sizeof(complex)*length);
	if (i<N)
	{
	for (j=0;j<i;j++){
		sgnt[j]=sgntemp[j];
	}
	for (j=i;j<N;j++){
		sgnt[j]=0.0;
	}
	}
	else 
	{
		for (j=0;j<N;j++)
		{
			sgnt[j]=sgntemp[j];
		}
	}
	FFT();
}

void Dsignal::fitForLenW(int N){
	complex *sgntemp;
	int i,j;
	sgntemp=(complex *)malloc(sizeof(complex)*length);
	for (j=0;j<length;j++)
		sgntemp[j]=sgnw[j];
	i=length;
	length=N;
	sgnt=(double *)malloc(sizeof(double)*length);
	sgnw=(complex *)malloc(sizeof(complex)*length);
	if (i<N)
	{
	for (j=0;j<i;j++){
		sgnw[j]=sgntemp[j];
	}
	for (j=i;j<N;j++){
		sgnw[j]=complex(0,0);
	}
	}
	else 
	{
		for (j=0;j<N;j++)
			sgnw[j]=sgntemp[j];
	}
	IFFT();
}

Dsignal operator+ (Dsignal s1,Dsignal s2){
	int length;
	if (s1.length>s2.length){
		length=s1.length;
		s2.fitForLen(length);
	}
	else {
		length=s2.length;
		s1.fitForLen(length);
	}
	Dsignal s3;
	s3.length=length;
	s3.sgnt=(double *)malloc(sizeof(double)*length);
	s3.sgnw=(complex *)malloc(sizeof(complex)*length);
	int i;
	for(i=0;i<length;i++)
		s3.sgnt[i]=s1.sgnt[i]+s2.sgnt[i];
	s3.FFT();
	return s3;
}

Dsignal operator- (Dsignal s1,Dsignal s2){
	int length;
	if (s1.length>s2.length){
		length=s1.length;
		s2.fitForLen(length);
	}
	else {
		length=s2.length;
		s1.fitForLen(length);
	}
	Dsignal s3;
	s3.length=length;
	s3.sgnt=(double *)malloc(sizeof(double)*length);
	s3.sgnw=(complex *)malloc(sizeof(complex)*length);
	int i;
	for(i=0;i<length;i++)
		s3.sgnt[i]=s1.sgnt[i]-s2.sgnt[i];
	s3.FFT();
	return s3;
}

Dsignal operator* (Dsignal s1,Dsignal s2){
	int length;
	if (s1.length>s2.length){
		length=s1.length;
		s2.fitForLen(length);
	}
	else {
		length=s2.length;
		s1.fitForLen(length);
	}
	Dsignal s3;
	s3.length=length;
	s3.sgnt=(double *)malloc(sizeof(double)*length);
	s3.sgnw=(complex *)malloc(sizeof(complex)*length);
	int i;
	for(i=0;i<length;i++)
		s3.sgnt[i]=s1.sgnt[i]*s2.sgnt[i];
	s3.FFT();
	return s3;
}

Dsignal operator^ (Dsignal s1,Dsignal s2){
	return Dsignal::conv(s1,s2);
}

Dsignal& Dsignal::ChangeToDB(void){
	double first;
	first=(*this).sgnt[0];
	for(int i=0;i<(*this).length;i++)
	{
		(*this).sgnt[i]=20*log10((*this).sgnt[i]/first);
	}
	return (*this);
}

Dsignal Dsignal::conv(Dsignal s1,Dsignal s2){
	int length;
	length=s1.length+s2.length-1;
	Dsignal s3;
	s3.length=length;
	s3.sgnt=(double *)malloc(sizeof(double)*length);
	s3.sgnw=(complex *)malloc(sizeof(complex)*length);
	s1.fitForLen(length);
	s2.fitForLen(length);
	s1.FFT();
	s2.FFT();
	int i;
	for (i=0;i<length;i++)
		s3.sgnw[i]=s1.sgnw[i]*s2.sgnw[i];
	s3.IFFT();
	return s3;
}

void Dsignal::print(void){
	int i;
	for (i=0;i<length;i++){
		cout<<sgnt[i]<<"    ";
		sgnw[i].print();
	}
}
bool operator==(const Dsignal s1, const Dsignal s2) {
  int flag=1;
	if (s1.length!=s2.length) flag=0;
  else {
	  int i,length;
	  length=s1.length;
	  i=0;
	  while (i<length){
		  if (fabs(s1.sgnt[i]-s2.sgnt[i])>0.00000001){
			  flag=0;
			  break;
		  }
		  i++;
	  }
	  
  }
  if (flag) return true;
  else return false;
}

Dsignal& Dsignal::GetLowpassHn(double wp,double ws,double attenuation)
{
	Dsignal hdn;
	Dsignal wn;
	wn=SelectWindowStyle(wp,ws,attenuation);
	hdn.fitForLen(wn.length);
	hdn.setStyle((wp+ws)/2.0,0);
	(*this)=hdn*wn;
	return *this;
}

Dsignal& Dsignal::GetLowpassHdn(double wp,double ws,double attenuation)
{
	Dsignal wn;
	wn=SelectWindowStyle(wp,ws,attenuation);
	(*this).fitForLen(wn.length);
	(*this).setStyle((wp+ws)/2.0,0);
	return (*this);
}

Dsignal& Dsignal::SelectWindowStyle(double wp,double ws,double attenuation)
{
	int windowstyle=0;
	double excessboard=0.9;
	double wc;
	double detaw;
	wc=(wp+ws)/2;
	detaw=ws-wp;
	int N;
	if (attenuation<=21) {windowstyle=RectangleW;excessboard=0.9;}
	else if (attenuation<=25) {windowstyle=Triangle;excessboard=2.1;}
	else if (attenuation<=44) {windowstyle=Hanning;excessboard=3.1;}
	else if (attenuation<=53) {windowstyle=Hamming;excessboard=3.3;}
	else if (attenuation<=74) {windowstyle=Blackman;excessboard=5.5;}
	N=(int)ceil(excessboard*2*M_PI/detaw);
	if (N%2==0) N++;
	(*this).fitForLen(N);
	setStyle(wc,windowstyle);
	return *this;
}

void Dsignal::setStyle(double wc,int style)
{	int n;
	double temp2;
	double temp1;
	switch(style)
	{
	case LowpassHd:
		for (n=0;n<length;n++)
		{	temp1=wc*(n-(length-1)/2.0);
			if (fabs(temp1)<0.0001)
				temp2=1.0;
			else temp2=sin(temp1)/(temp1);
			sgnt[n]=wc/M_PI*temp2;
		}
		break;
	case RectangleW:
		for (n=0;n<length;n++)
		{
			sgnt[n]=1.0;
		}
		break;
	case Triangle:
		for (n=0;n<=(length-1)/2;n++)
		{
			sgnt[n]=2.0*n/(length-1);
		}
		for (;n<length;n++)
		{
			sgnt[n]=2-2.0*n/(length-1);
		}
		break;
	case Hanning:
		for (n=0;n<length;n++)
		{
			sgnt[n]=1/2.0*(1-cos(2.0*M_PI*n/(length-1)));
		}
		break;
	case Hamming:
		for(n=0;n<length;n++)
		{
			sgnt[n]=0.54-0.46*cos(2*M_PI*n/(length-1));
		}
		break;
	case Blackman:
		for (n=0;n<length;n++)
		{
			sgnt[n]=0.42-0.5*cos(2*M_PI*n/(length-1))+0.08*cos(4*M_PI*n/(length-1));
		}
		break;
	default: break;
	}
}

Dsignal& Dsignal::FreqChangeToTime(Dsignal sgn)
{
	(*this).fitForLen(sgn.length);
	int i;
	for (i=0;i<sgn.length;i++)
	{
		(*this).sgnt[i]=sgn.sgnw[i].abs();
	}
	return (*this);
}

void Dsignal::SignalToCString(CString &str)
{
	str.Empty();
	CString strtemp;
	for (int i=0;i<length;i++)
	{
		strtemp.Format("%4.4f",sgnt[i]);
		str=str+strtemp+' ';
	}

}

Dsignal& Dsignal::CStringToSignal(CString str)
{
	CString temp;
	double array[100]={0};
	    for (int i=0,aim=0;i<str.GetLength();i++)
		{
			if (str.Mid(i,1)!=" "&&str.Mid(i,1)!="\0")
				temp+=str.Mid(i,1);
			else 
			{
				array[aim++]=atof(temp);
				temp.Empty();
			}
		}
		(*this).fitForLen(aim);
		for (i=0;i<aim;i++)
		(*this).sgnt[i]=array[i];
		(*this).FFT();
	return *this;
}

Dsignal& Dsignal::Mid(int left,int right)
{
	if (left<0)
		left=0;
	if (left>(*this).length-1)
		left=(*this).length-1;
	if (right<0)
		right=0;
	if (right>(*this).length-1)
		right=(*this).length-1;
	if (left<=right&&(*this).length>right)
	{
	Dsignal sgntemp;
	sgntemp=(*this);

	(*this).fitForLen(right-left+1);
	for (int i=0;i<right-left+1;i++)
	{
		(*this).sgnt[i]=sgntemp.sgnt[left+i];
	}
	}
	return (*this);
}

double Dsignal::maxnum(){
	double max=sgnt[0];
	double temp;
	for (int i=0;i<length;i++){
		temp=sgnt[i];
		if (temp>max) max=temp;
		
	}
	return max;
}

double Dsignal::minnum(){
	double min=sgnt[0];
	double temp;
	for (int i=0;i<length;i++){
		temp=sgnt[i];
		if (temp<min) min=temp;
		
	}
	return min;
}

Dsignal::~Dsignal(){
}


/*debug
void main(){
	
	Dsignal sgn1,sgn2,sgn3;
	sgn1.fitForLen(101);
	sgn2.fitForLen(101);
	sgn1.setStyle(M_PI/4,0);
	sgn2.setStyle(M_PI/2,5);
	sgn3=sgn1*sgn2;
	for (int i=0;i<sgn3.length/2;i++)
		cout<<sgn3.sgnw[i].abs()<<endl;
}
//debug*/