Hilbert Transform in ARM M4

DSP RTOS
Viewed: - times 05 Apr 2017 纠错

Hilbert Transform in ARM M4

这里我们主要讨论一下M4中如何实现 Hilbert Transform。
当然不光是M4中，其他同等级MCU一个道理。

理论

根据Hilbert Transform 定义可以得知频率响应可以由傅立叶变换给出：

$ H(\omega)=F\{h\}(\omega)=-i \bullet sgn(\omega) $

而 $sng(\omega) = \begin{cases} 1, &\quad for \;\omega > 0, \\ 0, &\quad for \;\omega = 0, \\ -1, &\quad for \;\omega < 0. \\ \end{cases}$

简单来说，对于一组信号$x(t)$ 假设长度为1024,经过1024点FFT后，
把结果 $x\_fft(f)$ 分成两组，
前半和后半 $a+bi,c+di$ ,

index从1开始，1和513的位置的值，乘$0$
在2～512的 $(a+bi)\bullet -i = b-ai$
在514～1024的 $(c+di) \bullet i = -d +ci$

到这的结果，进行IFFT，
得到的就是 $x\_h(t)$, 也就是matlab中hilbert()结果的虚部，
想获取这组信号的hilbert magnitude值还要进行abs()
所以matlab中所得到magnitude结果为 $abs(x(t)+x\_h(t)i)$
以上就是按照理论的过程。

但在matlab中 open('hilbert.m')打开hilbert函数
你会发现，并不是按照刚刚介绍的步骤来的，
而是，用了一种更为简单的方法。可以参考matlab hilbert
大概意思就如下：

$x\_fft(f) \bullet h(i)$ 得到的结果，进行IFFT就可以直接得到我们所需的magnitude值。

而 $h(i) = \begin{cases} 1, &\quad for \; i=1,(1024/2+1), \\ 2, &\quad for \; i=2,3,...,(1024/2), \\ 0, &\quad for \; i=(1024/2)+2,...,1024. \\ \end{cases} \quad \text{i is index from 1~1024}$

实践

按照后面介绍的matlab方法，就有了一下在M4上的hilbert实现。
有想参考的朋友要注意一下point的size
代码中有使用到arm_math.h

void hilbertTF(ValueFloat_t *z, uint16_t n, ValueFloat_t *output){
    status = arm_cfft_radix4_init_f32( &cfftinst_,n,0,1);      //forward fft
    arm_cfft_radix4_f32(&cfftinst_,z);      //z: input & output
    arm_fill_f32(0.0,checkpoint,n);

    arm_fill_f32(2.0,checkpoint,n/2+1);
    checkpoint[0]=1;
    checkpoint[513]=1;

    arm_cmplx_mult_real_f32(z,checkpoint,output,n);

    status = arm_cfft_radix4_init_f32( &cfftinst_,n,1,1);      //inverse fft
    arm_cfft_radix4_f32(&cfftinst_,output);      //z: input & output

    return;
}

结果

最后附上结果验证

先是fft结果比较：

突然误差为 $10^{-5}$，还算是在可接受范围之内
接下来是最终结果比较：

Reference:

matlab hilbert

TOP

Next: 博客补全 Previous: UDOO x86 unboxing and ubuntu installation