denoise.c

基于ccs2.0的小波消噪程序

#include "f2407_c.h"
#include "stdlib.h"
#include "stdio.h"
#include "wavelet_signal.h"
#include "math.h"

#define N0 64


double g0[8]={0.230378,0.714847, 0.630881, -0.027984 ,-0.187035, 0.030841, 0.032883, -0.010597};
double g1[8]={-0.010597, -0.032883, 0.030841, 0.187035, -0.027984, -0.630881, 0.714847, -0.230378};
double h0[8]={-0.010597, 0.032883, 0.030841, -0.187035, -0.027984, 0.630881, 0.714847, 0.230378};
double h1[8]={-0.230378, 0.714847, -0.630881, -0.027984, 0.187035, 0.030841, -0.032883, -0.010597};

double data_denoise[N0];

void system_init();

void inline Disable()
{
   asm(" setc INTM");
}                    
void inline Enable()
{
   asm(" clrc INTM");
}

int sign(double x)
{
	int m;
	if(x>0)
		m=1;
	else if(x==0)
		m=0;
	else
		m=-1;
	
	return(m);		
}

int signal_extension( double * original, int n, double * extended, int m)
										// original为原始信号,n为原始信号长度
										// extend为延拓信号,m为滤波器系数数组长度
{
	for( i = 0; i < m-1; i++ )
		extended[(m-1)-i-1] = original[i];
	for ( i = 0; i < n; i++ )
		extended[m-1+i] = original[i];
	for ( i = 0; i < m-1; i++ )
		extended[m-1+n+i] = original[i];
	
	return(2*m+n-2);		
}  

//  卷积计算  data为原始数据，n为数据长度，filter为滤波器系数，
//  m为滤波器系数数组长度, result为卷积结果,返回值为卷积之后的长度（m＋n－1）
 
int convolution(double *data,int n,double *filter,int m,double *result)
{
	int mn;
	double *data1=NULL;
	double sum;
	int i, j;

	mn=m+n-1;

	data1=malloc(mn*sizeof(double)); 
  
	
	if(data1 == NULL) 
		return(0);

	for(i=0;i<m-1;i++)
		data1[i]=0;
 	for(j=0;j<n;j++)
		data1[m-1+j]=data[j];
	
	
	if(mn>0)
	{
		for(i=0;i<mn;i++)
		{
			sum=0;
			for(j=0;j<m;j++)
				sum=data1[(i+j)%mn]*filter[m-1-j]+sum;
			result[i]=sum;
		}
	
	}
	
	free(data1);             //  

	return(mn);
}

//  二抽取
int downsampling(double *data,int n,double *downdata)
{
	int i,m;
	
	m=n/2;
	
	for(i=0;i<m;i++)
	{
		downdata[i]=data[i*2];
	}
	
	return(m);
}

//  二插值
int interpolation(double *data,int n,double *intpdata)
{
	int i;

	for(i=0;i<n;i++)
	{
		intpdata[2*i]=data[i];
		intpdata[2*i+1]=0;
	}

	return(2*n);
}

int wkeep(double *data,int n,int len,double *result)
{
	int i,m;
	
	m=(n-len)/2;
	
	for(i=0;i<len-1;i++)
	{
		result[i]=data[m+i];
	}
	
	return(len);
}

//  一级小波分解
int wavelet_decomposition(double *data,int n,   // 待分解信号
					  	   double *Hi,int n_Hi,  // 高通滤波器系数
				      	   double *Lo,int n_Lo,  // 低通滤波器系数
				      	   double *approximation, // 小波分解平滑信号
				      	   double *detail        // 小波分解细节信号
				       	  )       
{
	int m,mn;
	int i,k;
	double *approximation1,*detail1,*approximationtemp, *detailtemp;
	
	mn=signal_extension(data,n,data1,n_Hi);

	approximation1=malloc(mn*sizeof(double));
	detail1=malloc(mn*sizeof(double));
	approximationtemp=malloc(mn*sizeof(double));
	detailtemp=malloc(mn*sizeof(double));
	
	m=convolution(data1,mn,Hi,n_Hi,detail1);
	for(i=0;i<m-2*n_Hi-2;i++)
	{
		detailtemp[i]=detail1[n_Hi+i-1];
	}
	
	m=convolution(data1,mn,Lo,n_Lo,approximation1);
	for(i=0;i<m-2*n_Lo-2;i++)
	{
		approximationtemp[i]=approximation1[n_Lo+i-1];
	}
	
	k=downsampling(detailtemp,m-2*n_Hi-2,detail);
	k=downsampling(approximation1,m-2*n_Lo-2,approximation);
	
	free(approximation1);
	free(detail1);
	free(approximationtemp);
	free(detailtemp);
	
	return(k);
}

//  一级小波重建

int wavelet_synthesis( double *approximation, // 小波重建平滑信号
				       double *detail,int n, // 小波重建细节信号,n为信号长度
					   double *Hi,int n_Hi,  // 高通滤波器系数
				       double *Lo,int n_Lo,  // 低通滤波器系数
				       double *data,          // 小波重建信号
				       )   
{
	int m,i;
	double *approximation1, *detail1, *data1;


	approximation1=malloc(2*n*sizeof(double));
	detail1=malloc(2*n*sizeof(double));
	data1=malloc(2*n*sizeof(double));
	m=interpolation(approximation,n,approximation1);
	m=interpolation(detail,n,detail1);

	

	m=convolution(detail1,m,Hi,n_Hi,data);
	m=convolution(approximation1,m,Lo,n_Lo,data1);

	
	for(i=0;i<m;i++)
		data[i]=data1[i]+data[i];
	
	free(approximation1);
	free(detail1);
	free(data1);

	return(m);
}

void main()
{
	int m,i,j,k;
	int level;
	int nn,n;
	double resultdata[N0];
	double temp;
	double threshold;
	
	system_init();
	
	
	level=log(N0)/2+1;
	
	int *length;
	double *l, *h, *data,*result;
    n=N0;
                   
	l=malloc(n*sizeof(double));                  // 小波分解的平滑信号
	h=malloc(n*sizeof(double));                  // 小波分解的细节信号
	data=malloc(n*sizeof(double));               
	result=malloc((n+levle*14)*sizeof(double));
	
	for(i=0;i<N0;i++) 
		l[i]=data_noise[i];
	
	nn=N0;
	for(i=0;i<level;i++)
	{
		for (j=0;j<nn;j++) 
			data[j]=l[j];
		
		length[i]=wavelet_decomposition(data,nn,h1,8,h0,8,l,h);
		
		for(k=0;k<length[i]-1;k++)
		{
			result[k]=l[k];
		}

		for(k=0;k<length[i]-1;k++)
		{
			result[k+length[i]]=h[k];
		}

		nn=length[i];

	}
//  小波分解	
//	m=multilevel_decomposition_1D(data_noise,N0,h1,8,h0,8,resultdata,level);

//  消噪	
	threshold=sqrt(2*log(N0));
	
	temp[0]=0;
	for(i=1;i<=level;i++)
	{
		temp[i]=length[level-i]+temp[i-1];
	}

	first[0]=length[level-1];
	for(i=1;i<=level;i++)
	{
		first[i]=first[0]+temp[i];
	}

	first[level+1]=N0;


	for(i=m;i++;i<N0)
	{
		temp=abs(resultdata[i])-threshold;
		if(temp<0)
			resultdata[i]=0;
		else
			resultdata[i]=sign(resultdata[i])*(abs(resultdata[i])-threshold);	
	}
	
//  小波重建
	j=multilevel_reconstruction_1D(resultdata,N0,g1,8,g0,8,data_denoise,level);
	
	while(1);		
}



void system_init()
{
    asm(" setc	SXM");					// 抑制符号位扩展
	asm(" clrc	OVM");					// 累加器中结果正常溢出
	asm(" clrc	CNF");					// B0被配置为数据存储空间
	asm(" setc	INTM");					// 禁止所有中断 
	   
	SCSR1 = 0x8142;					    // CLKIN=10M，CLKOUT=4*CLKIN=40M,使能SCICLK
	WDCR = 0x0E8;				        // 不使能看门狗		      						    
	IMR = 0x0010;					    // 开中断优先级1           
	IFR = 0x0FFFF;					    // 清除全部中断标志，"写1清0"
}


void interrupt nothing() 
{
 	return;
}