2019-12-24

离散哈特莱变换（DHT）

author: RoWe98(RexRowe)

离散哈特莱变换（DHT）

摘要

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

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

二、设计平台

Linux平台
GCC编译器
vim
windows
Visual Studio 2017

三、设计原理

DHT的定义
设$\mathrm{x}(\mathrm{n}), \mathrm{n}=0,1, \ldots, \mathrm{N}-1 \mathrm{x}(\mathrm{n}), \mathrm{n}=0,1, \ldots, \mathrm{N}-1$为一实序列，其DHT定义为：

$$
\begin{aligned} X_{H}(k) &=\operatorname{DHTEx(n)} ] \ &=\sum_{n=0}^{N-1} x(n) \operatorname{cas}\left(\frac{2 \pi}{N} k n\right), \quad k=0,1, \cdots, N-1 \end{aligned}
$$

式中$\operatorname{cas}(a)=\cos (a)+\sin (a) \operatorname{cas}(a)=\cos (a)+\sin (a)$逆变换$(\mathrm{IDHT})$为：

$$
\begin{aligned} x(n) &=\operatorname{IDHT}\left[X_{H}(k)\right] \ &=(1 / N) \sum_{k=0}^{N-1} X_{H}(k) \operatorname{cas}\left(\frac{2 \pi}{N} k n\right), \quad n=0,1, \cdots, N-1 \end{aligned}
$$

DHT的正交证明：

$$
\sum_{k=0}^{N-1} \operatorname{cas}\left(\frac{2 \pi}{N} k n\right) \operatorname{cas}\left(\frac{2 \pi}{N} k m\right)=\sum_{k=0}^{N-1}\left[\cos \left(\frac{2 \pi}{N} k(n-m)\right)+\sin \left(\frac{2 \pi}{N} k(n+m)\right)\right]=\left{\begin{array}{ll}{N,} & {k=0} \ {0,} & {k \neq 0}\end{array}\right.
$$

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

$X_{\mathrm{H}}(k)=X_{\mathrm{He}}(k)+X_{\mathrm{Ho}}(k)$

其中:

$$
\begin{array}{c}{X_{H e}(k)=\frac{1}{2}\left[X_{H}(k)+X_{H}(N-k)\right]}{X_{H o}(k)=\frac{1}{2}\left[X_{H}(k)-X_{H}(N-k)\right]} \ {X(k)=X_{H e}(k)-j X_{H o}(k)} \ {X_{H}(k)=\operatorname{Re}[X(k)]-\operatorname{Im}[X(k)]} \ {X(k)=\frac{1}{2}\left[X_{H}(k)+X_{H}(N-k)\right]-\frac{1}{2} j\left[X_{H}(k)-X_{H}(N-k)\right]}\end{array}
$$

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.线性性

2.逆序列x(N-n)的DHT
$x(N-n) \leftrightarrow X_{\mathrm{H}}(\mathrm{N}-k)$

且
$$
\mathrm{X}{\mathrm{H}}(\mathrm{N}-\mathrm{k})=\sum{n=0}^{\mathrm{N}-1} x(n)\left[\cos \left(\frac{2 \pi}{\mathrm{N}} k n\right)-\sin \left(\frac{2 \pi}{N} k n\right)\right], \quad k=0,1, \cdots, N-1
$$

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

3.循环位移的性质
$$
\begin{array}{l}{x\left(\left(n-n_{0}\right){N} R{N}(n) \leftrightarrow X_{\mathrm{H}}(k) \cos \left(\frac{2 \pi}{N} k n_{0}\right)+X_{\mathrm{H}}(N-k) \sin \left(\frac{2 \pi}{N} k n_{0}\right)\right)} \ {x\left(\left(n+n_{0}\right){N} R{N}(n) \leftrightarrow X_{\mathrm{H}}(k) \cos \left(\frac{2 \pi}{N} k n_{0}\right)-X_{\mathrm{H}}(N-k) \sin \left(\frac{2 \pi}{N} k n_{0}\right)\right)}\end{array}
$$

四、实现代码

【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;
}

运行效果

请输入数据，格式(实部 虚部) : 
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
=====================================
FFT:    i          实部       虚部 
        1        36.0000      0.0000 
        2        -4.0000      9.6568 
        3        -4.0000      4.0000 
        4        -4.0000      1.6569 
        5        -4.0000      0.0000 
        6        -4.0000     -1.6568 
        7        -4.0000     -4.0000 
        8        -4.0001     -9.6569 
===================================== 
DHT:    i        结果 
        1         36.0000 
        2        -13.6568 
        3         -8.0000 
        4         -5.6568 
        5         -4.0000 
        6         -2.3432 
        7         -0.0000 
        8          5.6568 
=====================================

Rex's Blog

Code and Passion

离散哈特莱变换（DHT）

离散哈特莱变换（DHT）

摘要

二、设计平台

三、设计原理

四、实现代码

运行效果

End

离散哈特莱变换（DHT）

摘 要

二、设计平台

三、设计原理

四、实现代码

运行效果

End

摘要