signalprocess.cpp

FFT和Hilbert变换工具包

// SignalProcess.cpp : Defines the entry point for the DLL application.
//

#include "stdafx.h"
#include "SignalProcess.h"

BOOL APIENTRY DllMain( HANDLE hModule, 
                       DWORD  ul_reason_for_call, 
                       LPVOID lpReserved
					 )
{
    return TRUE;
}
/*-----------------------------------------------  
           FFT函数
    data:指向数据序列地指针
	a   :指向data的DFT序列的指针
	L   :2的L次方为FFT的点数
--------------------------------------------------*/
int fft(int *data,complex <double> *a,int L)
{
	complex <double> u;
	complex <double> w;
	complex <double> t;
	unsigned n=1,nv2,nm1,k,le,lei,ip;
	int i,j,m,number;
	double tmp;
	n<<=L;
	for(number = 0; number<n; number++)
	{
		a[number] = complex <double> (data[number],0);
	}
	nv2=n>>1;
	nm1=n-1;
	j=0;
	for(i=0;i<nm1;i++)
	{
		if(i<j)
		{
			t=a[j];
			a[j]=a[i];
			a[i]=t;
		}
		k=nv2;
		while(k<=j)
		{
			j-=k;
			k>>=1;
		}
		j+=k;
	}
	le=1;
	for(m=1;m<=L;m++)
	{
		lei=le;
		le<<=1;
		u=complex<double>(1,0);
		tmp=PI/lei;
		w=complex<double>(cos(tmp),-sin(tmp));
		for(j=0;j<lei;j++)
		{
			for(i=j;i<n;i+=le)
			{
				ip=i+lei;
				t=a[ip]*u;
				a[ip]=a[i]-t;
				a[i]+=t;
			}
			u*=w;
		}
	}

	return 0;
}
/*-----------------------------------------------  
           IFFT函数
    data:指向数据序列地指针
	a   :指向data的DFT序列的指针
	L   :2的L次方为FFT的点数
--------------------------------------------------*/

int ifft(complex <double> *a,int *data,int L)
{
	complex <double> u;
	complex <double> w;
	complex <double> t;
	unsigned n=1,nv2,nm1,k,le,lei,ip;
	int i,j,m,number;
	double tmp;
	n<<=L;
	nv2=n>>1;
	nm1=n-1;
	j=0;
	for(i=0;i<nm1;i++)
	{
		if(i<j)
		{
			t=a[j];
			a[j]=a[i];
			a[i]=t;
		}
		k=nv2;
		while(k<=j)
		{
			j-=k;
			k>>=1;
		}
		j+=k;
	}
	le=1;
	for(m=1;m<=L;m++)
	{
		lei=le;
		le<<=1;
		u=complex<double>(1,0);
		tmp=PI/lei;
		w=complex<double>(cos(tmp),sin(tmp));
		for(j=0;j<lei;j++)
		{
			for(i=j;i<n;i+=le)
			{
				ip=i+lei;
				t=a[ip]*u;
				a[ip]=a[i]-t;
				a[i]+=t;
			}
			u*=w;
		}
	}
	for(number = 0; number<n; number++)
	{
		data[number] = ceil(a[number].real())/n;
		a[number] = a[number]/complex<double>(n,0);
	}
	return 0;
}

/*--------------------------------------------------------------
                     Hilbert变换函数
	data：指向信号序列的指针
	filterdata：指向包络序列的指针
	dn：信号序列的点数
----------------------------------------------------------------*/
int  __stdcall hilbert(int * data , int *filterdata,int dn)
{
	int i = 0, j = 0, k = 0,l = 0,N = 0;
	complex<double> *zk;
	int *ldata;
	l = (int)(log(dn)/log(2))+1;
	N =(int) pow(2,l);
	zk = (complex<double>*)malloc(N*sizeof(complex<double>));
	ldata = (int *)malloc(N*sizeof(int));
	memcpy(ldata,data,dn*sizeof(int));
	for(i=dn;i<N;i++)
	{
		ldata[i] = 0;
	}
	fft(ldata,zk,l);
	for(i=0;i<N;i++)
	{
		if(i>=1 && i<=N/2-1)
		{
			zk[i] = complex<double>(2,0)*zk[i];
		}
		if(i>=N/2 && i<=N-1)
		{
			zk[i]= complex<double> (0,0);
		}
	}
	ifft(zk,ldata,l);
    for(i = 0 ;i<dn;i++)
	{
		filterdata[i] = (int)sqrt(pow(zk[i].imag(),2)+pow(zk[i].real(),2));
	}
	free(zk);
	free(ldata);
	return 0;
}

int __stdcall conv(int *h,int *data,int *result,int hn,int dn)
{
	int l,i,j,k,N;
	complex<double> *hk,*datak;
	l = (int)(log(hn+dn-1)/log(2))+1;
	N =(int) pow(2,l);
	int *lh,*ldata;
	lh =(int*)malloc(N*sizeof(int));
	ldata =(int*)malloc(N*sizeof(int));
	memcpy(lh,h,hn*sizeof(int));
	memcpy(ldata,data,dn*sizeof(int));
	for(i=hn;i<N;i++)
	{
		lh[i] = 0;
	}
	for(j=dn;j<N;j++)
	{
		ldata[j] = 0;
	}
    hk = (complex <double> *) malloc(N*sizeof(complex<double>));
    datak = (complex <double> *) malloc(N*sizeof(complex<double>));
	fft(lh,hk,l);
	fft(ldata,datak,l);
	for(k=0;k<N;k++)
	{
		datak[k] = datak[k]*hk[k];
	}
	ifft(datak,result,l);

	free(lh);
	free(ldata);
	free(hk);
	free(datak);
	return 0;
}
//////////////////////////////////////////////////////////////////
int fft_f(double *data,complex <double> *a,int L)
{
	complex <double> u;
	complex <double> w;
	complex <double> t;
	unsigned n=1,nv2,nm1,k,le,lei,ip;
	int i,j,m,number;
	double tmp;
	n<<=L;
	for(number = 0; number<n; number++)
	{
		a[number] = complex <double> (data[number],0);
	}
	nv2=n>>1;
	nm1=n-1;
	j=0;
	for(i=0;i<nm1;i++)
	{
		if(i<j)
		{
			t=a[j];
			a[j]=a[i];
			a[i]=t;
		}
		k=nv2;
		while(k<=j)
		{
			j-=k;
			k>>=1;
		}
		j+=k;
	}
	le=1;
	for(m=1;m<=L;m++)
	{
		lei=le;
		le<<=1;
		u=complex<double>(1,0);
		tmp=PI/lei;
		w=complex<double>(cos(tmp),-sin(tmp));
		for(j=0;j<lei;j++)
		{
			for(i=j;i<n;i+=le)
			{
				ip=i+lei;
				t=a[ip]*u;
				a[ip]=a[i]-t;
				a[i]+=t;
			}
			u*=w;
		}
	}

	return 0;
}


int conv_f(double *h,int *data,int *result,int hn,int dn)
{
	int l,i,j,k,N;
	complex<double> *hk,*datak;
	l = (int)(log(hn+dn-1)/log(2))+1;
	N =(int) pow(2,l);
	int *ldata;
	double *lh;
	lh =(double*)malloc(N*sizeof(double));
	ldata =(int*)malloc(N*sizeof(int));
	memcpy(lh,h,hn*sizeof(double));
	memcpy(ldata,data,dn*sizeof(int));
	for(i=hn;i<N;i++)
	{
		lh[i] = 0;
	}
	for(j=dn;j<N;j++)
	{
		ldata[j] = 0;
	}
    hk = (complex <double> *) malloc(N*sizeof(complex<double>));
    datak = (complex <double> *) malloc(N*sizeof(complex<double>));
	fft_f(lh,hk,l);
	fft(ldata,datak,l);
	for(k=0;k<N;k++)
	{
		datak[k] = datak[k]*hk[k];
	}
	ifft(datak,result,l);
	free(lh);
	free(ldata);
	free(hk);
	free(datak);
	return 0;
}

/*int __stdcall firwin_e(int n,int band,int fl,int fh,int fs,int wn,int *h)
{
	int i,n2,mid;
	double s,wc1,wc2,beta,delay;
	beta=0.0;
	double fln = (double)fl / fs;
	double fhn = (double)fh / fs;
//	if(wn==7)
//	{
//		printf("input beta parameter of Kaiser window(3<beta<10)\n");
//		scanf("%lf",&beta);
//	}
	beta = 6;
	if((n%2)==0)
	{
		n2=n/2-1;
		mid=1;
	}
	else
	{
		n2=n/2;
		mid=0;
	}
	delay=n/2.0;
	wc1=2.0*PI*fln;
	if(band>=3) wc2=2.0*PI*fhn;
	switch(band)
	{
	case 1://低通
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(int)((sin(wc1*s)/(PI*s))*window(wn,n+1,i,beta));
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=(int)(wc1/PI);
			break;
		}

	case 2: //高通
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(int)((sin(PI*s)-sin(wc1*s))/(PI*s));
				*(h+i)=*(h+i)*(int)(window(wn,n+1,i,beta));
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=(int)(1.0-wc1/PI);
			break;
		}
	case 3: //带通
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(int)((sin(wc2*s)-sin(wc1*s))/(PI*s));
				*(h+i)=*(h+i)*(int)(window(wn,n+1,i,beta));
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=(int)((wc2-wc1)/PI);
			break;
		}
	case 4: //带阻
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(int)((sin(wc1*s)+sin(PI*s)-sin(wc2*s))/(PI*s));
				*(h+i)=*(h+i)*(int)(window(wn,n+1,i,beta));
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=(int)((wc1+PI-wc2)/PI);
			break;
		}
	}
	return 0;
}*/

double window(int type,int n,int i,double beta)
{
	int k;
	double w=1.0;
	switch(type)
	{
	case 1:
		{
			w=1.0;
			break;
		}
	case 2:
		{
			k=(n-2)/10;
			if(i<=k) w=0.5*(1.0-cos(i*PI/(k+1)));
			if(i>n-k-2) w=0.5*(1.0-cos((n-i-1)*PI/(k+1)));
			break;
		}
	case 3:
		{
			w=1.0-fabs(1.0-2*i/(n-1.0));
			break;
		}
	case 4:
		{
			w=0.5*(1.0-cos(2*i*PI/(n-1)));
			break;
		}
	case 5:
		{
			w=0.54-0.46*cos(2*i*PI/(n-1));
			break;
		}
	case 6:
		{
			w=0.42-0.5*cos(2*i*PI/(n-1))+0.08*cos(4*i*PI/(n-1));
			break;
		}
	case 7:
		{
			w=kaiser(i,n,beta);
			break;
		}
	}
	return(w);
}

double kaiser(int i,int n,double beta)
{
	double a,w,a2,b1,b2,beta1;
	b1=bessel0(beta);
	a=2.0*i/(double)(n-1)-1.0;
	a2=a*a;
	beta1=beta*sqrt(1.0-a2);
	b2=bessel0(beta1);
	w=b2/b1;
	return(w);
}

double bessel0(double x)
{
	int i;
	double d,y,d2,sum;
	y=x/2.0;
	d=1.0;
	for(i=1;i<=25;i++)
	{
		d=d*y/i;
		d2=d*d;
		sum=sum+d2;
		if(d2<sum*(1.0e-8)) break;
	}
	return(sum);
}


int __stdcall firwin_e(int n,int band,int fl,int fh,int fs,int wn, int *data,int *result,int dn)
{
	int i,n2,mid;
	double *h;
	double s,wc1,wc2,beta,delay;
	beta=0.0;
	double fln = (double)fl / fs;
	double fhn = (double)fh / fs;
//	if(wn==7)
//	{
//		printf("input beta parameter of Kaiser window(3<beta<10)\n");
//		scanf("%lf",&beta);
//	}
	h = (double *)malloc((n+1)*sizeof(double));
	beta = 6;
	if((n%2)==0)
	{
		n2=n/2-1;
		mid=1;
	}
	else
	{
		n2=n/2;
		mid=0;
	}
	delay=n/2.0;
	wc1=2.0*PI*fln;
//	FILE *fp;
	if(band>=3) wc2=2.0*PI*fhn;
	switch(band)
	{
	case 1://低通
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(sin(wc1*s)/(PI*s))*window(wn,n+1,i,beta);
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=wc1/PI;
			break;
		}

	case 2: //高通
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(sin(PI*s)-sin(wc1*s))/(PI*s);
				*(h+i)=*(h+i)*window(wn,n+1,i,beta);
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=1.0-wc1/PI;
			break;
		}
	case 3: //带通
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(sin(wc2*s)-sin(wc1*s))/(PI*s);
				*(h+i)=*(h+i)*window(wn,n+1,i,beta);
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=(wc2-wc1)/PI;
			break;
		}
	case 4: //带阻
		{
			for (i=0;i<=n2;i++)
			{
				s=i-delay;
				*(h+i)=(sin(wc1*s)+sin(PI*s)-sin(wc2*s))/(PI*s);
				*(h+i)=*(h+i)*window(wn,n+1,i,beta);
				*(h+n-i)=*(h+i);
			}
			if(mid==1) *(h+n/2)=(wc1+PI-wc2)/PI;
			break;
		}
	}
//	fp=fopen("h.dat","w");
//	for(int p= 0;p<n+1;p++)
//	{
//		fprintf(fp, "%f\n", h[p]);
//	}
//	fclose(fp);
	conv_f(h,data,result,n+1,dn);
	free(h);
	return 0;
}