离散哈特莱变换和快速哈特莱变换的实现（MATLAB）

因为当初数字信号处理的课题设计不知道做什么题目，而且网络上也没有看到用MATLAB的实现代码，所以今天想起来想发到博客上，希望可以帮到大家，这是本人第一次发表博客，由于不会在博客上编辑公式，所以理论部分不发表了，只把实现代码发出来，希望大家捧场。

绘制频域采样信号的幅频特性图

function mstem(Xk)
%绘制频域采样序列向量Xk的幅频特性图（文件名mstem.m）
M = length(Xk);
k = 0:M-1;wk = 2*k/M; %产生M点DFT对应的采样点频率（关于pi归一化值）
stem(wk,abs(Xk),'.');box on; %绘制M点DFT的幅频特性图
xlabel('w/\pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(Xk))]);

这个函数是用来绘制信号的幅频响应。在下面两个文件中会用到这个函数。

DHT的实现和与DFT的比较

%离散哈特莱变换（文件名DHT.m）
N = 8;%采样点数
n = 0:N-1;
Fs = 2; %采样频率（Hz）
T = 1/Fs; %采样周期（s）
ft = cos(pi*n*T/4)+cos(pi*n*T/2); %采样信号，基频0.125Hz，最高频率0.25Hz
F1t = fft(ft,N);
F2t = zeros(1,N);
for i = n
    for j = 0:N-1 %哈特莱变换
        F2t(i+1) = F2t(i+1) + ft(j+1)*cos(2*pi/N*j*i) + ft(j+1)*sin(2*pi/N*j*i);
    end
end
ft = F2t;
for i = 0:N-1 %恢复傅里叶变换
    if(i == 0)
        F2t(1) = ft(1);
    else
        F2t(i+1) = ft(i+1)/2 + ft(N+1-i)/2 - 1j/2 * (ft(i+1) - ft(N+1-i));
    end
end
subplot(2,1,1);mstem(F1t);
title('(1a) 8点DFT[ ft ]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(F1t))]);
subplot(2,1,2);mstem(F2t);
title('(2a) 8点DHT[ ft ]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(F2t))]);

可以修改ft，n，Fs来实现不同的信号，采样点数量，采样频率。不过要符合采样定理。

FHT的实现和与FFT算法的比较

%快速哈特莱算法（文件名FHT.m）
N = 8; %8点采样
M = log2(N);
n = 0:N-1;
Fs = 2; %采样频率（Hz）
T = 1/Fs; %采样周期（s）
ft = cos(pi*n*T/4)+cos(pi*n*T/2); %采样信号，最高频率0.25Hz
FFt = fft(ft,N); %快速傅里叶变换
b = zeros(1,N); 
for i = 0:N-1 %奇偶顺序变换
    for k = 1:M
        b(i+1) = b(i+1)+mod(floor(i/2^(k-1)),2)*2^(M-k);
    end
end
Ft = ft;
for i = 0:N-1
    ft(i+1) = Ft(b(i+1)+1);
end
for i = 1:M %蝶形运算阶数
    for j = 0:2^(M-i)-1 %每一阶蝶形运算次数
        for k = 0:2^(i-1)-1 %k的范围
            if k == 0
                Ft(j*2^i+1) = ft(j*2^i+1) + ft(j*2^i+2^(i-1)+1);
                Ft(j*2^i+2^(i-1)+1) = ft(j*2^i+1) - ft(j*2^i+2^(i-1)+1);
            else
                Ft(k+j*2^i+1) = ft(k+j*2^i+1) + ft(k+j*2^i+2^(i-1)+1)*cos(2*pi/2^i*k) + ft(2^i*(j+1)+1-k)*sin(2*pi/2^i*k);
                Ft(k+j*2^i+2^(i-1)+1) = ft(k+j*2^i+1) - ft(k+j*2^i+2^(i-1)+1)*cos(2*pi/2^i*k) - ft(2^i*(j+1)+1-k)*sin(2*pi/2^i*k);
            end
        end
    end
    ft = Ft;
end
ft = Ft;
for i = n
    if(i == 0)
        Ft(1) = ft(1);
    else
        Ft(i+1) = ft(i+1)/2 + ft(N+1-i)/2 - 1j/2 * (ft(i+1) - ft(N+1-i));
    end
end
subplot(2,1,1);mstem(Ft);
title('(1a) 8点FHT[ft]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(Ft))]);
subplot(2,1,2);mstem(FFt);
title('(2a) 8点FFT[ft]');xlabel('ω/π');ylabel('幅度');
axis([0,2,0,1.2*max(abs(FFt))]);

重点，和DHT一样，可以修改ft，n，Fs来实现不同的信号，采样点数量，采样频率。不过要符合采样定理。
欢迎指导。