希尔伯特变换原理公式

MATLAB官方帮助文档中希尔伯特变换算法

hilbert uses a four-step algorithm:
Calculate the FFT of the input sequence, storing the result in a vector x.
Create a vector h whose elements h(i) have the values:
1 for i = 1, (n/2)+1
2 for i = 2, 3, ... , (n/2)
0 for i = (n/2)+2, ... , n
Calculate the element-wise product of x and h.
Calculate the inverse FFT of the sequence obtained in step 3 and returns the first n elements of the result.

常量和结构体定义

// ConstParam.h
#ifndef CONSTPARAM
#define CONSTPARAM
#include "fftw3.h"const float PI = 3.1415926535897932;
const float TWO_PI = 6.2831853071795864;const int REAL = 0;
const int IMAG = 1;const __int64 BufSize = 1048576;float CosVal[5] = { 1,0.309016994374947,-0.809016994374947 ,-0.809016994374948,0.309016994374947 };
float SinVal[5] = { 0,0.951056516295154, 0.587785252292473, -0.587785252292473, -0.951056516295154 };const float ACC = 40.0;
const float BIGNO = 1.0e10;
const float BIGNI = 1.0e-10;#define max(a,b) (a>b?a:b)typedef struct _Resampler_
{int _upRate;int _downRate;float *_transposedCoefs;float *_state;float *_stateEnd;int _paddedCoefCount;  // ceil(len(coefs)/upRate)*upRateint _coefsPerPhase;    // _paddedCoefCount / upRate__int64 _t; // "time" (modulo upRate)__int64 _xOffset;
}sResampler;
#endif

C语言实现(FFTW库的float版，double版类似)

// hilbert.h
#include <memory.h>
#include "ConstParam.h"
/*
in - 输入数据指针
out - 希尔伯特变换结果数组指针
N - 数据长度
*/
void hilbert(float* in, fftwf_complex* out,int N)
{// copy the data into the complex arrayfor (int i = 0; i < N; ++i){out[i][REAL] = in[i];out[i][IMAG] = 0;}//create a DFT plan and execute itfftwf_plan plan = fftwf_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);fftwf_execute(plan);//destroy the plan to prevent a memory leakfftwf_destroy_plan(plan);int hN = N >> 1;			// half of the length (N /2)int numRem = hN;		// the number of remaining elements// multiply the appropriate values by 2// (those that should be multiplied by 1 are left intact because they wouldn't change)for (int i = 1; i < hN; ++i) // 1,2,...,N/2 - 1 的项乘以2{out[i][REAL] *= 2;out[i][IMAG] *= 2;}// if the length is even, the number of remaining elements decreases by 1if (N % 2 == 0)numRem--;// if it's odd and greater than 1, the middle value must be multiplied by 2else if (N > 1)  // 奇数非空{out[hN][REAL] *= 2; out[hN][IMAG] *= 2;}// set the remaining values to 0// (multiplying by 0 gives 0, so we don't care about the multiplicands)memset(&out[hN + 1][REAL], 0, numRem * sizeof(fftwf_complex));// create an IDFT plan and execute itplan = fftwf_plan_dft_1d(N, out, out, FFTW_BACKWARD, FFTW_ESTIMATE); fftwf_execute(plan);// do some cleaningfftwf_destroy_plan(plan); fftwf_cleanup();// scale the IDFT output for (int i = 0; i < N; ++i){out[i][REAL] /= N;out[i][IMAG] /= N;}
}