devbox
IT小白
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

快速哈特莱变换（FHT）介绍和C语言实现

作者：IT小白 | 2024-07-13 05:37:28

赞

踩

fht

过年整理资料的时候，发现了之前导师介绍的，一个号称专门针对离散实序列的变换，经分析总运算量为普通FFT的几乎一半(O(nLog(n))，而且完全没有复数。这么强的吗？之前也是一知半解，于是花了一个上午，重新验证了以下，顺便在这里把这个东西稍微普及一下，不知大家是否能用得上...

预备知识

1. DFT的意义

2. FFT实现

3. C语言编程

原理部分可以参考西电的《数字信号处理》，或者参考https://www.cnblogs.com/RRRR-wys/p/10090007.html

一般处理流程

1. Matlab 生成测试数据

  1. Adc=2;  %直流分量幅度
  2. A1=5;   %频率F1信号的幅度
  3. A2=1.5;   %频率F2信号的幅度
  4. F1=50;  %信号1频率(Hz)
  5. F2=75;  %信号2频率(Hz)
  6. Fs=256; %采样频率(Hz)
  7. P1=-30; %信号1相位(度)
  8. P2=90;  %信号相位(度)
  9. N=256;  %采样点数
  10. t=[0:1/Fs:N/Fs]; %采样时刻
  11. %信号
  12. S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180) + 0.2 * rand();

2. 加窗滤波

  根据采样定理，假设采样频率Fs，则要求原采样信号的最高频率小于Fs/2，否则傅里叶变换后的信号会出现频谱叠加的现象，实际操作时，由于采样精度和噪声的因素，采样得到的离散信号还会包含了高频信号（大于Fs/2），所以需要用加窗滤波的方式去除，否则计算后的数据不准确，这里略过。

2. 倒序计算

原理参考《数字信号处理》，主要是体现为以下规律

 3. DIT-FHT蝶形变换信号流图

256点的类似

 4. 运算后数据处理

    N点数据处理后，结果还是N点（N=2^n），需要进行一下处理可得到原信号的频率和幅度信息
    x(i) = x(i) / N;    i∈（0 ~N-1）
    y(i) = x(i) ^2 + x(N - i)^2;     i∈（0 ~N/2）
    y(i) = y(i) * 2;     i∈（1 ~N/2）
    y(i) = sqrt(y(i));   i∈（0 ~N/2）

   即得到 0Hz -  Fs / 2 对应的幅度
   每个点对应频率  f(i) = i * Fs / N ; （i ∈ 0 - N/2）
   对应的幅度        y(i) = y(i);           （i ∈ 0 - N/2）

源代码如下，供大家参考 

  1. /*
  2.  HFT验证
  3. Auther:Terry
  4. Date:  2020.1.30
  5. */
  6. #include "stdio.h"
  7. #include "math.h"
  8.  
  9. #define pi  3.1416
  10. #define N   256
  11.  
  12. void daoxu(float * datas, int n);
  13. void fht(float *Data,int Log2N);
  14. void fftchange(float *f,float*p);
  15.  
  16. /*matlab生成的模拟采样数据  Fs = 256 N = 256*/
  17. float mm[256] = {
  18. 6.4659,4.5027,1.1457,-1.8293,-0.7943,5.7235,7.8800,0.6097,-4.0687,0.9870,
  19. 6.1950,5.2453,1.9558,-1.2726,-1.6702,4.0584,8.3512,2.7280,-3.9023,-0.7336,
  20. 5.5067,5.8766,2.7658,-0.5706,-2.1237,2.3254,8.1404,4.7758,-3.0589,-2.2715,
  21. 4.4053,6.3156,3.5751,0.2061,-2.1793,0.7180,7.3058,6.5113,-1.6160,-3.4260,
  22. 2.9543,6.4730,4.3728,1.0107,-1.9020,-0.6085,5.9799,7.7358,0.2676,-4.0300,
  23. 1.2770,6.2685,5.1273,1.8208,-1.3777,-1.5535,4.3462,8.3204,2.3729,-3.9779,
  24. -0.4539,5.6509,5.7824,2.6308,-0.6944,-2.0770,2.6104,8.2214,4.4500,-3.2442,
  25. -2.0367,4.6156,6.2596,3.4405,0.0738,-2.1953,0.9700,7.4834,6.2529,-1.8923,
  26. -3.2679,3.2159,6.4700,4.2416,0.8759,-1.9678,-0.4120,6.2276,7.5737,-0.0682,
  27. -3.9730,1.5654,6.3304,5.0067,1.6857,-1.4784,-1.4249,4.6312,8.2703,2.0171,
  28. -4.0343,-0.1701,5.7834,5.6832,2.4957,-0.8159,-2.0189,2.8981,8.2846,4.1166,
  29. -3.4123,-1.7920,4.8157,6.1962,3.3057,-0.0577,-2.2016,1.2293,7.6465,5.9810,
  30. -2.1553,-3.0951,3.4703,6.4572,4.1096,0.7412,-2.0263,-0.2050,6.4659,7.3940,
  31. -0.3967,-3.8979,1.8514,6.3808,4.8836,1.5507,-1.5743,-1.2846,4.9124,8.2011,
  32. 1.6617,-4.0716,0.1169,5.9041,5.5796,2.3607,-0.9348,-1.9493,3.1877,8.3294,
  33. 3.7765,-3.5624,-1.5384,5.0051,6.1257,3.1708,-0.1880,-2.1979,1.4952,7.7943,
  34. 5.6964,-2.4042,-2.9083,3.7171,6.4350,3.9768,0.6068,-2.0771,0.0119,6.6937,
  35. 7.1974,-0.7167,-3.8051,2.1341,6.4201,4.7584,1.4157,-1.6650,-1.1326,5.1888,
  36. 8.1129,1.3079,-4.0896,0.4060,6.0129,5.4718,2.2257,-1.0508,-1.8681,3.4782,
  37. 8.3556,3.4309,-3.6944,-1.2768,5.1836,6.0487,3.0358,-0.3171,-2.1839,1.7670,
  38. 7.9263,5.4000,-2.6382,-2.7082,3.9555,6.4037,3.8433,0.4728,-2.1199,0.2384,
  39. 6.9102,6.9843,-1.0273,-3.6952,2.4127,6.4483,4.6313,1.2806,-1.7501,-0.9691,
  40. 5.4595,8.0057,0.9568,-4.0886,0.6964,6.1099,5.3602,2.0908,-1.1635,-1.7751,
  41. 3.7687,8.3629,3.0810,-3.8078,-1.0082,5.3509,5.9655,2.9008,-0.4447,-2.1592,
  42. 2.0440,8.0419,5.0928,-2.8567,-2.4956,4.1851,6.3638,3.7094,0.3391,-2.1540,
  43. 0.4740,7.1145,6.7554,-1.3274,-3.5686,2.6863
  44. };
  45.  
  46. float fdata[N/2 + 1];
  47. int main(void)
  48. {
  49.     int i;
  50.  
  51.     daoxu(mm,256);
  52.     fht(mm,8);
  53.     fftchange(fdata,mm);
  54.     
  55.     for(i = 0; i < N/2;i++)
  56.     {
  57.         if(i % 10 == 0)
  58.             printf("\n");
  59.         printf("%.3f  ",fdata[i]);
  60.     }
  61.  
  62.     getchar();
  63.     return 0;
  64. }
  65.  
  66.  
  67. /*=====================================================================
  68. 倒序计算  datas 为数据  n 为 长度 
  69. ================================================================*/
  70. void daoxu(float * datas, int n)
  71. {
  72. 	int nv2 = n / 2;
  73. 	int nm1 = n-1;
  74. 	int j = 0, i = 0, k;
  75. 	float t;
  76. 	do
  77. 	{
  78. 		if(i < j) 
  79. 		{
  80. 			t = datas[i];
  81. 			datas[i] = datas[j];
  82. 			datas[j] = t;
  83. 		}
  84. 		k = nv2;
  85. 		while(k <= j)
  86. 		{
  87. 			j = j - k;
  88. 			k = k / 2;
  89. 		}
  90.  
  91. 		j = j + k;
  92. 		i = i + 1;
  93. 	}while(i < nm1);
  94. }
  95. /*哈特莱变换后数据处理，用来计算每个频率的幅度曲线*/
  96. void fftchange(float *f,float*p)   
  97. {
  98. 	int i;
  99. 	for( i = 0;i < N; i++)
  100. 		p[i] = p[i] / N ;
  101.  	for( i = 0 ; i < (N/2) ; i++)
  102.      {
  103.         printf("\t%.3f,\t%.3f\n",p[i],p[N-i]);
  104. 		f[i] = p[i] * p[i] + p[N-i] *p[N-i];
  105.         if(i > 0)
  106.             f[i] = f[i] * 2 ;
  107.         
  108.         f[i] = sqrt(f[i]);  /*数据开方，可得到每个频率点对应的幅度*/
  109.      }
  110. }
  111.  
  112. /*===========================================================
  113. 哈特莱变换  data 为滤波后的数据  Log2N为阶数
  114. =============================================================*/
  115. void fht(float *Data,int Log2N)
  116. {
  117.     int length,i,j,k,step,i0,i1,i2,i3;
  118.     float ck,sk,temp,temp0,temp1,temp2,temp3;
  119.  
  120.     length = 1<<Log2N;   
  121.     for(i = 0; i < length; i += 2) 
  122.     {
  123.         temp = Data[i];
  124.         Data[i] = temp + Data[i+1];
  125.         Data[i+1] = temp - Data[i+1];
  126.     }
  127.  
  128.     for(i = 2; i <= Log2N; i++)        
  129.     {
  130.         step = 1 << i;                       
  131.  
  132.         for(k = 0; k < length/step; k++)
  133.         {
  134.             i0=k * step; 
  135.             i1 = k * step + step/2; 
  136.             i2 = k * step + step/4;
  137.             i3=k*step+step*3/4;
  138.  
  139.             temp0 = Data[i0] + Data[i1];
  140.             temp1 = Data[i0] - Data[i1];
  141.             temp2 = Data[i2] + Data[i3];
  142.             temp3 = Data[i2] - Data[i3];
  143.             Data[i0] = temp0;
  144.             Data[i1] = temp1;
  145.             Data[i2] = temp2;
  146.             Data[i3] = temp3;
  147.         }
  148.         for(j = 1; j < step/4; j++)
  149.         {
  150.             ck = cos(2.0 * pi * j/step);
  151.             sk = sin(2.0 * pi * j/step);
  152.             for(k = 0;k < length/step;k++)
  153.             {
  154.                 i0=k*step+j;i1=k*step+step/2+j;i2=k*step+step/2-j;i3=k*step+step-j;
  155.                 temp0 = Data[i0]+ck*Data[i1]+sk*Data[i3];
  156.                 temp1 = Data[i0]-ck*Data[i1]-sk*Data[i3];
  157.                 temp2 = Data[i2]+sk*Data[i1]-ck*Data[i3];
  158.                 temp3 = Data[i2]-sk*Data[i1]+ck*Data[i3];
  159.                 Data[i0] = temp0;
  160.                 Data[i1] = temp1;
  161.                 Data[i2] = temp2;
  162.                 Data[i3] = temp3;
  163.             }
  164.         }
  165.     }
  166. }
  167.

 计算后的数据

Matlab FFT计算结果来检查HFT的计算数据

  1. Ayy =
  2.  
  3.   1 至 10 列
  4.  
  5.     2.1357    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  6.  
  7.   11 至 20 列
  8.  
  9.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  10.  
  11.   21 至 30 列
  12.  
  13.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  14.  
  15.   31 至 40 列
  16.  
  17.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  18.  
  19.   41 至 50 列
  20.  
  21.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  22.  
  23.   51 至 60 列
  24.  
  25.     5.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  26.  
  27.   61 至 70 列
  28.  
  29.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  30.  
  31.   71 至 80 列
  32.  
  33.     0.0000    0.0000    0.0000    0.0000    0.0000    1.5000    0.0000    0.0000    0.0000    0.0000
  34.  
  35.   81 至 90 列
  36.  
  37.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  38.  
  39.   91 至 100 列
  40.  
  41.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  42.  
  43.   101 至 110 列
  44.  
  45.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  46.  
  47.   111 至 120 列
  48.  
  49.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  50.  
  51.   121 至 130 列
  52.  
  53.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  54.  
  55.   131 至 140 列
  56.  
  57.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  58.  
  59.   141 至 150 列
  60.  
  61.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
  62.  
  63.   151 至 156 列
  64.  
  65.     0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

 C语言计算数据

 绘制曲线

 

原信号：
      Adc=2;  %直流分量幅度 加 rand(0,0.2)的随机噪声
      A1=5;   %频率F1信号的幅度
      A2=1.5;   %频率F2信号的幅度
HFT变换得到的数据
     直流分量    2.136，且近似等于Matlab的计算值2.1357
     F1 的幅度  为 5，等于 A1 和 FFT计算数据
     F2 的幅度 为 1.5，等于 A2 和 FFT计算数据

由此可见，HTF计算结果和FFT计算结果一致，且准确反映了原信号的特征，计算正确。

本文内容由网友自发贡献，转载请注明出处：https://www.wpsshop.cn/w/IT小白/article/detail/817829
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号