DCT-傅立叶变换的局限性

接上一篇DCT的文章，DCT只提取了整个信号域的频率信息，但是没有时间信息，或者说没有位置信息，也就是说不知道这个频率是发生在哪个位置。
如果在图像处理的语境下的话，那就是不知道这个像素变换发生在什么位置。
我们来写段代码作一个实验：

    img1 = np.zeros((100, 100, 3), dtype="uint8")cv2.rectangle(img1, (60, 20), (70, 70), (255, 255, 255), 1)cv2.rectangle(img1, (80, 20), (90, 70), (255, 255, 255), 1)img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)dft1 = cv2.dft(np.float32(img1), flags=cv2.DFT_COMPLEX_OUTPUT)print(dft1)cv2.imshow("img_ori1", img1)img2 = np.zeros((100, 100, 3), dtype="uint8")cv2.rectangle(img2, (20, 20), (30, 70), (255, 255, 255), 1)cv2.rectangle(img2, (40, 20), (50, 70), (255, 255, 255), 1)img2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)dft2 = cv2.dft(np.float32(img2), flags=cv2.DFT_COMPLEX_OUTPUT)print(dft2)cv2.imshow("img_ori2", img2)

图像显示如下：

两个dft的输出结果一模一样，都是下面的一串：
[[[ 6.1200000e+04 0.0000000e+00]
[ 0.0000000e+00 4.7354695e+04]
[-1.5694755e+04 -1.2207031e-04]
…
[ 1.1026633e-03 -1.1961035e+04]
[-1.5694755e+04 1.2207031e-04]
[ 0.0000000e+00 -4.7354695e+04]]

从上面的这个小实验就可以看出来，两个方框造成的图像像素值变换的频率是一样的，但是出现的位置是不一样的，但是在普通的dft中无法得出（前一篇文章中提到，dft的两个维度，一个是频率，另外一个是像素值，但是没有出现的位置）

STFT 短时傅里叶变换

所谓的“短时”，就是把整个信号域划分成一个一个的小段，就可以得到每一个段的频率和相位谱，再加上每个时间窗这个“时间”维度，就知道了这个“位置-信号中的位置”数据了。当然，还有一个好处就是，可以细分一个时间窗口到某一个稳定信号变换域，就不会受频率变换的影响了。

这里面有任何一个滑动窗口函数涉及到的问题，就是这个窗口到底划多宽合适。同样，写点代码来试试看。

代码是用的python中的scipy框架库里的signal子库：

  scipy.signal.stft（x，fs = 1.0，window =‘hann’，nperseg = 256，noverlap = None，nfft = None，detrend = False，return_oneside = True，boundary =‘zeros’，padded = True，axis = -1 ）

• 输入参数

  • x： STFT变换的时域信号
  • fs： 时域信号的采样频率
  • window： 时域信号分割需要的窗函数，可以自定义窗函数（但是这个方面没有尝试，需要自定义的话请自己尝试）
  • nperseg： 窗函数长度
  • noverlap： 窗函数重叠数，默认为50%。
  • nfft： FFT的长度，默认为nperseg。如大于nperseg会自动进行零填充
  • return_oneside ： True返回复数实部，None返回复数。
    剩下的参数一般不会涉及，采用默认的参数。

• 输出参数

  • f: 频率
  • t： 时间
  • Zxx： STFT时频数据

从另一个角度来理解图像的“时域”数据

我在网上找了很多资料，都没有解释上面几个入参和出参的具体物理含义，也就是这些参数在整个变换过程中起到了什么作用。所以我就先从简单一维数组信号开始作实验，来仔细说说这些参数到底代表什么含义。

看看fs和t这两个参数

一段最基本的代码，在这段代码中，信号就是一个直流信号，没有任何变化。

x = [1, 1, 1, 1]
fs = 1f, t, Zxx = signal.stft(x, fs, nperseg=1)

我们来看一下输出：

首先来看一下这个fs这个入参和t这个出参。

还是把x这个信号看成是一个完整信号的采样结果，代码中的[1, 1, 1, 1]就是对这个信号采样结果，我们这个例子中的 fs = 1，表示这一个数据表示一个信号周期，而t这个出参就是这四个数据是在第几个时间周期内。

如果fs=1，那么上面就是4个周期：1，2，3，4；
如果fs=2，那么上面就是2个周期：1，2。因为每个周期里面有两次采样，也就是采样频率是2；
以此类推

我们改一下参数fs试试看上面的这段描述。

x = [1, 1, 1, 1]
fs = 2f, t, Zxx = signal.stft(x, fs, nperseg=1)

在这个demo中，因为采样频率为2，所以总共才两个周期，4个数据也就是总共4个采样点，分布在2个信号周期内，t就变成[0, 0.5, 1, 1.5]了。

再看看怎么划分窗口

stft对dft的最大的改进就是对信号进行了窗口的划分，对每个窗口内的“短时”信号进行单独的dft变换。
控制窗口的主要由两个参数来控制：
nperseg：表示多少个采样点来组成一个窗口。
noverlap：表示窗口与窗口之间有多少个采样点重合，默认值就是nperseg / 2，这里要做整除。

窗口输出后，同样会在出参t上体现。

x = [1, 1, 1, 1, 1, 1, 1, 1]
fs = 1f, t, Zxx = signal.stft(x, fs, nperseg=4)

我理解，t这个出参实际上是记录了每个窗口的起始周期(时间维度)，那么在这个小demo中，总共分成了4个窗口，共5个轴(axis，5个线形成4个区段，每个轴就是一个时间点)。这是因为窗口与窗口之间必须有2个采样点的重合，所以5个轴的时间点就是0， 2， 4， 6， 8。

为什么说是其实周期呢，我们把fs改成2试一下。

x = [1, 1, 1, 1, 1, 1, 1, 1]
fs = 2f, t, Zxx = signal.stft(x, fs, nperseg=4)

t的输出结果就变成了[0. 1. 2. 3. 4.]，因为2个采样点在一个周期内，所以结果就是原来的一半。

另外，如果分窗的时候不能整除，那么剩下的采样点部分就需要进行补充。
补充的规则由另外两个参数来确定：

scipy官网上的描述如下：

• boundarystr or None, optional
  Specifies whether the input signal is extended at both ends, and how to generate the new values, in order to center the first windowed segment on the first input point. This has the benefit of enabling reconstruction of the first input point when the employed window function starts at zero. Valid options are [‘even’, ‘odd’, ‘constant’, ‘zeros’, None]. Defaults to ‘zeros’, for zero padding extension. I.e. [1, 2, 3, 4] is extended to [0, 1, 2, 3, 4, 0] for nperseg=3.

• paddedbool, optional
  Specifies whether the input signal is zero-padded at the end to make the signal fit exactly into an integer number of window segments, so that all of the signal is included in the output. Defaults to True. Padding occurs after boundary extension, if boundary is not None, and padded is True, as is the default.

另外，补充边缘还需要满足下列的条件。

• In order to enable inversion of an STFT via the inverse STFT in istft, the signal windowing must obey the constraint of “Nonzero OverLap Add” (NOLA), and the input signal must have complete windowing coverage (i.e. (x.shape[axis] - nperseg) % (nperseg-noverlap) == 0). The padded argument may be used to accomplish this.

最后看另外两个出参

以下面这demo为例：

x = [1, 1, 1, 1, 1, 1, 1, 1]
fs = 1f, t, Zxx = signal.stft(x, fs, nperseg=4)

先看下t，这个数组上面说到了，这就是多少个时间窗，这个例子中就是：
[0. 2. 4. 6. 8.]
表示总共有5个时间窗。

先看看频率，在signal的库中，f是一个数组，我自己的理解，这里的频率是指一个时间周期内的样本梯度。

输出的f为：
[0. 0.25 0.5 ]，这表示总共有三种频率。我们拆开了看这个过程：
nperseg = 4的话，那么overlap就等于2，上面看到了总共是5个时间窗，那么我理解补充完边界的结构是：
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]
时间窗就是：
[0, 0, 1, 1]
[1, 1, 1, 1]
[1, 1, 1, 1]
[1, 1, 1, 1]
[1, 0, 0, 0]

上面几个时间窗中的样本梯度分布就是：0.5，0.25，0

Zxx返回结构

Zxx结构就是一个二维数组：
[[ 0.75+0.j 1. +0.j 1. +0.j 1. +0.j 0.25+0.j ]
[-0.5 +0.25j -0.5 +0.j -0.5 +0.j -0.5 +0.j 0. -0.25j]
[ 0.25+0.j 0. +0.j 0. +0.j 0. +0.j -0.25+0.j ]]

横着的就是时间窗，竖着的就是梯度分布。

plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud')plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

使用上面的代码就可以画出一个时域-频域的分布图。每个交叉点上就是Zxx二维数组中的频率+相位的值。

图像的stft

图像是个二维数据，而上面的过程的输入都是一个一维数组，二维数组的是stft是对每行数据作一次stft，然后保存在zxx数据结构中。

# 分别作stft
f1, t1, Zxx1 = signal.stft(img1, fs=1, nperseg=100)
f2, t2, Zxx2 = signal.stft(img2, fs=1, nperseg=100)# 第二十行有区别
plt.pcolormesh(t1, f1, np.abs(Zxx2[20]), vmin=0, vmax=200, shading='gouraud')plt.pcolormesh(t2, f2, np.abs(Zxx2[20]), vmin=0, vmax=200, shading='gouraud')

结果如下：

从分布上就能看出来，stft可以按照窗口来区分不同的频率，从而达到区分两者的目的。