Hilbert的C语言实现

利用C语言实现Hilbert变换，以cos（）函数为例，对其进行hilbert变换，计算相角

/*   main.c 
/******************************************/ 
  
#include  
#include  
  
#define N 500 
#define PI 3.14159265 
  
#define sq（X） （（X）*（X）） 

void hilbert（int double[] double[]）;   

int     main（） 
{
	FILE *fp1*fp2*fp3*fp4;
	fp1=fopen（“input.txt““w+“）;
	fp2=fopen（“output.txt““w+“）;
	fp3=fopen（“phase.txt““w+“）;
	fp4=fopen（“angel.txt““w+“）;

        int             i; 
        double          x; 
        double          delta[N]; 
        double          kappa[N]; 
        double          y; 
        double          xmin = -10.; 
        double          xmax = 10.; 
        double          cd = -1.; 
        double          w = 2.; 
        double          h = （xmax - xmin） / N; 
  
        for （i = 0; i < N; i++） 
        { 
                x = 2. * （xmin + i * h - cd） / w; 
                if （x < -1.） 
                        delta[i] = 0.; 
                else if （x < 1.） 
                        delta[i] = sqrt（1. - sq（x））; 
                else 
                        delta[i] = 0.; 
//				fprintf（fp1“%.3f  “delta[i]）;
        } 
		for （i = 0; i < N; i++）
		{
			delta[i]=cos（4.*PI*i/N）;
			fprintf（fp1“%.3f  “delta[i]）;
		}
  
        hilbert（N delta kappa）; 
  
        for （i = 0; i < N; i++） 
        {
                x = 2. * （xmin + i * h - cd） / w; 
                if （x < -1.） 
                        y = x + sqrt（sq（x） - 1.）; 
                else if （x < 1.） 
                        y = x; 
                else 
                        y = x - sqrt（sq（x） - 1.）; 

				fprintf（fp2“%.3f  “kappa[i]）;
				fprintf（fp3“%.3f  “xmin + i * h）;
				fprintf（fp4“%.3f  “y）;

                （void） printf（“%.3f %.3f %.3f %.3f\n“ xmin + i * h delta[i] kappa[i] y）; 
        } 
  
        return 0; 

		fclose（fp1）;
		fclose（fp2）;
		fclose（fp3）;
		fclose（fp4）;
} 
  
/******************************************/ 
/*   hilbert.c 
/******************************************/ 
  

  
/******************************************/ 
/*   hilbert.c 
/******************************************/ 

void hilbert（int n double delta[] double kappa[]） 
{ 
        double d1 d2 d3 d4; 
        int i1 i2; 
  
        for （i1 = 0; i1 < n; i1++） 
        { 
                kappa[i1] = 0.; 
                for （i2 = 1; i2 < n; i2++） 
                { 
                        d1 = （i1+i2                        d2 = （i1-i2>=0）? delta[i1-i2]: 0.; 
                        d3 = （i1+i2+1                        d4 = （i1-i2-1>=0）? delta[i1-i2-1]: 0.; 
  
                        kappa[i1] -= 0.5 * （d1-d2） / i2 + 0.5 * （d3 - d4） / （i2+1）; 
                } 
                kappa[i1] /= PI; 
        } 
} 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        2830  2005-03-03 10:57  Hilbert.cpp