离散哈特莱变换(DHT)用C语言实现

ゞ 浴缸里的玫瑰 2022-02-16 00:37 321阅读 0赞

离散哈特莱变换(DHT)

摘 要

离散哈特莱变换(DHT)是一种与傅里叶变换相关的转换。类似于离散傅里叶变换。与傅里叶变换在信号处理及其他相关领域有相似的应用。
本设计介绍了DHT的定义以及使用C语言实现其算法。

关键字:傅里叶变换 哈特莱变换

二、设计平台

Linux平台、GCC编译器、VIM、windows、Visual Studio 2017

三、设计原理

DHT的定义
设x(n),n=0,1,…,N−1x(n),n=0,1,…,N−1,为一实序列,其DHT定义为:

avatar

式中cas(a)=cos(a)+sin(a)cas(a)=cos(a)+sin(a)逆变换(IDHT)为:

avatar

DHT的正交证明:

avatar

DHT和DFT关系
用X(k)表示实序列x(n)的,用XH(k)表示x(n)的DHT,分别用XHe(k), XHo(k)表示XH(k)的偶对称分量与奇对称分量,即:

avatar

其中:

avatar

DHT的优点

  1. DHT为实值,避免了复数运算
  2. DHT正反变换形式基本一致
  3. DHTDFT的转换容易实现

DHT的性质

DHT的性质与DFT的性质类似,但由于DHT是实序列间的变换,有些性质有具体的表达形式。这里只给出结论。
设x(n)、y(n)的DHT分别为Xh(k)、Yh(k)。用符号x(n)↔Xh(k)表示Xh(k)=DHT(x(n))。

1.线性性
avatar

2.逆序列x(N-n)的DHT
avatar


avatar

当k=0时,可得Xh(N)=Xh(0)。

3.循环位移的性质
avatar

四、实现代码

【C语言】

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. const int N = 1024;
  5. const float PI = 3.1416;
  6. inline void swap(float &a, float &b)
  7. {
  8. float t;
  9. t = a;
  10. a = b;
  11. b = t;
  12. }
  13. void bitrp(float xreal[], float ximag[], int n)
  14. {
  15. // 位反转置换 Bit-reversal Permutation
  16. int i, j, a, b, p;
  17. for (i = 1, p = 0; i < n; i *= 2)
  18. {
  19. p++;
  20. }
  21. for (i = 0; i < n; i++)
  22. {
  23. a = i;
  24. b = 0;
  25. for (j = 0; j < p; j++)
  26. {
  27. b = (b << 1) + (a & 1); // b = b * 2 + a % 2;
  28. a >>= 1; // a = a / 2;
  29. }
  30. if (b > i)
  31. {
  32. swap(xreal[i], xreal[b]);
  33. swap(ximag[i], ximag[b]);
  34. }
  35. }
  36. }
  37. void FFT(float xreal[], float ximag[], int n)
  38. {
  39. // 快速傅立叶变换,将复数 x 变换后仍保存在 x 中,xreal, ximag 分 // 别是 x 的实部和虚部
  40. float wreal[N / 2], wimag[N / 2], treal, timag, ureal, uimag, arg;
  41. int m, k, j, t, index1, index2;
  42. bitrp(xreal, ximag, n);
  43. // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
  44. arg = -2 * PI / n;
  45. treal = cos(arg);
  46. timag = sin(arg);
  47. wreal[0] = 1.0;
  48. wimag[0] = 0.0;
  49. for (j = 1; j < n / 2; j++)
  50. {
  51. wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
  52. wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
  53. }
  54. for (m = 2; m <= n; m *= 2)
  55. {
  56. for (k = 0; k < n; k += m)
  57. {
  58. for (j = 0; j < m / 2; j++)
  59. {
  60. index1 = k + j;
  61. index2 = index1 + m / 2;
  62. t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中 //的下标为 t
  63. treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
  64. timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
  65. ureal = xreal[index1];
  66. uimag = ximag[index1];
  67. xreal[index1] = ureal + treal;
  68. ximag[index1] = uimag + timag;
  69. xreal[index2] = ureal - treal;
  70. ximag[index2] = uimag - timag;
  71. }
  72. }
  73. }
  74. }
  75. void FFT_test()
  76. {
  77. float xreal[N] = {}, ximag[N] = {};
  78. int n = 8;
  79. int i = 0;
  80. printf("请输入数据,格式(实部 虚部) : \n");
  81. for (i = 0; i < 8; i++)
  82. {
  83. scanf("%f%f",xreal+i,ximag+i);
  84. }
  85. n = i; // 要求 n 为 2 的整数幂
  86. while (i > 1)
  87. {
  88. if (i % 2)
  89. {
  90. printf("%d is not a power of 2! ", n);
  91. }
  92. i /= 2;
  93. }
  94. FFT(xreal, ximag, n);
  95. printf("=====================================\n");
  96. printf("FFT: i 实部 虚部 \n");
  97. for (i = 0; i < n; i++)
  98. {
  99. printf(" %4d %8.4f %8.4f ", i+1, xreal[i], ximag[i]);
  100. printf("\n");
  101. }
  102. printf("===================================== \n");
  103. printf("DHT: i 结果 \n");
  104. for (i = 0; i < n; i++)
  105. {
  106. printf(" %4d %8.4f ", i+1, xreal[i]-ximag[i]);
  107. printf("\n");
  108. }
  109. printf("=====================================\n ");
  110. }
  111. int main()
  112. {
  113. FFT_test();
  114. system("pause");
  115. return 0;
  116. }

END

发表评论

表情:
评论列表 (有 0 条评论,321人围观)

还没有评论,来说两句吧...

相关阅读

    相关 兄弟的家庭教育

    莱特兄弟的三个设计经验 你可以从他们如何使得一辆自行车飞行中学到很多东西 克莱夫·汤普森 大家都知道莱特兄弟是第一个实现动力飞行的人——他们的飞机于1903年12月17日在

    相关 C语言实现radon变换

    因为实验室的要求,需要用C语言实现radon变换,对于刚刚我这个大二刚接触这种纯理论纯算法的人来说,真心不是一般的纠结。。。网络上搜索了好久,也没有找到。最后在CSDN论坛见过