希尔伯特变换c程序例子 [复制链接]

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"FFT1.h"
const int L=4;//信号长度
//x初始序列，也存变换结果序列
void HT(float x[],int num)
{
    float *xi,*x1r, *x1i, *sgn;
    int i,k;
    k=log(num)/log(2);
                     xi=(float *)calloc(num,sizeof(float));
    x1r=(float *)calloc(num,sizeof(float));
    x1i=(float *)calloc(num,sizeof(float));
    sgn=(float *)calloc(num,sizeof(float));
    //apply FFT to x[]
    fft(x,xi,k,1);
                     //construct sgn[]
    for( i=1; i<num/2; i++)
        sgn=1.0;
    for(i=num/2+1; i<num; i++)
        sgn=-1.0;
    //X(k)sgn(k)
    for(i=0; i<num; i++)
    {
        x1r=x*sgn;
        x1i=xi*sgn;
    }
    //apply IFFT to x1
    fft(x1r,x1i,k,-1);
    for(i=0; i<num; i++)
          x= x1i;
                      free(x1r);
    free(x1i);
    free(sgn);
    free(xi);
}
void main()
{
    float x[L]={1,2,3,4};
    float xx[L];
    int i;
    for(i=0; i<L; i++)
            xx=x;
    HT(x,L);
    printf("Hilbert变换得到的解析信号为:\n");
    for(i=0; i<L; i++)
             printf("%f+i(%f)\n",xx,x);
}

挺好的 比较简单 适合学习

好资料
谢谢

不写程序好多年了！

学习了。

#include"FFT1.h"
fft()
这两个程序呢？

"FFT1.h" 就是快速傅里叶变换的函数：

//sr-存初始序列，也存结果的实部;inv=1,正变换,-1,反变换;
//反变换后sr存实序列
void fft(float sr[],float sx[],int m0,int inv)
{
int i,j,lm,li,k,lmx,np,lix,mm2;
double t1,t2,c,s,cv,st,ct;
double PI=3.1415926;
if(m0<0)
  return;
lmx=1;
for(i=1;i<=m0;++i)
        lmx+=lmx;   //form 2**m0
    cv=2.0*PI/(double)lmx;
ct=cos(cv); st=-inv*sin(cv);
np=lmx;mm2=m0-2;
/* fft butterfly numeration */
for(i=1;i<=mm2;++i)
{
  lix=lmx;lmx/=2;
  c=ct;s=st;
  for(li=0;li<np;li+=lix)
  {
   j=li;k=j+lmx;
   t1=sr[j]-sr[k];t2=sx[j]-sx[k];
   sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2;
   ++j;++k;
   t1=sr[j]-sr[k];t2=sx[j]-sx[k];
   sr[j]+=sr[k];sx[j]+=sx[k];
   sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
  }
  for(lm=2;lm<lmx;++lm)
  {
   cv=c;c=ct*c-st*s;s=st*cv+ct*s;
   for(li=0;li<np;li+=lix)
   {
    j=li+lm;k=lmx+j;
    t1=sr[j]-sr[k];t2=sx[j]-sx[k];
    sr[j]+=sr[k];sx[j]+=sx[k];
    sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
   }
  }
  cv=ct;ct=2.0*ct*ct-1.0;st=2.0*st*cv;
}
/* 4 points DFT */
if(m0>=2)
for(li=0;li<np;li+=4)
{
  j=li;k=j+2;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];
  sx[j]+=sx[k];
  sr[k]=t1;sx[k]=t2;
  ++j;++k;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];sx[j]+=sx[k];
  sr[k]=inv*t2;sx[k]=-inv*t1;
}
/* 2 points DFT */
for(li=0;li<np;li+=2)
{
  j=li;k=j+1;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];sx[j]+=sx[k];
  sr[k]=t1;sx[k]=t2;
}
/* sort according to bit reversal */
lmx=np/2;j=0;
for(i=1;i<np-1;++i)
{
  k=lmx;
  while(k<=j)
  {
   j-=k;k/=2;
  }
  j+=k;
  if(i<j)
  {
   t1=sr[j];sr[j]=sr;sr=t1;
   t1=sx[j];sx[j]=sx;sx=t1;
  }
}
/*  if Inverse FFT, multiply 1.0/np  */
if(inv!=-1)
  return;
t1=1.0/np;
for(i=0;i<np;++i)
{
  sr*=t1;sx*=t1;
}
}

学习了，谢谢楼主分享