Hilbert-c.rar_Hilbert DSP实现_Hilbert变换C_c++ Hilbert_dsp hilbert_h

#include <stdio.h> #include <math.h> #include <stdlib.h> #define N 1000 struct complex {float real,imag;}; /*定义一个复数结构*/ void fft(complex*x,int n); void ifft(complex *x,int n); void hilbert(float *x,int n,float *y) /*输入信号为x，输出信号为y，信号长度为n*/ { struct compx *z; int i,k; float temp; z=(struct compx *)malloc(sizeof(struct compx)*n); for(i=0;i<n;i++) { z[i].real=x[i]; z[i].imag=0; } fft(z,n); k=n/2; z[0].real=0; z[0].imag=0; z[k].real=0; z[k].imag=0; for(i=1;i<k;i++) /*进行加权*/ { temp=z[i].real; z[i].real=-z[i].imag; z[i].imag=temp; } for(i=k+1;i<n;i++) /*进行加权*/ { temp=-z[i].real; z[i].real=z[i].imag; z[i].imag=temp; } ifft(z,n); for(i=0;i<n;i++) x[i]=z[i].real; free(z); } void fft(){ int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i< log(size_x)/log(2) ;i++){ /*一级蝶形运算*/ l=1<<i; for(j=0;j<size_x;j+= 2*l ){ /*一组蝶形运算*/ for(k=0;k<l;k++){ /*一个蝶形运算*/ mul(x[j+k+l],W[size_x*k/2/l],&product); add(x[j+k],product,&up); sub(x[j+k],product,&down); x[j+k]=up; x[j+k+l]=down; } } } } /*快速傅里叶逆变换*/ void ifft(){ int i=0,j=0,k=0,l=size_x; complex up,down; for(i=0;i< (int)( log(size_x)/log(2) );i++){ /*一级蝶形运算*/ l/=2; for(j=0;j<size_x;j+= 2*l ){ /*一组蝶形运算*/ for(k=0;k<l;k++){ /*一个蝶形运算*/ add(x[j+k],x[j+k+l],&up); up.real/=2;up.img/=2; sub(x[j+k],x[j+k+l],&down); down.real/=2;down.img/=2; divi(down,W[size_x*k/2/l],&down); x[j+k]=up; x[j+k+l]=down; } } } change(); }