使用Matlab进行信号和图像压缩的Walsh Hadamard变换
Walsh-Hadamard变换是一种非正弦波的正交变换，广泛用于信号和图像处理。在这种变换中，信号被分解成一组基函数（类似于傅里叶的谐波）。
这些基础函数是沃尔什函数，即带有+1或-1的矩形或方波。此外，像快速傅里叶变换（FFT）一样，沃尔什-哈达马德变换有一个快速版本，称为快速沃尔什-哈达马德变换（FWHT）。
Walsh-Hamard在科学和工程领域有广泛的应用。这些应用包括图像处理、语音处理、信号和图像压缩、功率谱分析、扩频分析等。
本教程将使用Matlab研究信号和图像压缩。我们将研究其理论和一维和二维信号的计算。
前提条件
要跟上本教程，你需要。
- 安装[MATLAB]。
- 对[MATLAB]基础知识的正确理解。
正向和逆向WHT
正向（FWHT）和反向（IFWHT）是对称的，两者使用相同的计算过程。N 我们对一个长度为x(t) 的信号的FWHT和IWHT定义为。
FWHT: y_n=\frac{1}{N}\sum x_i=\sum y_n WAL(n,i)$ for i = 1,2,...,N
其中。
WAL(n, i) 是沃尔什函数。
要计算WHT，信号N 的长度必须是2(i.e., N=2^2) 的幂。如果长度不是2的幂，就用零来使其达到2的幂。因此，你通过矩阵乘法来计算FWHT的方程式。
x(t) ={x_1, x_2, ..., x_n}
y = \frac{1}{n}(H_n, x)
为了详细说明这一点，我们有。
\begin{vmatrix}y_1\\y_2\\:\\:\\y_N\end{vmatrix}=\frac{1}{N}\begin{vmatrix}H_{11}&H_{12}....&H_{1N}\\H_{21}&H_{22}....&H_{2N}\\:\\:\\H_{N1}&H_{N2}....&H_{NN}\end{vmatrix}\begin{vmatrix}x_1\\x_2\\:\\:\\x_N\end{vmatrix}
这里，H_N 是Hadamard矩阵，其中N =2^n 。同样地，反变换可以通过以下方式得到。
x = H_n^{-1}=H_ny
H_n^{-1}=H_n^T=H_n
这是因为H_n 是正交的和内卷的。正交意味着有一个垂直的矢量。因此，当我们算出这一点时，我们将得到。
\begin{vmatrix}x_1\\x_2\\:\\:\\x_N\end{vmatrix}=\frac{1}{N}\begin{vmatrix}H_{11}&H_{12}....&H_{1N}\\H_{21}&H_{22}....&H_{2N}\\:\\:\\H_{N1}&H_{N2}....&H_{NN}\end{vmatrix}\begin{vmatrix}y_1\\y_2\\:\\:\\y_N\end{vmatrix}
哈达玛德矩阵
它是2x2,4x4,8x8, e.t.c. 的方形矩阵。尺寸由2^N 给出。一阶Hadamard矩阵(H_1)由以下公式给出。
H_1=
\left(\begin{array}{cc}
1 & 1\\
1 & -1
\end{array}\right)
第二阶Hadamard矩阵由以下公式给出。
H_2=
\left(\begin{array}{cc}
1 & 1 & 1 & 1\\
1 & -1 & 1 & -1\\
1 & 1 & -1 & -1\\
1 & -1 & -1 & 1
\end{array}\right)
矩阵H_2 是由H_1 使用克朗克乘积得到的。
H_2=H_1 * H_1
其中
* 是克朗克乘积。克朗克积的意思是，如果我们有两个矩阵A和B，那么它就是矩阵A和B元素的乘积，如下图所示。
A*B=
\left(\begin{array}{cc}
A(1,1).B & A(1,2).B & ... & A(1,N).B\\
: & : & & : &\\
: & : & & :\\
A(N,1).B & A(N,2).B & ... & A(N,N).B
\end{array}\right)
高阶哈达玛德矩阵H_N ，可以通过下面的迭代规则得到。
H_N=H_1 * H_{N-1}
我们想看看如何从Hadamard得到Walsh矩阵。
考虑一下哈达玛德矩阵（H_2）。
如果我们看第一行，我们会发现元素的符号从第一列到最后一列都没有变化。然而，符号从1变为-1，然后变为1，再变为-1，这在第二行中加起来有三个变化，然后你再转到第三行和第四行。在这里，我们考虑符号如何变化。如果我们按照符号变化的数量递增来排列各行，就变成了沃尔什矩阵，如下图所示。
图像的Walsh-Hadamard变换
对于二维信号（N×N的图像），你可以通过以下矩阵方程得到正向沃尔什法变换。
Y=\frac{1}{N^2}(H_N*H_N)
其中
x是输入图像。
y是输出变换矩阵。
H_N 是Hadamard或Walsh矩阵。
请注意，如果你的输入图像的尺寸是MxN，那么你将得到\frac{1}{N^2} ，而不是\frac{1}{MN^2} 。这是因为M和N必须是2的整数倍。
获得反变换的方法与正向变换的操作类似。
X=H_NYH_N
Matlab中一维信号的WHT实例
let $$x(n)=
\left(\begin{array}{cc}
1 & 5 & -3 & 2\end{array}\right)
使用N（哈达玛德矩阵）为4，我们可以找到矩阵的反变换和正变换。这可以通过下面的代码来完成。
x = [1;5;-3;2]; % signal
H = hadamard(4); %hadamard matrix
y = (1/4)*H*x; %forward transform
x_r = H*y; %Inverse transform formula
在上面的代码中，我们使用了寻找矩阵的正向和反向变换的公式。我们用1/4*H*x ，因为它是一个一维信号，而对于反变换，公式是H*y 。因此，h 是Hadamard矩阵，而x 是我们的矩阵。
另外，matlab有内置的函数fwht ，用于正向变换，ifwht ，用于反向变换。为了在我们的例子中应用这一点，我们将有下面的代码。
Yy = fwht(x,4, 'hadamard'); %Forward WHT
xr = ifwht(Yy,4, 'hadamard'); %Inverse WHT
使用Matlab的二维信号的WHT实例
let $$x(n)=
\left(\begin{array}{cc}
250 & 128 & 100 & 25\\
40 & 102 & 95 & 48\\
75 & 64 & 120 & 97\\
10 & 255 & 0 & 55\end{array}\right)
就像我们对一维信号所做的那样，我们也将在这里做这件事。我们找到我们矩阵的正变换和反变换。Hadamard矩阵仍然是4。对于这个二维信号，公式将是1/16*H*x 。
I = [250, 128, 100, 25;40 102 95 48;75, 64, 120, 97;10, 225, 0, 55]; % signal
H = hadamard(4); %hadamard matrix
y = (1/16)*H*I*H; %forward transform
x_r = H*y*H; %Inverse transform formula
快速WHT
这个方法的复杂度为O(NlogN)。它只使用加法和减法。为了解其工作原理，我们使用了下图所示的蝴蝶结构。
蓝色箭头表示加法，而红色箭头表示减法。它使计算更快。Matlab的内置函数fwht 和ifwht 是基于这种快速的WHT算法的。
Matlab使用FWHT实现信号过滤的代码
这个程序首先在一个.mat文件中加载心电图信号。为了加载数据，我们使用load 命令。
为使其发挥作用，确保该文件在当前目录下。
%program for 1D hadamard transform
load ecg; %loading ECG signal
然后我们提取信号的长度，对应于2^n$。由于我们的信号是5000个样本，而5000不是2的幂，我们取4096。我们使用 "地板 "函数来确保样本数为4096。
n = log(length(x))/log(2); %Adjusting length to make it power of 2
n = floor(n);
x = x(1:2^n);
让我们加入噪音，使这个信号变得嘈杂，并找到WHT系数。
x = x + 0.1.*randn(1, length(x)); %Adding noise
y = fwht(x); %Walsh hadamard transform(WHT)
randn 函数向信号添加随机噪声，fwt 是Matlab内置的函数，用于寻找前向沃尔什哈达玛变换。
由于我们想将其可视化，我们应该备份我们的系数。我们将系数复制到一个不同的变量，yo ，以获得系数的备份。
y0=y; % backup
因为我们只需要系数的首字母，所以我们丢弃这些系数。这意味着我们只取小于500的系数，用它们来重新构建我们的信号，使用ifwht 函数。
y(500:end) = 0; %Removing higher coefficients
xr = ifwht(y); %signal reconstruction using inverse WHT
现在绘制信号图，使输出可视化。
subplot(221) %plot for the noisy ecg signal
plot(x); title('noisyECG signal');
subplot(222) %plot for the hadamard coefficients
plot(abs(y0), 'k'); title('Hadamard transform co-efficient')
subplot(223) % Truncated hadamard transform coefficients
plot(abs(y), 'k'); title('Truncated Hadamard transform co-efficient')
subplot(224) % filtered ecg signal
plot(x); title('filtered ECG signal');
让我们用下面的代码为原始信号和过滤后的信号在同一图中创建一个图，以更清楚地看到这一点。
figure %plot of original and filtered ecg signals.
plot(x); hold on;
plot(xr, 'r', 'LineWidth', 1)
title('Original and filtered ECG signal')
x 是原始信号，而 是过滤后的信号。 表示过滤后的信号被绘制成红色，线宽为1。xr r
输出结果如下图所示。
使用FWHT实现图像压缩的Matlab代码
这种方法具有非常好的能量压缩特性。能量压实特性意味着极少数系数具有最大部分的能量。因此，你可以从变换矩阵中忽略或清空大量的协程。这就是为什么你可以实现输入图像的压缩。
%program for 2D Hadamard transform
[filename, pathname] = uigetfile('*.*', 'Select your greyscale image');
filewithpath = strcat(pathname, filename);
img = imread(filewithpath);
[r, c] = size(img); %getting image size
然后，我们使用double 函数使我们的图像成为双格式。
imgg = double(img);
现在用Matlab的内置函数对信号进行WHT处理。Matlab没有能力对二维信号进行转换。所以我们用一维变换的方法进行列和行的变换。
%Forward WHT
yc = fwht(imgg); %Column wise operation
yr = fwht(yc'); %Row wise operation
y = yr'; %WHT coefficients
yo=y; %Co-efficient backup
yc 给出了列变换， ，给出了等于WHT系数的行变换。yr
丢弃你的系数。它是通过使用下面的代码完成的。
y(256:r, 256:c) = 0; %truncating WHT coefficients
然后你对丢弃的系数进行反WHT，得到最终的信号。另外，在这里，Matlab不能对二维信号进行反WHT运算。因此，我们以列和行的方式进行，如下图所示。
%inverse WHT
Ir1 = ifwht(y); %column wise operation
Ir2 = ifwht(Ir1'); %Row wise operation
imgr = Ir2'; %recorvered image
由于反WHT的输出图像是二维格式的，我们将其转换为uint8 ，并将图像写入当前目录，名称为imgcompressed.jpg 。
imgr8 = uint8(imgr);
imwrite(imgr8, 'imgcompressed.jpg'); % writing image
现在比较原始图像和压缩后的图像的峰值信噪比（PSNR）。这使你能够了解质量的情况。
%Calculate PSNR
mse = sum(sum((imgr-imgg).^2))/(r*c);
maxp = max(max(imgr(:)), max(imgg(:)));
PSNR = 10*log10(double((maxp^2)/mse));
fprintf('\n PSNR=%f\n', PSNR);
mse 是信号的均方根，而 是图像的最大像素。因此，找到 的公式是 。maxp PSNR 10log10(maxp^2)/mse
为了使它更明显，我们用下面的代码绘制输出。
%displaying results
subplot(221) %plot of the original image
imshow(img); title('Original image')
subplot(222) %plot of the hadamard transform coefficients
imshow(imadjust(yo)); title('Hadamard Transform co-efficients')
subplot(223) %plot of the truncated hadamard transform coefficients
imshow(imadjust(y)); title('Truncated Hadamard transform co-efficients')
subplot(224) % Plot of the finally compressed image
imshow(imgr8); title('compressed image');
输出结果如下图所示。
我们可以看到，输入的图像和输出的图像在外观上没有区别。这就是这个方法的工作效率。
要看你的图像是否被压缩了，你可以在文件夹中比较输入和输出的大小。图像从113kb压缩到32kb，这是一个很好的压缩。
结论
Walsh-Hadamard算法是Walsh和Hadamard两种算法的结合，正如我们所看到的，它是一种有效的信号过滤和图像压缩方法。这是因为这些算法在矩阵的基础上工作。
由于信号和图像是矩阵的组合，这种方法发现很容易有效和快速地处理它们。另外，Matlab中的内置函数使整个过程变得简单，因为在没有Matlab的情况下编写整个算法的代码是很累的。
