离散哈特莱变换（DHT）用C语言实现

ゞ浴缸里的玫瑰 2022-02-16 00:37 224阅读 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][avatar 1]

DHT的正交证明：

![avatar][avatar 2]

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

![avatar][avatar 3]

其中:

![avatar][avatar 4]

DHT的优点

DHT为实值，避免了复数运算
    
    DHT正反变换形式基本一致
    
    DHT与DFT的转换容易实现

DHT的性质

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

1.线性性  
![avatar][avatar 5]

2.逆序列x(N-n)的DHT  
![avatar][avatar 6]

且  
![avatar][avatar 7]

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

3.循环位移的性质  
![avatar][avatar 8]

# 四、实现代码 #

【C语言】

#include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    
    const int N = 1024;
    const float PI = 3.1416;
    
    inline void swap(float &a, float &b)
    {
    	float t;
    	t = a;
    	a = b;
    	b = t;
    }
    
    void bitrp(float xreal[], float ximag[], int n)
    {
    	// 位反转置换 Bit-reversal Permutation
    	int i, j, a, b, p;
    
    	for (i = 1, p = 0; i < n; i *= 2)
    	{
    		p++;
    	}
    	for (i = 0; i < n; i++)
    	{
    		a = i;
    		b = 0;
    		for (j = 0; j < p; j++)
    		{
    			b = (b << 1) + (a & 1);    // b = b * 2 + a % 2;
    			a >>= 1;        // a = a / 2;
    		}
    		if (b > i)
    		{
    			swap(xreal[i], xreal[b]);
    			swap(ximag[i], ximag[b]);
    		}
    	}
    }
    
    void FFT(float xreal[], float ximag[], int n)
    {
    	// 快速傅立叶变换，将复数 x 变换后仍保存在 x 中，xreal, ximag 分		//	别是 x 的实部和虚部
    	float wreal[N / 2], wimag[N / 2], treal, timag, ureal, uimag, arg;
    	int m, k, j, t, index1, index2;
    
    	bitrp(xreal, ximag, n);
    
    	// 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * 		wimag [j] , j = 0, 1, ... , n / 2 - 1
    	arg = -2 * PI / n;
    	treal = cos(arg);
    	timag = sin(arg);
    	wreal[0] = 1.0;
    	wimag[0] = 0.0;
    	for (j = 1; j < n / 2; j++)
    	{
    		wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
    		wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
    	}
    
    	for (m = 2; m <= n; m *= 2)
    	{
    		for (k = 0; k < n; k += m)
    		{
    			for (j = 0; j < m / 2; j++)
    			{
    				index1 = k + j;
    				index2 = index1 + m / 2;
    				t = n * j / m;    // 旋转因子 w 的实部在 wreal [] 中					//的下标为 t
    				treal = wreal[t] * xreal[index2] - wimag[t] * 						ximag[index2];
    				timag = wreal[t] * ximag[index2] + wimag[t] * 						xreal[index2];
    				ureal = xreal[index1];
    				uimag = ximag[index1];
    				xreal[index1] = ureal + treal;
    				ximag[index1] = uimag + timag;
    				xreal[index2] = ureal - treal;
    				ximag[index2] = uimag - timag;
    			}
    		}
    	}
    }
    
    
    
    void FFT_test()
    {
    	float xreal[N] = {}, ximag[N] = {};
    
    	int n = 8;
    	int i = 0;
    	printf("请输入数据，格式(实部 虚部) : \n");
    	for (i = 0; i < 8; i++)
    	{
    		scanf("%f%f",xreal+i,ximag+i);
    	}
    
    	n = i;    // 要求 n 为 2 的整数幂
    	while (i > 1)
    	{
    		if (i % 2)
    		{
    			printf("%d is not a power of 2! ", n);
    		}
    		i /= 2;
    	}
    
    	FFT(xreal, ximag, n);
    	printf("=====================================\n");
    	printf("FFT:    i	      实部       虚部 \n");
    	for (i = 0; i < n; i++)
    	{
    		printf("     %4d       %8.4f    %8.4f ", i+1, xreal[i], 			ximag[i]);
    		printf("\n");
    	}
    	printf("===================================== \n");
    
    	
    	printf("DHT:    i	    结果 \n");
    	for (i = 0; i < n; i++)
    	{
    		printf("     %4d        %8.4f ", i+1, xreal[i]-ximag[i]);
    		printf("\n");
    	}
    	printf("=====================================\n ");
    }
    
    int main()
    {
    	FFT_test();
    	system("pause");
    	return 0;
    }

# END #

[avatar]: http://62.234.50.12/DHT/1.png
[avatar 1]: http://62.234.50.12/DHT/2.png
[avatar 2]: http://62.234.50.12/DHT/3.png
[avatar 3]: http://62.234.50.12/DHT/4.png
[avatar 4]: http://62.234.50.12/DHT/5.png
[avatar 5]: http://62.234.50.12/DHT/6.png
[avatar 6]: http://62.234.50.12/DHT/7.png
[avatar 7]: http://62.234.50.12/DHT/8.png
[avatar 8]: http://62.234.50.12/DHT/9.png