Python频率函数 python频域分析

文章目录

• 什么是域

• 空间域
• 频域

• 空间域到频域的转换

• DCT(离散余弦转换)

• DCT的基本逻辑：
• 简单一维数组的DCT DEMO
• 换一组频率变化大和变化小的一维数组

• DCT逆变换
• 二维DCT
• 真实图像
• DFT，离散傅里叶变换

• 实例

• DFT

在各种数字图像处理的资料中，频繁出现如下的几个概念：

• 空间域
• 频域
• 时域
• 傅里叶变换 & 傅里叶级数
• DCT
• DFT

特别是几个概念搅和在一起，分分钟劝退。琢磨了半天，在在这个系列中插播一点理论方面的理解。还是那句老话，我是用比较通俗的话说出我自己对这些概念的理解，不一定很严谨，但是基本正确，能更好的理解。毕竟我们是程序员，不是专门做数学研究的。

什么是域

域是一个数学上的概念。还记得数学中的定义域和值域么？就是表示变量的数字在哪个范围里面。在数学中，一般表示某种映射关系的话，都可以使用X-Y坐标系来表示，X轴(横轴)表示自变量，Y轴(数轴)表示因变量。

表示某种数字信号或者函数与其因变量之间的映射关系。

空间域

空间域就是说，自变量是图像中的像素点位置，而值域就是针对这个位置所标记的像素的值进行一些操作和变换。这个就是空间域，这些操作也就叫做空间域操作。
$y = f(x,y)$

其中的x, y就表示图像中像素的x坐标和y坐标，y就表示这个像素的像素值，f表示这样一种映射关系。因为这个映射关系的自变量是x和y坐标，这是一个空间上的概念，所以这样一种映射关系就可以称作图像的空间域表达。

频域

如果把图像中的一行抽取出来，比如下面这张图(512 * 512的尺寸)：

图一

我们提取出其中的第100行数据：

img = cv2.imread("e:\\lena.jpg", 0)
line = img[100]
print(line)

输出的结果是这样的一个数组：

[ 93 88 85 85 85 84 88 92 92 92 92 92 92 92 92 92 92 92 …]

用plt画出来的曲线是这样的：

图二

这个图像实际上可以描述这一行图像像素数据的灰度变换程度，某种程度上就可以被成为图像的频率，针对图像这样频率的处理的方式就可以称为频率域。

空间域到频域的转换

DCT(离散余弦转换)

一般的空间域到频域的转换都是从DFT，也就是离散傅里叶变换开始说起，但是我觉得DCT更好理解一点，所以先说说DCT。

DCT的基本逻辑：

1. 和频率域的基本逻辑一样，把图像中的某一行拿出来，其中的每一个像素点可以被描述成(x, y) -> 像素值的映射关系，其中x还是固定的。如果用像素的y做横坐标，像素值做纵坐标，可以形成一个二维坐标系。这些像素数据或者形成的二维曲线可以被成为原始信号。

1. 根据傅里叶变换的概念，任何函数曲线都可以用正弦或者余弦去模拟，那么上面这个图二肯定也是可以的。
2. 傅里叶变化是可以用很多个cos或者sin函数去逼近的，那么这里用多少个呢。因为上面的曲线是由若干个点组成的，可以理解是对这个未知曲线做的一个采样。
3. 而dct中的逻辑是使用采样点个cos余弦函数去逼近，也就是这个函数被认为是

$f(x) = X_0cos(0)+X_1cos(x)+X_2cos(2x)+...+X_Ncos(Nx)$

要想确定这个函数，就是要找出这N个权重系数X，用N个余弦波是为了后面矩阵计算的方便与可行。

1. 关于上面这一条，这里我看了一个视频，我对这个地方的理解是，第k个样本点可以替代频率为kx的余弦函数在这个函数中占的比例，总共N个样本点，所以是N个余弦函数。

放上视频的链接：

1. 使用下面的dct公式，来通过样本点计算这些权重系数X。

$X_k = \sum_{n=0}^{N-1}x_ncos[\frac{(2n+1)\pi k}{2N}]$

公式一

1. $x_n$可以称作原始信号采样点，而cos可以被看做是这个样本点代表余弦频率的余弦波的采样点。如果把上面的公式一写成向量形式：

$X=Cx$

实际上就是表示第k个样本点所代表的kx的余弦波的N个采样点的值向量C，然后由N个这样的行向量组成了矩阵C，这个矩阵与原始信号向量x进行点乘，形成最后的权重向量X。

可以参考下面两张图。

简单一维数组的DCT DEMO

假设有一个一维数组：(1, 1, 1, 1, 1)

那么根据dct的公式：

因为k=0，所以代入dct公式中，cos部分都等于1了，所以结果就是5。
$X_0=1 + 1 + 1 + 1 + 1 = 5$
因为所有的$x_n$都是1，所以公式就只剩下余弦部分了。
$X_1 = cos[\frac{(2*0+1)*1\pi}{10}] + cos[\frac{(2*1+1)*1\pi}{10}] \\ + cos[\frac{(2*2+1)*1\pi}{10}] + cos[\frac{(2*3+1)*1\pi}{10}] + cos[\frac{(2*4+1)*1\pi}{10}] \\ = cos(\frac{\pi}{10}) + cos(\frac{3\pi}{10}) + cos(\frac{\pi}{2}) + cos(\frac{7\pi}{10}) + cos(\frac{9\pi}{10}) \\ =0$
其他三个X大家有兴趣也可以算一算，都是等于0。所以输出就是（5,0,0,0,0），当然opencv输出的shape是一个(5, 0)的数组，而且opencv在公式一的基础上还增加了一点东西，矩阵$C$加了一个前缀：
$C_{ij}^{(N)}=\sqrt{\alpha_j/N}cos(\frac{\pi (2k+1)j}{2N}) \\ \alpha_0=1, \alpha_j=2, j>0$
所以输出值不是5，没有去细究了，大体的逻辑肯定是一样的。

line = (1, 1, 1, 1, 1)
dft = cv2.dft(np.float32(line))

print(dft.shape)
print(dft)

输出的dft如下：

• 物理含义，假设(1, 1, 1, 1, 1)是像素值，那么这张图片的像素之间是没有任何变化的，也可以理解没有频率，和直流电一样，没有变化，最终变换到频率域之后，就只有常量$X_0$上有作用，很好的表示了图像在频率上的描述。
• 而且常量值可以很好的表达像素值，如果像素值比较大，这个常量值$X_0$也会比较大。
• 在上面的向量公式中，矩阵C中的不同行向量之间是互相正交的，也就是点积为0。在这个小例子中，原始信号就是直流信号，只有$X_0$会对其影响，所以其他的行向量进行点乘后，结果都为0。

换一组频率变化大和变化小的一维数组

换一个一维数组：(1, 100, 200, 2, 4, 250, 23, 198, 2, 251)。

这个数组的输出值为：

再看看变化小的一个数组：(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)

从上面两个数组的dct输出可以看出，频率变化大和变化小的像素信号经过DCT转换后由很大的区别，那么就很好对这样的一些信号去进行区分和处理了。

DCT逆变换

根据上面提到的DCT的向量表达信息，实际上只需要求出矩阵$C$的逆矩阵$C^{-1}$，再用这个矩阵去乘以权重向量X就可以得到原来的数组了。
$x = C^{-1}X$

demo_line=(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
dct_demo = cv2.dct(np.float32(demo_line))
iLine = cv2.idct(dct_demo)
print(iLine)

就可以得到原来的数据了。

二维DCT

从一维往二维拓展一下，比如一个8 * 8的矩阵，先对每一行做一次DCT操作，最终可以得到一个 8 * 8的DCT矩阵，然后再针对每一列做一次DCT，还是获得一个8 * 8的DCT矩阵。如下图所示。

真实图像

我针对上面的图一用opencv的dct函数做了一次转换：

dct_img = cv2.imread("e:\\lena.jpg", 0)
dct = cv2.dct(np.float32(dct_img), flags=cv2.DFT_COMPLEX_OUTPUT)

cv2.imwrite("e:\\aaa.jpg", dct)

打开图片一看，是这么个东西。

我的理解如下：回忆一下二维DCT的过程，在X方向上，从左往右的数据分别是
$cos(0x), cos(x), cos(2x), cos(3x) ...... cos(Nx)$
的权重系数
$X_0,X_1,X_2......X_N$
而在图像中，只有左上角是比较黑的，所以可以说这幅图像在X方向上的频率(也就是像素值的变化率)集中在低频区。换言之，在Y轴方向上的频率也是集中在低频区。这样就可以采用一些低通，高通等不同的滤波算法进行区分和操作了。

当然，我个人觉得这里的频率变换可能也是相对的，我感觉最右边的可能也到了$cos(100x)$的高频了吧。

DFT，离散傅里叶变换

在图像处理中，用的更多的是这个DFT，也就是离散傅里叶变换。应该说DCT是DFT的一种特例，DCT说完了，就再说说DFT。

从DFT的公式说起：
$X_k = \sum_{n=0}^{N-1}x_ne^{-2\pi j\frac{kn}{N}}$

然后根据欧拉公式：
$e^{jw}=cos(w)-jsin(w)$
所以DFT公式可以展开为：
$X_k = \sum_{n=0}^{N-1}x_n[cos(2\pi kn/N) - jsin(2\pi kn/N)]$
相对于DCT，DFT有两个部分，实数部分与虚数部分，其中实数部分是和DCT一样的。那么虚数部分代表啥呢，图像中的像素是一个一个真实的像素点，而且是0-255之间的一个整数。

其实我的理解是在opencv中，傅里叶变换就是把这两个部分做了图像两个维度的变换，一个叫频率谱(实数部分)，一个叫相位谱(虚数部分)。不要把它和数学中的虚数一一对应，就是两个维度的数字就好。

其中频率谱部分在DCT的部分已经说过了，就是代表了图像两个方向的像素变换情况。

那么相位谱又代表了什么呢？我在网上找了蛮久的资料，没有找到准确合适的描述，下面这段话是我自己的理解，去理解图像的相位谱到底代表啥意思。

• 首先，图像的数据是离散的，在信号处理领域，如果信号上的点是离散的，相当于对信号进行采样：

图中有8个采样点，每个点的相位就是他在X轴上的位置，比如第一个点的相位就是$\pi /8$。

• 那么回到图像中，一行图像的数据可以被认为是某种信号，这个上面说过了。

那么这个信号可以拆成若干个频率不同的余弦波。那么我的理解是

傅里叶的虚数部分，也就是相位谱，就是这个采样点在频率为x或者2x的不同频率中的波信号中所对应的相位。

比如第一个点，可能在 cos(x)这个信号中的相位是$\pi /9$，cos(2x)这个频率的信号中的相位是$\pi /4$，类似这种的。

• 那么在图像中的物理意义，我是这么理解，频率谱是表达了像素之间的变化，而相位包含了在某个频率信号中的位置，这个位置可以得到这个信号中具体的取值，可以对应到具体的信号幅度，某种意义上来说可以代表在图像中的像素值。

实例

DFT

python代码：

import numpy as np 
import cv2 as cv
from matplotlib import pyplot as plt

# 1 读取图像

img = cv.imread('e:\\lena.jpg',0)

# 2 傅里叶变换

# 2.1 正变换

dft = cv.dft(np.float32(img),flags = cv.DFT_COMPLEX_OUTPUT) 

# 2.2 频谱中心化

dft_shift = np.fft.fftshift(dft)

# 2.3 计算频谱和相位谱

mag, angle = cv.cartToPolar(dft_shift[:,:,0], dft_shift[:,:,1], angleInDegrees=True)
mag=20*np.log(mag)

# 3 傅里叶逆变换

# 3.1 反变换

img_back = cv.idft(dft)

# 3.2 计算灰度值

img_back = cv.magnitude(img_back[:,:,0],img_back[:,:,1])

# 4 图像显示

plt.figure(figsize=(10,8))
plt.subplot(221),plt.imshow(img, cmap = 'gray') 
plt.title('输入图像'), plt.xticks([]), plt.yticks([]) 
plt.subplot(222),plt.imshow(mag, cmap = 'gray')
plt.title('频谱'), plt.xticks([]), plt.yticks([]) 
plt.subplot(223),plt.imshow(angle, cmap = 'gray')
plt.title('相位谱'), plt.xticks([]), plt.yticks([])
plt.subplot(224),plt.imshow(img_back, cmap = 'gray')
plt.title('逆变换结果'), plt.xticks([]), plt.yticks([])
plt.show()

输出图像为：