P某的备忘录

传统方法下图像处理的数学原理整理 (1)

传统方法下图像处理的数学原理整理 (1)

封面图片 credit: From Digital Image Processing 4E, Global Edition

众所周知在图像处理领域,神经网络还没有大规模运用时,支撑这个领域半边天的是各种数学原理(虽然现在数学也很重要)。这一篇文章里面将会滚动整理目前本P接触到的图像处理数学原理方法。

在讨论数学问题之前,我们需要明白,到底发明许许多多的变换方法的意义何在。明明在空域/时域足够直观的东西,为什么要额外发明一个变换,让人们乍一看不知道是啥。

从现代数学的眼光来看,各种变换都是将满足一定条件的某个函数(图片)表示成基函数的线性组合或者积分。这些变换大抵在数学中都属于「调和分析」种类。「分析」二字,可以解释为深入的研究。从字面上来看,「分析」二字,实际就是「条分缕析」而已。它通过对函数的「条分缕析」来达到对复杂函数的深入理解和研究。从哲学上看,「分析主义」和「还原主义」,就是要通过对事物内部适当的分析达到增进对其本质理解的目的。比如近代原子论试图把世界上所有物质的本源分析为原子,而原子不过数百种而已,相对物质世界的无限丰富,这种分析和分类无疑为认识事物的各种性质提供了很好的手段。

也就是说,人们发明如此多种变换方式,是想更加清晰地「提取」原来图像中的各种成分。各种变换区别在于对图像的各种「成分」描述,以及「基元」的不同。

Disclaimer: 以下讨论均局限于二维情况,FF表示变换后的函数,ff表示变换前的函数

离散傅里叶变换(Discrete Fourier Transform, DFT)

DFT的正/逆变换如下:

F(u,v)=x=0M1y=0N1f(x,y)ej2π(ux/M+vy/N)F(u, v)=\sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) e^{-j 2 \pi(u x / M+v y / N)}

f(x,y)=1MNu=0M1v=0N1F(u,v)ej2π(ux/M+vy/N)f(x, y)=\frac{1}{M N} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) e^{j 2 \pi(u x / M+v y / N)}

通过以上变换,我们就将图像拆为了频域成分。

在离散傅里叶变换中,最重要的是如何采样,于是引申出采样定理

奈奎斯特频率(英语:Nyquist frequency)是离散信号系统采样频率的一半,因瑞典裔美国工程师哈里·奈奎斯特(Harry Nyquist)或奈奎斯特-香农采样定理得名。采样定理指出,只要离散系统的奈奎斯特频率高于被采样信号的最高频率或带宽,就可以避免混叠现象。

从理论上说,即使奈奎斯特频率恰好大于信号带宽(但不可相等),也足以通过信号的采样重建原信号。但是,重建信号的过程需要以一个低通滤波器或者带通滤波器将在奈奎斯特频率之上的高频分量全部滤除,同时还要保证原信号中频率在奈奎斯特频率以下的分量不发生畸变,而这是不可能实现的。在实际应用中,为了保证抗混叠滤波器的性能,接近奈奎斯特频率的分量在采样和信号重建的过程中可能会发生畸变。因此实际应用中信号带宽并不能无限接近奈奎斯特频率,具体的情况要看所使用的滤波器的性能。
——来自 Wikipedia

由上所知,采样周期要小于原始周期的12\frac{1}{2},采样频率应该大于原始频率的2倍。否则将发生混叠(又称相位/频率模糊,或者叫锯齿更接地气一点?)。

为什么在数字图象已经离散的情况下,我们仍需要采样定理?由于大多数高精度(远小于人眼分辨率)图像可以认为是「连续」的,如果从这种图上面进行一些压缩或者采集小图,确实需要精确设置取样点。

讲完了这么多,以下奉上离散傅里叶变换的对称 cheetsheet 一张:

从图中可见,若原函数为纯实函数,则对应傅里叶变换函数为偶函数,且实部为偶函数,虚部为奇函数;若原函数为纯虚函数,则对应傅里叶变换为奇函数,且实部为奇函数,虚部为偶函数。后续性质可类推。

除去对称性质,DFT还有以下性质:
平移性质:

f(x,y)ej2π(u0x/M+v0y/N)F(uu0,vv0)f(xx0,yy0)F(u,v)ej2π(x0u/M+y0v/N)\begin{array}{c} f(x, y) e^{j 2 \pi\left(u_{0} x / M+v_{0} y / N\right)} \Leftrightarrow F\left(u-u_{0}, v-v_{0}\right) \\ f\left(x-x_{0}, y-y_{0}\right) \Leftrightarrow F(u, v) e^{-j 2 \pi\left(x_{0} u / M+y_{0} v / N\right)} \end{array}

可见无论是在哪个域平移,均不会改变对应域中的幅值。

旋转性质:

f(r,θ+θ0)F(ω,φ+θ0)f\left(r, \theta+\theta_{0}\right) \Leftrightarrow F\left(\omega, \varphi+\theta_{0}\right)

可见无论是在哪个域旋转,另一个域均会旋转相应角度。

周期性(显然,此处略)

但一个利用周期性把将频域信号“居中”的有趣过程是进行如下变换:

f(x,y)(1)x+yF(uM/2,vN/2)f(x, y)(-1)^{x+y} \Leftrightarrow F(u-M / 2, v-N / 2)

为啥好好的图要居中呢?「For Displaying Purpose」,好咯(


讲完数学原理,来到了「应用实例」模块。

我们学了这么多傅里叶变换,自然是要拿来在频域滤波的,但在滤波之前,首先要解决周期性DFT带来的问题:

假设我们不管DFT带来的周期性问题,那么贸然卷积将会造成「想当然地」周期延拓,这点对于图像尤其重要。因为我们知道大多数情况下处理图像均为「一张」独立的照片,如果不适当限制,傅立叶卷积会默认如右上角图的情况。可实际情况下我们只有粗体那一段有效距离,所以「padding」问题应运而生。目前通用的做法是「padding」两倍如下:

P=2MQ=2N\begin{array}{c} P=2 M \\ Q=2 N \end{array}

其中M,NM,N指的是原始图片尺寸,P,QP,Q指的是padding的总尺寸大小。综上,总频域滤波步骤如下:

采用MATLAB内置的rice.png,使用Butterworth高通滤波器H(u,v)=11+[D0/D(u,v)]2nH(u, v)=\frac{1}{1+\left[D_{0} / D(u, v)\right]^{2 n}}进行滤波操作,其中D0,nD_0, n为自定义参数,D(u,v)D(u,v)指的是flter kernel某点到填充0后的图像中心的欧几里德距离。我们选取D0=60,n=2D_0=60, n=2这一组参数编写对应MATLAB代码如下:

% MATLAB Code | Butterworth High Pass Filter 
clear      
% Reading input image : input_image 
input_image = imread('rice.png'); 
  
% Saving the size of the input_image in pixels- 
% M : no of rows (height of the image) 
% N : no of columns (width of the image) 
[M, N] = size(input_image); 

% Zero padding to avoid aliasing
P = 2*M;
Q = 2*N;
fp = zeros(P,Q); % initialize zero padding matrix 
fp(1:M, 1:N) = input_image;

% Getting Fourier Transform of the input_image 
% using MATLAB library function fft2 (2D fast fourier transform) 
FT_img = fft2(double(fp)); 
 
% Assign the order value 
n = 2; % one can change this value accordingly 
  
% Assign Cut-off Frequency 
D0 = 10; % one can change this value accordingly 
  
% Designing filter and filtering mask
Hp = zeros(P,Q);
for u = 1:P
    for v = 1:Q
        D = sqrt((u-P/2).^2+(v-Q/2).^2); % Calculating Euclidean Distance 
        Hp = 1./(1 + (D0./D).^(2*n)); 
    end
end

  
% Convolution between the Fourier Transformed  
% image and the mask 
G = Hp.*FT_img; 
  
% Getting the resultant image by Inverse Fourier Transform  
% of the convoluted image using MATLAB library function   
% ifft2 (2D inverse fast fourier transform)    
Gp = real(ifft2(double(G)));  
output_image = Gp(1:M,1:N);
% Displaying Input Image and Output Image  
subplot(2, 1, 1), imshow(input_image),title('original')  
subplot(2, 1, 2), imshow(output_image, [ ]), title('filtered'); 

结果如图:


傅里叶变换是如此实用,但多少存在些针对性不足的问题。比如图像上输入均为实数,这对数学上的傅立叶变换而言有种「大炮打蚊子」的感觉。再者,航天器上的计算终端由于目前条件限制可不能携带高计算性能的处理器(指频繁求和差分等操作)所以接下来的几种变换应运而生。需要注意的是,以下操作都是基于矩阵分块运算,与之前的傅立叶变换不一致,这也是为了起到节省运算内存的目的。这里我们插一段矩阵运算数学基础。

为了说明简单,我们把一维傅立叶变换写成如下形式:

T(u)=x=0N1f(x)r(x,u)T(u)=\sum_{x=0}^{N-1} f(x) r(x, u)

这里的f(x)f(x)是原函数,xx是坐标变量(亦称空间变量),uu是相对应频域内的变量,r(x,u)r(x,u)指的是前向变换掩模(forward transform kernel)

易知逆变换形式如下:

f(x)=u=0N1T(u)s(x,u)f(x)=\sum_{u=0}^{N-1} T(u) s(x, u)

有了这些准备后,我们将式中三个符号写成矩阵形式:

f=[f(0)f(1)f(N1)]=[f0f1fN1]t=[T(0)T(1)T(N1)]=[t0t1tN1]su=[s(0,u)s(1,u)s(N1,u)]=[su,0su,1su,N1] for u=0,1,,N1\begin{array}{c} \mathbf{f}=\left[\begin{array}{c} f(0) \\ f(1) \\ \vdots \\ f(N-1) \end{array}\right]=\left[\begin{array}{c} f_{0} \\ f_{1} \\ \vdots \\ f_{N-1} \end{array}\right] \\ \mathbf{t}=\left[\begin{array}{c} T(0) \\ T(1) \\ \vdots \\ T(N-1) \end{array}\right]=\left[\begin{array}{c} t_{0} \\ t_{1} \\ \vdots \\ t_{N-1} \end{array}\right] \\ \mathbf{s}_{u}=\left[\begin{array}{c} s(0, u) \\ s(1, u) \\ \vdots \\ s(N-1, u) \end{array}\right]=\left[\begin{array}{c} s_{u, 0} \\ s_{u, 1} \\ \vdots \\ s_{u, N-1} \end{array}\right] \text { for } u=0,1, \ldots, N-1 \end{array}

那么以上的一维傅立叶变换即可用线性代数中的内积符号<><>表示:

T(u)=su,f for u=0,1,,N1T(u)=\left\langle\mathbf{s}_{u}, \mathbf{f}\right\rangle \quad \text { for } u=0,1, \ldots, N-1

为了使得表达式简洁,我们把su,xs_{u,x}中的uu略去,于是构造变换矩阵AA

A=[s0Ts1TsN1T]=[s0s1sN1]T\mathbf{A}=\left[\begin{array}{c} \mathbf{s}_{0}^{T} \\ \mathbf{s}_{1}^{T} \\ \vdots \\ \mathbf{s}_{N-1}^{T} \end{array}\right]=\left[\begin{array}{llll} \mathbf{s}_{0} & \mathbf{s}_{1} & \ldots & \mathbf{s}_{N-1} \end{array}\right]^{T}

又于是乎,上述表示可以完全变成矩阵形式:

亦有:

t=Af\mathbf{t}=\mathbf{A} \mathbf{f}

观察得到:

于是我们得到了一维下的矩阵变换表达式:

t=Aff=ATt\begin{array}{l} \mathbf{t}=\mathbf{A} \mathbf{f} \\ \mathbf{f}=\mathbf{A}^{T} \mathbf{t} \end{array}

同理可得,二维下的表达式有:

T=AFATF=ATTA\begin{array}{l} \mathbf{T}=\mathbf{A} \mathbf{F} \mathbf{A}^{T} \\ \mathbf{F}=\mathbf{A}^{T} \mathbf{T} \mathbf{A} \end{array}

以上即为矩阵分块处理的数学推导。从以下方法开始,阐述的方法都是描述上述提到的变换矩阵AA


离散余弦变换(Discrete Cosine Transform, DCT)

首先抛出M×NM\times N矩阵 AA 的二维 DCT及其逆变换的定义:

BpqB_{pq} 称为 AADCT 系数,此处的NN即是指DCT矩阵的边长。

这种变化实质上取了傅立叶变化里面的实数成分,再硬造了一个让原函数偶延拓的系数,为了讨论简便,我们在一维情况下推导。


首先用欧拉公式展开DFT的exp\exp,祭出DFT的原始表达式:

X[k]=n=0N1x[n](cos2πknN)jn=0N1x[n]sin(2πknN)X[k]=\sum_{n=0}^{N-1} x[n]\left(\cos \frac{2 \pi \mathrm{kn}}{N}\right)-j \sum_{n=0}^{N-1} x[n] \sin \left(\frac{2 \pi k n}{N}\right)

显而易见的,DFT变换的结果,实数部分由n=0N1x[n](cos2πknN)\sum_{n=0}^{N-1} x[n]\left(\cos \frac{2 \pi \mathrm{kn}}{N}\right)承包了,而虚数部分则由jn=0N1x[n]sin(2πknN)j \sum_{n=0}^{N-1} x[n] \sin \left(\frac{2 \pi k n}{N}\right) 负责,为了描述简便,记cos(2πknN)=cos(kt),sin(2πknN)=sin(kt)\cos \left(\frac{2 \pi \mathrm{kn}}{N}\right)=\cos (k t), \sin \left(\frac{2 \pi \mathrm{kn}}{N}\right)=\sin (k t)

那么DFT变换的实部有:

Re[k]=n=0N1x[n]cos(kt)\operatorname{Re}[k]=\sum_{n=0}^{N-1} x[n] \cos (k t)

虚部有:

Im[k]=n=0N1x[n]sin(kt)\operatorname{Im}[k]=\sum_{n=0}^{N-1} x[n] \sin (k t)

显然的,cos\cos是一个偶函数,sin\sin是一个奇函数,因此有:

Re[k]=Re[k],Im[k]=Im[k]\operatorname{Re}[k]=\operatorname{Re}[-k], \operatorname{Im}[k]=-\operatorname{Im}[k]

因此,当 x[n]x[n]是一个实数函数时,其频域的实部是偶函数,虚部是一个奇函数。

进一步地,假如原信号x[n]x[n]是一个全是实数的偶函数信号,那么显然,偶函数乘偶函数还是偶函数,奇函数乘偶函数还是奇函数,因此:

Im[k]=n=0N1x[n]sin(kt)=0\operatorname{Im}[k]=\sum_{n=0}^{N-1} x[n] \sin (k t)=0

于是虚部消失了。因此,当原时域信号是一个实偶信号时,我们就可以把DFT写成:

X[k]=\sum_{n=0}^{N-1} x[n]\left(\cos \frac{2 \pi \mathrm{kn}}{N}\right)\right)

上式就是DCT变换的核心,但这还不能投入生产运用,毕竟实际中的信号(图片)虽然是实数,但不可能那么巧全是偶信号。于是发明DCT的人通过数学手段“造”了一个偶信号出来。

设一长度为NN的实数离散信号{x[0],x[1]x[N1]}\{x[0], x[1] \ldots \ldots x[N-1]\} ,首先,我们先将这个信号长度扩大成原来的两倍,即长度变为2N2N,定义新信号x[m]x^{\prime}[m]为:

x[m]=x[m](0mN1)x[m]=x[m1](Nm1)\begin{array}{l} x^{\prime}[m]=x[m](0 \leq m \leq N-1) \\ x^{\prime}[m]=x[-m-1](-N \leq m \leq-1) \end{array}


也就是说,信号变成了这个样子。可还是不对,这玩意也不关于原点对称啊,于是我们把信号向右进行平移12\frac{1}{2}个单位长度:

于是,这个信号相应DFT变成了:

X[k]=m=N+12N12x[m12]ej2πmk2NX[k]=\sum_{m=-N+\frac{1}{2}}^{N-\frac{1}{2}} x^{\prime}\left[m-\frac{1}{2}\right] e^{\frac{-j 2 \pi m k}{2 N}}

再用欧拉公式展开,丢掉虚数部分,我们有:

X[k]=m=N+12N12x[m12]ej2πmk2N=m=N+12N12x[m12]cos(2πmk2N)X[k]=\sum_{m=-N+\frac{1}{2}}^{N-\frac{1}{2}} x^{\prime}\left[m-\frac{1}{2}\right] e^{\frac{-j 2 \pi m k}{2 N}}=\sum_{m=-N+\frac{1}{2}}^{N-\frac{1}{2}} x^{\prime}\left[m-\frac{1}{2}\right] \cos \left(\frac{2 \pi m k}{2 N}\right)

到这一步,虽然信号问题解决了,但mm算出来居然是一个小数还带负数。因此,上式我们还需要进一步变形。

首先我们知道,这个序列是一个偶对称序列,因此根据对称性:

X[k]=m=N+12N12x[m12]cos(2πmk2N)=2m=12N12x[m12]cos(2πmk2N)X[k]=\sum_{m=-N+\frac{1}{2}}^{N-\frac{1}{2}} x^{\prime}\left[m-\frac{1}{2}\right] \cos \left(\frac{2 \pi m k}{2 N}\right)=2\sum_{m=\frac{1}{2}}^{N-\frac{1}{2}} x^{\prime}\left[m-\frac{1}{2}\right] \cos \left(\frac{2 \pi m k}{2 N}\right)

然后,记n=m12n=m-\frac{1}{2},并代入上式:

X[k]=2n=0N1x[n]cos(2π(n+12)k2N)=2n=0N1x[n]cos((n+12)πkN)X[k]=2 \sum_{n=0}^{N-1} x^{\prime}[n] \cos \left(\frac{2 \pi\left(n+\frac{1}{2}\right) k}{2 N}\right)=2 \sum_{n=0}^{N-1} x^{\prime}[n] \cos \left(\frac{\left(n+\frac{1}{2}\right) \pi k}{N}\right)

至此,我们距离常用DCT的标准式已经非常的接近了,现在问题就是那个标准式里的那个系数c(u)c(u)是个啥玩意。实际上这个c(u)c(u)如果在函数计算中,加不加都无所谓,实际上在DFT变换中,这个值也是存在的,看不出来是因为它常常取1,因此没有再进一步写出来,实际上,这个值因为一些工程学上的意义,DFT中也常常会取1N{\frac{1}{N}}1N\sqrt{\frac{1}{N}}

那么,DCT中它的出现,主要是为了在DCT变换变成矩阵运算的形式时,将该矩阵正交化以便于进一步的计算,那么,这个系数就应该取12N\sqrt{\frac{1}{2N}}

最终乘上这个系数,于是就有了DCT现在这个样子,推导完成。


从应用层面上,DCT变换较DFT变换具有更好的频域能量聚集度,那么对于那些不重要的频域区域和系数就能够直接裁剪掉,因此,DCT变换非常适合于图像压缩算法的处理,例如现在大名鼎鼎的jpeg就是使用了DCT作为图像压缩算法。如何体现呢?我们通过以下的MATLAB代码演示:

% MATLAB Code | Discrete Cosine Transform
clear      

I = imread('cameraman.tif');
I = im2double(I);

% Perform DCT
T = dctmtx(8);
dct = @(block_struct) T * block_struct.data * T'; %dct transform
B = blockproc(I,[8 8],dct); % 8x8 subregion processing
B_0 =B; %temporarliy store DCT matrix data


% Method1: Left out redundant information
B(abs(B)<0.3) = 0;
B_masked = B;

% Perform IDCT

invdct = @(block_struct) T' * block_struct.data * T;
I1 = blockproc(B_masked,[8 8],invdct);

% Method2: Directly process DCT matrix
mask = [1   1   1   1   0   0   0   0
        1   1   1   0   0   0   0   0
        1   1   0   0   0   0   0   0
        1   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0
        0   0   0   0   0   0   0   0];
masked = blockproc(B_0,[8 8],@(block_struct) mask .* block_struct.data);

% Perform IDCT
invdct = @(block_struct) T' * block_struct.data * T;
I2 = blockproc(masked,[8 8],invdct);

% Show processed image

figure, subplot(141), imshow(I,[]), title('Original image');
subplot(142), imshow(B,[]), colorbar, title('DCT result');
subplot(143), imshow(I1,[]), title('Method1 processed image');
subplot(144), imshow(I2,[]), title('Method2 processed image');

以上代码通过两种方式展示了DCT处理的方式,通过DCT变换结果可以发现,图片的“能量”,亦即响应,几乎都集中在了大于0的部分,且大于0的部分均集中在左上角(白点区域,注意原图是分块处理的)。于是方式1把响应小于0.3的部分直接清除,方式2改变了DCT处理的8×88\times 8矩阵,把右下角那一堆全部不要了。
结果如图:

可以发现图像进行了压缩。大概的轮廓仍然可以辨认,这验证了其压缩算法的实用性。

Walsh-Hadamard 变换

阿达马变换(Hadamard transform),或称沃尔什-阿达玛转换,是一种广义傅立叶变换(Fourier transforms),作为变换编码的一种在影片编码当中使用有很久的历史。这种变换干脆就放弃了积分/差分运算。由于计算简便,该变换经常被用在航天器上。

在开始变换之前,我们引入矩阵的克罗内克积


数学上,克罗内克积(英语:Kronecker product)是两个任意大小的矩阵间的运算,表示为\bigotimes。克罗内克积是外积从向量到矩阵的推广,也是张量积在标准基下的矩阵表示。

尽管没有明显证据证明德国数学家利奥波德·克罗内克是第一个定义并使用这一运算的人,克罗内克积还是以其名字命名。在历史上,克罗内克积曾以Johann Georg Zehfuss名字命名为Zehfuss矩阵。

如果AA是一个 m×nm × n的矩阵,而BB是一个 p×qp × q的矩阵,克罗内克积$ A\otimes B$则是一个 $mp × nq $的分块矩阵:

AB[a11Ba1nBam1BamnB]mk×nl\mathbf{A} \otimes \mathbf{B} \triangleq\left[\begin{array}{lll} a_{11} \mathbf{B} & \cdots & a_{1 n} \mathbf{B} \\ \cdots & \cdots & \cdots \\ a_{m 1} \mathbf{B} & \cdots & a_{m n} \mathbf{B} \end{array}\right]_{m k \times n l}

进一步地有:

AB=[a11b11a11b12a11b1qa1nb11a1nb12a1nb1qa11b21a11b22a11b2qa1nb21a1nb22a1nb2qa11bp1a11bp2a11bpqa1nbp1a1nbp2a1nbpqam1b11am1b12am1b1qamnb11amnb12amnb1qam1b21am1b22am1b2qamnb21amnb22amnb2qam1bp1am1bp2am1bpqamnbp1amnbp2amnbpq]A \otimes B=\left[\begin{array}{cccccccccc} a_{11} b_{11} & a_{11} b_{12} & \cdots & a_{11} b_{1 q} & \cdots & \cdots & a_{1 n} b_{11} & a_{1 n} b_{12} & \cdots & a_{1 n} b_{1 q} \\ a_{11} b_{21} & a_{11} b_{22} & \cdots & a_{11} b_{2 q} & \cdots & \cdots & a_{1 n} b_{21} & a_{1 n} b_{22} & \cdots & a_{1 n} b_{2 q} \\ \vdots & \vdots & \ddots & \vdots & & & \vdots & \vdots & \ddots & \vdots \\ a_{11} b_{p 1} & a_{11} b_{p 2} & \cdots & a_{11} b_{p q} & \cdots & \cdots & a_{1 n} b_{p 1} & a_{1 n} b_{p 2} & \cdots & a_{1 n} b_{p q} \\ \vdots & \vdots & & \vdots & \ddots & & \vdots & \vdots & & \vdots \\ \vdots & \vdots & & \vdots & & \ddots & \vdots & \vdots & & \vdots \\ a_{m 1} b_{11} & a_{m 1} b_{12} & \cdots & a_{m 1} b_{1 q} & \cdots & \cdots & a_{m n} b_{11} & a_{m n} b_{12} & \cdots & a_{m n} b_{1 q} \\ a_{m 1} b_{21} & a_{m 1} b_{22} & \cdots & a_{m 1} b_{2 q} & \cdots & \cdots & a_{m n} b_{21} & a_{m n} b_{22} & \cdots & a_{m n} b_{2 q} \\ \vdots & \vdots & \ddots & \vdots & & & \vdots & \vdots & \ddots & \vdots \\ a_{m 1} b_{p 1} & a_{m 1} b_{p 2} & \cdots & a_{m 1} b_{p q} & \cdots & \cdots & a_{m n} b_{p 1} & a_{m n} b_{p 2} & \cdots & a_{m n} b_{p q} \end{array}\right]


讲完了克罗内克积,我们现在开始定义 Hadamard 矩阵:

H112[1111]Hn=H1Hn1=12[Hn1Hn1Hn1Hn1]\begin{array}{c} \mathbf{H}_{1} \triangleq \frac{1}{\sqrt{2}}\left[\begin{array}{rr} 1 & 1 \\ 1 & -1 \end{array}\right] \\ \mathbf{H}_{n}=\mathbf{H}_{1} \otimes \mathbf{H}_{n-1}=\frac{1}{\sqrt{2}}\left[\begin{array}{cr} \mathbf{H}_{n-1} & \mathbf{H}_{n-1} \\ \mathbf{H}_{n-1} & -\mathbf{H}_{n-1} \end{array}\right] \end{array}

可以看到 Hadamard 矩阵是递归定义的,我们举n=2,3n=2, 3情况为例:

H2=H1H1=12[H1H1H1H1]=14[1111111111111111]\mathbf{H}_{2}=\mathbf{H}_{1} \otimes \mathbf{H}_{1}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc} \mathbf{H}_{1} & \mathbf{H}_{1} \\ \mathbf{H}_{1} & -\mathbf{H}_{1} \end{array}\right]=\frac{1}{\sqrt{4}}\left[\begin{array}{rrrr} 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & -1 & 1 \end{array}\right]

矩阵H3{\bf H}_3右边的第一列是N=8N=8行的索引标记,第二列代表每行正负号改变的次序,称序列。序列与频率相似但又不同,它衡量的是非周期性信号的变化率该序列亦称作 Walsh 序列,对应情况如下图(credit: Harvey Mudd College)

容易证明,Hadamard 矩阵为实正交对称矩阵。

另一种定义 Hadamard 矩阵第kkmm列矩阵值的方法是:

h[k,m]=(1)i=0n1kimi=i=0n1(1)kimi=h[m,k](k,m=0,1,,N1)h[k, m]=(-1) \sum_{i=0}^{n-1} k_{i} m_{i}=\prod_{i=0}^{n-1}(-1)^{k_{i} m_{i}}=h[m, k] \quad(k, m=0,1, \cdots, N-1)

其中:

k=i=0n1ki2i=(kn1kn2k1k0)2(ki=0,1)m=i=0n1mi2i=(mn1mn2m1m0)2(mi=0,1)\begin{array}{c} k=\sum_{i=0}^{n-1} k_{i} 2^{i}=\left(k_{n-1} k_{n-2} \cdots k_{1} k_{0}\right)_{2} \quad\left(k_{i}=0,1\right) \\ m=\sum_{i=0}^{n-1} m_{i} 2^{i}=\left(m_{n-1} m_{n-2} \cdots m_{1} m_{0}\right)_{2} \quad\left(m_{i}=0,1\right) \end{array}

注意,这里的n=log2Nn=\log_{2}^{N}。可以看见此处二进制定义了两个变量:kkmm,正好,二进制在机器语言里面非常常见,所以就很通用了。

我们举个例子:在8×88\times 8的 Hadamard 矩阵中,位于(3,5)(3,5)的坐标值为:

h[3,5]=i=02(1)kimi,k=011,m=101h[3,5]=\prod_{i=0}^{2}(-1)^{k_i m_i}, k=011, m=101

于是

h[3,5]=(1)(11)(1)(10)(1)(01)=111=1h[3,5]=(-1)^(1\cdot 1)\cdot (-1)^(1\cdot 0)\cdot (-1)^(0\cdot 1)=-1\cdot 1\cdot 1=-1


应用层面上,Walsh-Hadamard 变换亦常用于压缩,此处举一维心电图(Electrocardiography, ECG)信号压缩为例:

% MATLAB Code | Walsh-Hadamard Transform
clear

addpath('/path/to/WalshHadamardTransformForCompressionOfECGSignalsExample/') % addpath for ecg signal

% Input data
x1 = ecg(512);                    % Single ecg wave
x = repmat(x1,1,8);                 
x = x + 0.1.*randn(1,length(x));  % Noisy ecg signal

% Apply Walsh-Hadamard Transform
y = fwht(x);                      % Fast Walsh-Hadamard transform
figure
subplot(2,1,1)
plot(x)
xlabel('Sample index')
ylabel('Amplitude')
title('ECG Signal')
subplot(2,1,2)
plot(abs(y))
xlabel('Sequency index')
ylabel('Magnitude')
title('WHT Coefficients')

% Process out the low magnitude sequency component
y(1025:length(x)) = 0;            % Zeroing out the higher coefficients    
xHat = ifwht(y);                  % Signal reconstruction using inverse WHT  
figure
plot(x)
hold on
plot(xHat,'r')
xlabel('Sample index')
ylabel('ECG signal amplitude')
legend('Original Signal','Reconstructed Signal')

压缩信号与原始信号如下:

分别对应信息如下:

  Name                       Size              Bytes   Class     Attributes

  reconstructed             1x1025            8200   double              
  original                  1x4096            32768  double              

可见压缩率达到了惊人的4:1。

Powered by Gridea. CopyRight 2018-2020 MrPeterJin's Media Dept.