文章目錄

• 希爾伯特變換
• • 連續時間信號的希爾伯特變換
  • • 離散時間信號的希爾伯特變換
    • C程序
    • matlab程序
    • 結果圖

希爾伯特變換

連續時間信號的希爾伯特變換

定義：給定一連續的時間信號x(n)，其希爾伯特變換$x(t)$定義為

$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)?\frac{1}{pi?t}$

$x(t)$可以看成是x(t)通過一濾波器的輸出，該濾波器的單位沖激響應$h(t)=\frac{1}{\pi ?t}$

離散時間信號的希爾伯特變換

定義：設離散時間信號x(n)的希爾伯特變換是$x(n)$,希爾伯特變換器的單位抽樣響應為h(n)，由連續時間信號希爾伯特變換的性質及H（jΩ）和$H（e^{jw}）$的關系，得到

因此

$h(n)=\frac{1}{2?\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}$

及

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

則x(n)的解析信號

$z(n)=x(n)+jx(n)$

也可用DFT方便地求出一個信號x(n)的解析信號及希爾伯特變換，步驟如下

對輸入信號x(n)做 DFT，得X(k)，k=0，…，N-1，注意k=\frac{N}{2},…,N-1對應負頻率

令

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

可求得希爾伯特變換$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);figurestem(n,y)

結果圖

c的結果圖

matlab的結果

對比的結果差異不大。

總結

以上是生活随笔為你收集整理的hilbert C语言的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。