线性系统的基本理论与运算

线性系统与非线性系统

设系统的特性可以表示成对输入图像进行T运算, 并令$f_1(x,y)$, $f_2(x,y)$与$T[f_1(x,y)]$,$T[f_2(x,y)]$分别代表两对不同的输入和输出图像, 则当系统满足:

$T\left[f_{1}(x, y)\right]+T\left[f_{2}(x, y)\right]=T\left[f_{1}(x, y)+f_{2}(x, y)\right]$

关系时, 称系统具有叠加性.当系统满足:

$T\left[k f_{1}(x, y)\right]=k T\left[f_{1}(x, y)\right]$

关系时, 称系统具有齐次性.

同时满足叠加性和齐次性的系统称为线性系统.由于图像是二维的, 所以这样的线性系统称为二维线性系统, 由以上两个公式所定义的运算称为二维线性运算. 显然,二维线性系统应一般地满足:

$T\left[\sum k_{i} f_{i}(x, y)\right]=\sum k_{i} T\left[f_{i}(x, y)\right]$

凡是不满足叠加性和其次性的系统都属于非线性系统.

线性系统的基本理论

移不变系统

当系统的单位脉冲输入为$\boldsymbol{\delta}(\boldsymbol{x}-\boldsymbol{a}, \boldsymbol{y}-\boldsymbol{\beta})$, 也即输入的单位脉冲函数延时了$\alpha, \beta$单位时, 输出为$\boldsymbol{h}(\boldsymbol{x}-\boldsymbol{a}, \boldsymbol{y}-\beta)$, 即 输出结果性态不变, 仅在位置上延迟了$\alpha, \beta$单位, 则称,这样的系统为移不变系统

显然, 对于移不变系统来说, 系统的输出仅与输入函数的性态有关, 而与输入函数作用的起点无关.且:

$T[\delta(x-\alpha, y-\beta)]=h(x-\alpha, y-\beta)$

线性移不变系统

如果一个系统是线性系统, 又是移不变系统, 则称该系统是线性移不变系统.

对于一个二维线性移不变系统$h(x,y)$, 设其输入为$f(x,y)$, 输出为$g(x,y)$, 线性以不变系统的运算为T,则有:

$\begin{aligned} g(x, y) &=T[f(x, y)] \\ &=T\left[\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\alpha, \beta) \delta(x-\alpha, y-\beta) d \alpha d \beta\right] \\ &\left.=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} T[f(\alpha, \beta) \delta(x-\alpha, y-\beta)] d \alpha d \beta\right] \\ &=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\alpha, \beta) T[\delta(x-\alpha, y-\beta)] d a d \beta \\ &=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\alpha, \beta) h(x-\alpha, y-\beta) d \alpha d \beta \\ &=f(x, y) * h(x, y) \end{aligned}$

即: 线性移不变系统的输出等于系统的输入与系统的脉冲响应的卷积.

同理有:

$\begin{aligned} g(x, y) &=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x-\alpha, y-\beta) h(\alpha, \beta) d \alpha d \beta \\ &=h(x, y) * f(x, y) \end{aligned}$

也有:

$h(x, y) * f(x, y)=f(x, y) * h(x, y)$

所以, 二维线性移不变系统的输入, 输出和运算关系可以描述为:

离散傅立叶变换

一维离散傅立叶变换

设f(x)是在时域上的等距离采样得到的N点离散序列, x是离散实变量, u为离散频率变量, 则离散傅立叶变换对定义为:

$F(u)=\frac{1}{\sqrt{N}} \sum_{x=0}^{N-1} f(x) \exp \left[-\frac{j 2 \pi x u}{N}\right] \quad, \quad u=0,1, \cdots, N-1$

$f(x)=\frac{1}{\sqrt{N}} \sum_{u=0}^{N-1} F(u) \exp \left[\frac{j 2 \pi \alpha x}{N}\right] \quad, \quad x=0,1, \cdots, N-1$

其中, $F(u)$为正变换, $f(x)=f^{-1}[F(u)]$为反变换;

$e^{-j \frac{2 \pi}{N} x t}$是正变换核,$e^{j \frac{2 \pi}{N} w}$是反变换核.

根据欧拉公式$e^{\pm i x}=\cos x \pm i \sin x$, 有:

$\exp (-j 2 \pi x u)=\cos 2 \pi u x-j \sin 2 \pi u x$

所以, $F(u)$一般为复数, 并可以写成

$F(u)=R(u)+j I(u)$

其中, R(u)和I(u)分别为$F(u)$的实部和虚部, 指数形式为:

$F(u)=|R(u)|+\exp [-j \phi(u)]$

且:$|F(u)|=\sqrt{R^{2}(u)+I^{2}(u)}$, $\varphi(u)=\arctan \left[\frac{I(u)}{R(u)}\right]$

其中, $\boldsymbol{\Phi}(\boldsymbol{u})$称为相位角, 反映了$f(x)$的相频特性.

二维离散傅立叶变换

设$f(x,y)$是在空间域等间隔采样得到的 $M \times N$的二维离散信号, x和y是离散实变量, u和v是离散频率变量, 则二维离散傅立叶变换对一般地定义为:

$\begin{aligned} F(u, v) &=\sqrt{\frac{1}{M N}} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) \exp \left[-j 2 \pi\left(\frac{x u}{M}+\frac{y v}{N}\right)\right] \\ \end{aligned}$

$(\boldsymbol{u}=\mathbf{0}, \mathbf{1}, \ldots, \boldsymbol{M}-\mathbf{1} ; \boldsymbol{v}=\mathbf{0}, \mathbf{1}, \ldots, \boldsymbol{N}-\mathbf{1})$

$f(x, y)=\sqrt{\frac{1}{M N}} \sum_{u=0}^{M-1 N-1} F(u, v) \exp \left[j 2 \pi\left(\frac{u x}{M}+\frac{v y}{N}\right)\right]$

$(x=0,1, \dots, M-1 ; y=0,1, \dots, N-1)$

在图像处理中有时为了讨论方便, 取M=N, 并考虑到正变换与反变换的对称性, 就将二维离散傅立叶变换对定义为:

$\left.F(u, v)=\frac{1}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \exp F \frac{j 2 \pi(x u+y)}{N}\right]$

$f(x, y)=\frac{1}{N} \sum_{u=0}^{N-1} \sum_{v=0}^{N-1} F(u, v) \exp \left[\frac{j 2 \pi(u x+v y)}{N}\right]$

其中, $x, y, u, v=0,1, \dots, N-1$

与一维的情况类似. 可将二维离散傅立叶变换的频谱和相位角定义为:

$P(u, v) = \left.| F(u, v)\right|^{2}=R^{2}(u, v)+I^{2}(u, v)$

$\phi(u, v)=\arctan [I(u, v) / R(u, v)]$

将二维离散傅立叶变换频谱的平方定义为功率谱, 记为:

$P(u, v) = \left.| F(u, v)\right|^{2}=R^{2}(u, v)+I^{2}(u, v)$

反映了二维离散信号的能量在空间频率域上的分布情况.

图像傅立叶变换的意义

1. 简化计算, 也即傅立叶变换可将空间域中复杂的卷积运算转化为频率域中简单的乘积运算.
2. 对于某些在空间域中难以处理或处理起来比较复杂的问题, 利用傅里叶变换把空间域表示的图像映射到频率域, 再利用频率域滤波或频域分析的方法对其进行处理和分析 ,然后再把其在频域中处理和分析的结果变换回空间域, 从而达到简化处理和分析的目的.
3. 某些只能在频率域处理的特定应用需求, 比如在频率域进行图像特征提取, 数据压缩, 纹理分析, 水印嵌入等.

离散傅立叶变换的若干重要性质

由离散傅立叶变换的反变换式:

$f(x, y)=\frac{1}{N} \sum_{v=0}^{N-1} \sum_{u=0}^{N-1} F(u, v) \exp \left[\frac{j 2 \pi(x u+y v)}{N}\right]$

可知, 由于 u和v 均有0,1,…,N-1的N个可能取值, 所以$f(x,y)$由$N^2$个频率分量组成, 所以每个频率分量都与一个特定的(u,v)值相对应; 且对于某个确定的值来说, 当(x,y)取遍所有可能的值(x=0,1,…,N-1;y=0,1,…,N-1)时, 就可得到对应于该特定值(u,v)的一幅基图像.

$f_{u, v}=\frac{1}{N}\left[\begin{array}{ccc} \exp [j 2 \pi((0 u+0 v) / N)] & \cdots & \exp [j 2 \pi((0 u+(N-1) v) / N)] \\ \vdots & \ddots & \vdots \\ \exp [j 2 \pi(((N-1) u+0 v) / N)] & \cdots & \exp [j 2 \pi(((N-1) u+(N-1) 0 v) / N)] \end{array}\right]$

将二维离散傅立叶变换对写成如下分离形式:

$\left.F(u, v)=\frac{1}{N} \sum_{x=0}^{N-1} \exp \left[\frac{-j 2 \pi \alpha u}{N}\right] \sum_{y=0}^{N-1} f(x, y) \exp \left[\frac{-j 2 \pi y v}{N}\right]\right)$

$f(x, y)=\frac{1}{N} \sum_{u=0}^{N-1} \exp \left[\frac{j 2 \pi u x}{N}\right]\left(\sum_{v=0}^{N-1} F(u, v) \exp \left[\frac{j 2 \pi v y}{N}\right]\right)$

上述的可分离表示形式说明, 可以连续运用两次一维离散傅立叶变换来实现一个二维离散傅立叶变换.

以式$\left.F(u, v)=\frac{1}{N} \sum_{x=0}^{N-1} \exp \left[\frac{-j 2 \pi \alpha u}{N}\right] \sum_{y=0}^{N-1} f(x, y) \exp \left[\frac{-j 2 \pi y v}{N}\right]\right)$为例, 可先沿着y轴方向进行一维的列变换而求得:

$F(x, v)=\frac{1}{\sqrt{N}} \sum_{y=0}^{N-1} f(x, y) \exp \left[\frac{-j 2 \pi v y}{N}\right]$

然后在对$F(x,v)$沿着x方向进行一维的行变换得到最后的结果:

$F(u, v)=\frac{1}{N} \sum_{x=0}^{N-1} F(x, v) \exp \left[\frac{-j 2 \pi u x}{N}\right]$

$\bar{f}=\frac{1}{N^{2}} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y)$

如果将u=v=0代入式:$F(u, v)=\frac{1}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \exp \left[\frac{j 2 \pi(x u+y v)}{N}\right]$可得:

$F(0,0)=\frac{1}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y)$

所以, 一幅图像的灰度平均值可由DFT在原点处的值求得, 即:

$\bar{f}=\frac{1}{N} F(0,0)$

$F(u, v)=F(u+m M, v+n N)$

$(m, n=0,\pm 1,\pm 2, \cdots)$

$F(u, v)=F^{*}(-u,-v)$

$|F(u, v)|=|F(-u,-v)|$

图像傅立叶变换的频谱分析

图像傅立叶频谱关于(M/2,N/2)的对称性

设$f(x,y)$是一幅大小为$M\times N$的图像, 根据离散傅立叶变换的周期性公式:

$F(u, v)=F(u+m M, v+n N)$

有:

$|F(u, v)|=|F(u+M, v+N)|$

再根据离散傅立叶变换的共轭对称公式:

$|F(u, v)|=|F(-u,-v)|$

就可得:

$|F(u, v)|=| F(M-u, N-v)|$

对于u=0 :

• 当v=0时:$|F(0,0)|=F(M, N) |$
• 当v=1时:$|F(0,1)|=F(M, N-1) |$
• 当v=2时:$|F(0,2)|=F(M, N-2) |$
• …
• 当v=N/2时:$|F(0,N/2)|=F(M, N/2) |$

同理, 对于v=0:

• 当u=0时:$|F(0,0)|=F(M, N) |$
• 当u=1时:$|F(1,0)|=F(M-1, N) |$
• 当u=2时:$|F(2,0)|=F(M-2, N) |$
• …
• 当u=M/2时:$|F(M/2,0)|=F(M/2, N) |$

由此可得:

下图是原点坐标位于(0,0)的图像的傅立叶频谱关于(M/2,N/2)对称的两个例子.

1

fft2(I)

图像傅立叶频谱特性及其频谱图

由于图像的傅立叶频谱显示图像的能量分布在图像的四周, 为了便于观察分析, 采用如下的平移策略对图像的傅立叶频谱进行平移, 使得平移后的频谱图中能量集中在图像中央显示.在Matlab中可以使用fftshif(I)实现该操作.

对于式:$f(x, y) \exp \left[j 2 \pi\left(\frac{u_{0} x}{M}+\frac{v_{0} y}{N}\right)\right] \leftrightarrow F\left(u-u_{0}, v-v_{0}\right)$当$\mathbf{u}_{0}=\mathbf{M} / \mathbf{2}, \boldsymbol{v}_{\mathbf{0}}=\mathbf{N} / \mathbf{2}$时, 有:

$\exp \left[j 2 \pi\left(\frac{u_{0} x}{M}+\frac{v_{0} y}{N}\right)\right]=\exp \left[j 2 \pi\left(\frac{M}{2} \cdot \frac{x}{M}+\frac{N}{2} \cdot \frac{y}{N}\right)\right]$

$e^{j \pi(x+y)}=\left(e^{j \pi}\right)^{(x+y)}=(\cos \pi+j \sin \pi)^{(x+y)}=(-1)^{(x+y)}$

也即:

$f(x, y)(-1)^{(x+y)} \leftrightarrow F\left(u-\frac{M}{2}, v-\frac{N}{2}\right)$

也就是说, 上面的频谱图(a),(b)实质上是函数$f(x, y)(-1)^{(x+y)}$的傅立叶频谱图.

傅立叶变换在图像处理中的应用

1. 基本思路是:

   先用$(-1)^{x+y}$乘以图像得$f(x, y)(-1)^{(x+y)}$;然后对其进行傅立叶正变换得到原点在(M/2,N/2)之处的$F(u,v)$; 接着根据图像的频率特性, 利用有关的低通滤波器, 或者高通滤波器等, 对其进行滤波处理; 再将处理的结果进行傅立叶反变换; 最后给反变换的结果在乘以$(-1)^{x+y}$就可得到最终的结果.

2. 典型的应用有:

   • 去除图像噪声
   • 图像数据压缩
   • 图像识别
   • 图像重构和图像描述等

1
2
3
4
5
6
7

I=imread('f:\lena.jpg'); %读入图像
F=fft2(im2double(I)); %FFT
F=fftshift(F); %FFT频谱搬移
F=real(F);
T=log(F+1); % 频谱对数变换
subplot(1,2,1),imshow(I),title('原始图像');
subplot(1,2,2),imshow(T,[]),title('原始图像其频谱图');

图像的离散余弦变换

• 函数的偶对称性使DCT只有实数域变换结果, 不再涉及复数运算, 运算简单, 费时少;
• 又保持了变换域的频率特性
• 与人类视觉系统相适应
• 得到了更加广泛的应用

二维偶DCT

基本思想:

把一个$N\times N$的图像数据矩阵延拓成二维平面上的偶对称阵列. 延拓方式有两种:

1. 围绕图像边缘(但不重叠)将其折叠成对称形式而得到的变换称为偶离散余弦变换;
2. 通过重叠图像的第一列像素和第N-1行像素将其折叠成对称形式得到的变换称为奇离散余弦变换.

下面只介绍偶离散余弦变换.

设$f(x,y)$为一$N\times N$的图像数据阵列, 将$f(x,y)$围绕其左边缘和下边缘不重叠的折叠成偶对称图像, 即下图:

并表示为:

$f_{s}(x, y)=\left\{\begin{array}{ll} f(x, y) & x \geq 0, y \geq 0 \\ f(-x-1, y) & x<0, y \geq 0 \\ f(x,-y-1) & x \geq, y<0 \\ f(-x-1,-y-1) & x<0, y<0 \end{array}\right.$

可见, $2N\times 2N$的新图像的对称中心位于图像中细十字虚线的交叉处, 也就是(-1/2,-1/2)处.

对上述新图像$f_{s}(x, y)$取二维傅立叶变换可得:

$F_{s}(u, v)=\frac{1}{2 N} \sum_{x=-N}^{N-1} \sum_{y=-N}^{N-1} f_{s}(x, y) \exp \left(-\frac{j 2 \pi\left[u\left(x+\frac{1}{2}\right)+v\left(y+\frac{1}{2}\right)\right]}{2 N}\right)$

由于$f_{s}(x, y)$是实对称函数, 欧拉展开式后的正弦项为0值, 所以上式可简化为:

$F_{s}(u, v)=\frac{1}{2 N} \sum_{x=-N}^{N-1} \sum_{y=-N}^{N-1} f_{s}(x, y) \cdot \cos \left[\frac{\pi(2 x+1) u}{2 N}\right] \cdot \cos \left[\frac{\pi(2 y+1) v}{2 N}\right]$

由于该对称函数四个象限的变换结果完全相同,所以

$F_{s}(u, v)=\frac{2}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \cdot \cos \left[\frac{\pi(2 x+1) u}{2 N}\right] \cdot \cos \left[\frac{\pi(2 y+1) v}{2 N}\right]$

把上述变换矩阵定义成归一正交矩阵形式,可得$f_{s}(x, y)$的二维DCT为:

$F(u, v)=\frac{2}{N} K(u) K(v) \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \cdot \cos \left[\frac{\pi(2 x+1) u}{2 N}\right] \cdot \cos \left[\frac{\pi(2 y+1) v}{2 N}\right]$

其中

$K(u)=\left\{\begin{array}{ll} \frac{1}{\sqrt{2}} & u=0 \\ 1 & u=1,2, \cdots, N-1 \end{array}\right.$

$K(v)=\left\{\begin{array}{ll} \frac{1}{\sqrt{2}} & v=0 \\ 1 & v=1,2, \cdots, N-1 \end{array}\right.$

一种更直观的二维正DCT表示形式为:

$F(0,0)=\frac{1}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y)$

$F(u, 0)=\frac{\sqrt{2}}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \cdot \cos \left[\frac{\pi(2 x+1) u}{2 N}\right] \quad(u \neq 0)$

$F(0, v)=\frac{\sqrt{2}}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \cdot \cos \left[\frac{\pi(2 x+1) v}{2 N}\right] \quad(v \neq 0)$

$F(u, v)=\frac{2}{N} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \cdot \cos \left[\frac{\pi(2 x+1) u}{2 N}\right] \cdot \cos \left[\frac{\pi(2 y+1) v}{2 N}\right]$

其中,$f(x,y)$是二维空间向量元素,$F(u,v)$是变换系数矩阵之元素.

二维离散余弦变换的正反变换核是相同的,对称的,可分离的, 即:

$\begin{aligned} Q(x, y, u, v) &=\frac{2}{N} K(u) K(v) \cdot \cos \left[\frac{\pi(2 x+1) u}{2 N}\right] \cdot \cos \left[\frac{\pi(2 y+1) v}{2 N}\right] \\ &=q_{1}(x, u) \cdot q_{2}(y, v)=q_{1}(x, u) \cdot q_{1}(y, v) \end{aligned}$

并记:

$q=q_{1}(x, u)=\sqrt{\frac{2}{N}} K(u) \cdot \cos \left[\frac{\pi(2 x+1) u}{2 N}\right]$

二维DCT的正反变换空间矢量表示形式为:

$F=q \cdot f \cdot q^{T}$

$f=q^{T} F \cdot q$

变换矩阵的形式为:

$q^{T}=\sqrt{\frac{2}{N}}\left[\begin{array}{cccc} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \cdots & \frac{1}{\sqrt{2}} \\ \cos \frac{\pi}{2 N} & \cos \frac{3 \pi}{2 N} & \cdots & \cos \frac{(2 N-1) \pi}{2 N} \\ \vdots & \vdots & \ddots & \vdots \\ \cos \frac{(N-1) \pi}{2 N} & \cos \frac{3(N-1) \pi}{2 N} & \cdots & \cos \frac{(2 N-1)(N-1) \pi}{2 N} \end{array}\right]$

DCT变换的基函数和基图像

如前所述, DCT正变换核反变换可以描述为:

$F(u, v)=\sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y)\cdot Q(x, y, u, v)$

$f(x, y)=\sum_{u=0}^{N-1} \sum_{v=0}^{N-1} F(u, v) \cdot Q(x, y, u, v)$

其中, 正反变换核$Q(x,y,u,v)$也称二维DCT变换的基函数或基图像,$F(u,v)$称为变换系数.

当N=4时的二维DCT变换的基图像如下图:

DCT变换的实例

1
2
3
4
5
6
7
8
9

RGB=imread('C:\Users\Adminstrator');
figure(1);
imshow(RGB);
GRAY=rgb2gray(RGB);
figure(2);
imshow(GRAY);
DCT=dct2(GRAY);
figure(3);
imshow(log(abs(DCT)),[]);