Hilbert变换C源码

/******************************************/ 
/*   main.c 
/******************************************/ 
  
#include <stdio.h> 
#include <math.h> 
  
#define N 500 
  
#define sq(X) ((X)*(X)) 
  
int     main() 
{ 
        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[i1] = 0.; 
                else if (x < 1.) 
                        delta[i1] = sqrt(1. - sq(x)); 
                else 
                        delta[i1] = 0.; 
        } 
  
        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.); 
                (void) printf("%.3f %.3f %.3f %.3f\n", xmin + i * h, delta[i], k 
ppa[i], y); 
        } 
  
        return 0; 
} 
  
/******************************************/ 
/*   hilbert.c 
/******************************************/ 
  
void hilbert(int, double[], double[]); 
  
/******************************************/ 
/*   hilbert.c 
/******************************************/ 
#include <stdio.h> 
  
#include "hilbert.h" 
  
#define PI 3.14159265 
  
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; 
        } 
} 
/*********************************************** 
/*   说明 
/*********************************************** 
  
The Hilbert transform 
This package performs a numerical Hilbert transform. 
Download Compressed tar-file 
or 
hilbert.c and hilbert.h 
  
  
Header file 
The program must include the header file 
#include "hilbert.h" 
/*********************************************** 
/*   说明 
                        kappa[i1] -= 0.5 * (d1-d2) / i2 + 0.5 * (d3 - d4) / (i2 
 1); 
                } 
                kappa[i1] /= PI; 
        } 
} 
/*********************************************** 
/*   说明 
/*********************************************** 
  
The Hilbert transform 
This package performs a numerical Hilbert transform. 
Download Compressed tar-file 
or 
hilbert.c and hilbert.h 
  
  
Header file 
The program must include the header file 
#include "hilbert.h" 
  
and then call hilbert() to perform the transform. 
  
Usage 
If delta and kappa are arrays of n doubles, both arrays are allocated by the mai 
 program. 
The input data are in delta and the call hilbert(n, delta, kappa) will return th 
 Hilbert transform of delta in kappa. 
  
n and delta are not modified. 
  
Compilation 
Compile the package using a C-compiler. 
If you use a makefile and this makefile looks like this 
  
a.out: a.o b.o c.o 
       cc a.o b.o c.o 
... 
a.o: a.c d.h 
       cc -c a.c 
... 
  
and you use the package in a.c, change the makefile to look like this 
a.out: a.o b.o c.o hilbert.o 
       cc a.o b.o c.o hilbert.o 
... 
a.o: a.c d.h hilbert.h 
       cc -c a.c 
... 
hilbert.o: hilbert.c hilbert.h 
       cc -c hilbert.c 
  
Demo 
main.c uses the package to perform a Hilbert transform of a semi-ellipsis. 
Output from the program is a table. 
  
First column: x 
Second column: The semi-ellipsis 
Third column: The Hilbert transform calculated by hilbert(). 
Fourth column: The analytical Hilbert transform of the semi-ellipsis.