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

最新推荐文章于 2023-12-31 15:10:26 发布

盐厂的伙计

最新推荐文章于 2023-12-31 15:10:26 发布

阅读量6.3k

点赞数 4

分类专栏： matlab 文章标签：哈特莱变换 C语言实现 Matlab

本文链接：https://blog.csdn.net/Terrys0518/article/details/104114279

版权

matlab 专栏收录该内容

12 篇文章 1 订阅

订阅专栏

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

预备知识

1. DFT的意义

2. FFT实现

3. C语言编程

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

一般处理流程

1. Matlab 生成测试数据

Adc=2;  %直流分量幅度
A1=5;   %频率F1信号的幅度
A2=1.5;   %频率F2信号的幅度
F1=50;  %信号1频率(Hz)
F2=75;  %信号2频率(Hz)
Fs=256; %采样频率(Hz)
P1=-30; %信号1相位(度)
P2=90;  %信号相位(度)
N=256;  %采样点数
t=[0:1/Fs:N/Fs]; %采样时刻
%信号
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蝶形变换信号流图

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）

源代码如下，供大家参考

/*
 HFT验证
Auther:Terry
Date:  2020.1.30
*/
#include "stdio.h"
#include "math.h"

#define pi  3.1416
#define N   256

void daoxu(float * datas, int n);
void fht(float *Data,int Log2N);
void fftchange(float *f,float*p);

/*matlab生成的模拟采样数据  Fs = 256 N = 256*/
float mm[256] = {
6.4659,4.5027,1.1457,-1.8293,-0.7943,5.7235,7.8800,0.6097,-4.0687,0.9870,
6.1950,5.2453,1.9558,-1.2726,-1.6702,4.0584,8.3512,2.7280,-3.9023,-0.7336,
5.5067,5.8766,2.7658,-0.5706,-2.1237,2.3254,8.1404,4.7758,-3.0589,-2.2715,
4.4053,6.3156,3.5751,0.2061,-2.1793,0.7180,7.3058,6.5113,-1.6160,-3.4260,
2.9543,6.4730,4.3728,1.0107,-1.9020,-0.6085,5.9799,7.7358,0.2676,-4.0300,
1.2770,6.2685,5.1273,1.8208,-1.3777,-1.5535,4.3462,8.3204,2.3729,-3.9779,
-0.4539,5.6509,5.7824,2.6308,-0.6944,-2.0770,2.6104,8.2214,4.4500,-3.2442,
-2.0367,4.6156,6.2596,3.4405,0.0738,-2.1953,0.9700,7.4834,6.2529,-1.8923,
-3.2679,3.2159,6.4700,4.2416,0.8759,-1.9678,-0.4120,6.2276,7.5737,-0.0682,
-3.9730,1.5654,6.3304,5.0067,1.6857,-1.4784,-1.4249,4.6312,8.2703,2.0171,
-4.0343,-0.1701,5.7834,5.6832,2.4957,-0.8159,-2.0189,2.8981,8.2846,4.1166,
-3.4123,-1.7920,4.8157,6.1962,3.3057,-0.0577,-2.2016,1.2293,7.6465,5.9810,
-2.1553,-3.0951,3.4703,6.4572,4.1096,0.7412,-2.0263,-0.2050,6.4659,7.3940,
-0.3967,-3.8979,1.8514,6.3808,4.8836,1.5507,-1.5743,-1.2846,4.9124,8.2011,
1.6617,-4.0716,0.1169,5.9041,5.5796,2.3607,-0.9348,-1.9493,3.1877,8.3294,
3.7765,-3.5624,-1.5384,5.0051,6.1257,3.1708,-0.1880,-2.1979,1.4952,7.7943,
5.6964,-2.4042,-2.9083,3.7171,6.4350,3.9768,0.6068,-2.0771,0.0119,6.6937,
7.1974,-0.7167,-3.8051,2.1341,6.4201,4.7584,1.4157,-1.6650,-1.1326,5.1888,
8.1129,1.3079,-4.0896,0.4060,6.0129,5.4718,2.2257,-1.0508,-1.8681,3.4782,
8.3556,3.4309,-3.6944,-1.2768,5.1836,6.0487,3.0358,-0.3171,-2.1839,1.7670,
7.9263,5.4000,-2.6382,-2.7082,3.9555,6.4037,3.8433,0.4728,-2.1199,0.2384,
6.9102,6.9843,-1.0273,-3.6952,2.4127,6.4483,4.6313,1.2806,-1.7501,-0.9691,
5.4595,8.0057,0.9568,-4.0886,0.6964,6.1099,5.3602,2.0908,-1.1635,-1.7751,
3.7687,8.3629,3.0810,-3.8078,-1.0082,5.3509,5.9655,2.9008,-0.4447,-2.1592,
2.0440,8.0419,5.0928,-2.8567,-2.4956,4.1851,6.3638,3.7094,0.3391,-2.1540,
0.4740,7.1145,6.7554,-1.3274,-3.5686,2.6863
};

float fdata[N/2 + 1];
int main(void)
{
    int i;

    daoxu(mm,256);
    fht(mm,8);
    fftchange(fdata,mm);
    
    for(i = 0; i < N/2;i++)
    {
        if(i % 10 == 0)
            printf("\n");
        printf("%.3f  ",fdata[i]);
    }

    getchar();
    return 0;
}


/*=====================================================================
倒序计算  datas 为数据  n 为 长度 
================================================================*/
void daoxu(float * datas, int n)
{
	int nv2 = n / 2;
	int nm1 = n-1;
	int j = 0, i = 0, k;
	float t;
	do
	{
		if(i < j) 
		{
			t = datas[i];
			datas[i] = datas[j];
			datas[j] = t;
		}
		k = nv2;
		while(k <= j)
		{
			j = j - k;
			k = k / 2;
		}

		j = j + k;
		i = i + 1;
	}while(i < nm1);
}
/*哈特莱变换后数据处理，用来计算每个频率的幅度曲线*/
void fftchange(float *f,float*p)   
{
	int i;
	for( i = 0;i < N; i++)
		p[i] = p[i] / N ;
 	for( i = 0 ; i < (N/2) ; i++)
     {
        printf("\t%.3f,\t%.3f\n",p[i],p[N-i]);
		f[i] = p[i] * p[i] + p[N-i] *p[N-i];
        if(i > 0)
            f[i] = f[i] * 2 ;
        
        f[i] = sqrt(f[i]);  /*数据开方，可得到每个频率点对应的幅度*/
     }
}

/*===========================================================
哈特莱变换  data 为滤波后的数据  Log2N为阶数
=============================================================*/
void fht(float *Data,int Log2N)
{
    int length,i,j,k,step,i0,i1,i2,i3;
    float ck,sk,temp,temp0,temp1,temp2,temp3;

    length = 1<<Log2N;   
    for(i = 0; i < length; i += 2) 
    {
        temp = Data[i];
        Data[i] = temp + Data[i+1];
        Data[i+1] = temp - Data[i+1];
    }

    for(i = 2; i <= Log2N; i++)        
    {
        step = 1 << i;                       

        for(k = 0; k < length/step; k++)
        {
            i0=k * step; 
            i1 = k * step + step/2; 
            i2 = k * step + step/4;
            i3=k*step+step*3/4;

            temp0 = Data[i0] + Data[i1];
            temp1 = Data[i0] - Data[i1];
            temp2 = Data[i2] + Data[i3];
            temp3 = Data[i2] - Data[i3];
            Data[i0] = temp0;
            Data[i1] = temp1;
            Data[i2] = temp2;
            Data[i3] = temp3;
        }
        for(j = 1; j < step/4; j++)
        {
            ck = cos(2.0 * pi * j/step);
            sk = sin(2.0 * pi * j/step);
            for(k = 0;k < length/step;k++)
            {
                i0=k*step+j;i1=k*step+step/2+j;i2=k*step+step/2-j;i3=k*step+step-j;
                temp0 = Data[i0]+ck*Data[i1]+sk*Data[i3];
                temp1 = Data[i0]-ck*Data[i1]-sk*Data[i3];
                temp2 = Data[i2]+sk*Data[i1]-ck*Data[i3];
                temp3 = Data[i2]-sk*Data[i1]+ck*Data[i3];
                Data[i0] = temp0;
                Data[i1] = temp1;
                Data[i2] = temp2;
                Data[i3] = temp3;
            }
        }
    }
}

计算后的数据

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

Ayy =

  1 至 10 列

    2.1357    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  11 至 20 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  21 至 30 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  31 至 40 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  41 至 50 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  51 至 60 列

    5.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  61 至 70 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  71 至 80 列

    0.0000    0.0000    0.0000    0.0000    0.0000    1.5000    0.0000    0.0000    0.0000    0.0000

  81 至 90 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  91 至 100 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  101 至 110 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  111 至 120 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  121 至 130 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  131 至 140 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  141 至 150 列

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  151 至 156 列

    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计算结果一致，且准确反映了原信号的特征，计算正确。

盐厂的伙计

关注

4
点赞
踩
15

收藏

觉得还不错? 一键收藏
1
评论
快速哈特莱变换（FHT）介绍和C语言实现

过年整理资料的时候，发现了之前导师介绍的，一个号称专门针对离散实序列的变换，经分析总运算量为普通FFT的几乎一半(O(nLog(n))，而且完全没有复数。这么强的吗？之前也是一知半解，于是花了一个上午，重新验证了以下，顺便在这里把这个东西稍微普及一下，不知大家是否能用得上...预备知识1. DFT的意义2. FFT实现3. C语言编程原理部分可以参考西电的《数字信号处理》，或者...
复制链接

扫一扫

专栏目录