hilbert C语言 程序有问题

希尔伯特变换

程序有问题，详情参考评论

连续时间信号的希尔伯特变换

定义：给定一连续的时间信号x(n)，其希尔伯特变换$\hat{x}(t)$定义为

$\hat{x}(t)=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\frac{x(\tau )}{t-\tau }d(\tau )$

$=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\frac{x(t-\tau )}{\tau }d(\tau )=x(t)\ast \frac{1}{pi\ast t}$

$\hat{x}(t)$可以看成是x(t)通过一滤波器的输出，该滤波器的单位冲激响应$h(t)=\frac{1}{\pi \ast t}$

离散时间信号的希尔伯特变换

定义：设离散时间信号x(n)的希尔伯特变换是$\hat{x}(n)$,希尔伯特变换器的单位抽样响应为h(n)，由连续时间信号希尔伯特变换的性质及H（jΩ）和$H（e^{jw}）$的关系，得到

因此

$h(n)=\frac{1}{2\ast \pi }{\int }_{-\pi }^{\pi }H(e^{jw}{e}^{jwn}dw=\frac{1}{2\pi }{\int }_{-\pi }^{o}j{e}^{jwn}dw-2\pi {\int }^{\pi }j{e}^{jwn}dw$

求解上式得积分，可得

$h(n)=\frac{1-(-1{)}^n}{n\pi }=\{\begin{array}{ll}0,& \text{n为偶数}\\ \frac{2}{n\pi }& \text{n为奇数}\end{array}$

及

$\hat{x}(n)=x(n)\ast h(n)=\frac{2}{\pi }{\sum }_{m=-\infty }^{\infty }\frac{x(n-2m-1)}{2m+1}$

则x(n)的解析信号

$z(n)=x(n)+j\hat{x}(n)$

也可用DFT方便地求出一个信号x(n)的解析信号及希尔伯特变换，步骤如下

1. 对输入信号x(n)做 DFT，得X(k)，k=0，…，N-1，注意k=\frac{N}{2},…,N-1对应负频率

2. 令
   

3. 对Z(k)做 逆DFT，得到x(n)的解析信号z(n).
   dft:$X(k)={\sum }_{n=0}^{N-1}x(n)\ast (cos(2\ast \Pi /N\ast n\ast k)+i\ast sin(2\ast \Pi /N\ast n\ast k))$
   idft:$x(n)=\frac{1}{N}{\sum }_{k=0}^{N-1}X(k)exp(j\frac{2\pi }{N}nk)$

4. 可求得希尔伯特变换$\hat{x}(n)=-j[z(n)-x(n)]$

   C程序

   #include <stdio.h>
   #include <math.h>
   #define N 64
   #define pi  3.1415
   typedef struct
   {
       double real,imag;
   } complex;
   complex dft_out[100];
   complex dft_one[100];
   
   /*void idft(double z_r,double z_i)
   {
       int n,k;
       double a_r[100],a_i[100];
       double x_r[100],x_i[100];
       for(k=0; k<N; k++)//k循环
       {
           for(n=0; n<N; n++)//n循环
           {
   
               a_r[n]=z_r[n]*cos(2*pi/N*n*k);//实部信号
               a_i[n]=z_i[n]*sin(2*pi/N*n*k);//虚部信号
               x_r[k]+=x_r[n];
               x_i[k]+=x_i[n];//DFT后的实部，虚部相加
           }
           return x_r[k],x_i[k];
       }
   }*/
   /*void dftk(double x_r,double x_i)
   {
       int n;
       double z_r[100],z_i[100];
       for (n=0; n<N; n++)
       {
           if(n==0)
           {
               z_r[n]=x_r[n];
               z_i[n]=x_i[n];
   
           }
           if(1<=n&&n<N/2-1)
           {
               z_r[n]=2*x_r[n];
               z_i[n]=2*x_i[n];
   
           }
           if(N/2<=N&&n<=63)
           {
               z_r[n]=0;
               z_i[n]=0;
   
           }
       }
   }*/
   int main()
   {
       int n,k;
       double x_r[100],x_i[100];
       double z_r[100],z_i[100];
       double y_r[100],y_i[100];
       double s_r[100],s_i[100];
       double m_r[100],m_i[100];
       double xn,XK[100];
       double amp[100],xr[100],xi[100];
       for(k=0; k<N; k++)//DFT
       {
           for(n=0; n<N; n++)
           {
               xn=cos(n*pi/6);//原始信号
               dft_one[n].real=xn*cos(2*pi/N*n*k);//实部信号
               dft_one[n].imag=xn*sin(2*pi/N*n*k);//虚部信号
               dft_out[k].real+=dft_one[n].real;
               dft_out[k].imag+=dft_one[n].imag;//DFT后的实部，虚部相加
           }
           x_r[k]=dft_out[k].real;
           x_i[k]=dft_out[k].imag;
       }
       for (n=0; n<N; n++)//得到Z(K)
       {
           if(n==0)
           {
               z_r[n]=x_r[n];
               z_i[n]=x_i[n];
   
           }
           if(0<n&&n<N/2)
           {
               z_r[n]=2*x_r[n];
               z_i[n]=2*x_i[n];
   
           }
           if(N/2<=n&&n<64)
           {
               z_r[n]=0;
               z_i[n]=0;
   
           }
   
       }
       for(n=0; n<N; n++) //idft
       {
   
           for(k=0; k<N; k++)
           {
               m_r[n]+=z_r[k]*cos(2*pi/N*n*k)-z_i[k]*sin(2*pi/N*n*k);
               m_i[n]+=z_i[k]*cos(2*pi/N*n*k)+z_r[k]*sin(2*pi/N*n*k);
           }
           s_r[n]=1/N*m_r[n];
           s_i[n]=1/N*m_i[n];
       }
   
       for(n=0; n<N; n++)//输出
       {
           xr[n]=cos(n*pi/6);//原始信号
           y_r[n]=s_i[n];
           y_i[n]=xr[n]-s_r[n];
           amp[n]=sqrt(y_r[n]*y_r[n]+y_i[n]*y_i[n]);
           printf("%d    %f\n",n,y_i[n]);
       }
   }
   

matlab程序

N=64;

n=0:1:N;

n=0:1:63;

x=cos(n*pi/6);

y=hilbert(x);

figure

stem(n,y)

结果图

c的结果图

matlab的结果

对比的结果差异不大。