hilert.cpp

利用快速傅立叶变换间接计算的Hilbert变换。程序已测试

/*   main.c 
/******************************************/ 
  
#include <stdio.h> 
#include <math.h> 
  
#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<n)? delta[i1+i2]: 0.; 
                        d2 = (i1-i2>=0)? delta[i1-i2]: 0.; 
                        d3 = (i1+i2+1<n)? delta[i1+i2+1]: 0.; 
                        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; 
        } 
}