傅里叶变换->小波变化

傅里叶变换->小波变化

声明:文中大多数内容来自(知乎1335)[https://www.zhihu.com/people/zhi-yuan-ya-77],matlab源代码:1368069096/From_FT_to_WT_examples-,部分为个人理解阐明

傅里叶变换FT

基础知识

(FOURIER TRANSFORM,简称FT)

为什么傅里叶变换可以把一个信号从时域变换到频域?

先给出公式,傅里叶变换的形式为:\(X(w)=\int_{-\infty}^{+\infty} x(t) e^{-j w t} d t\)

PS:傅里叶变换还存在系数,有的文章写的是 \(\frac{1}{2 \pi}\) ,有的文章写的是\(\sqrt\frac{1}{2 \pi}\) ,两个系数只要满足正变换系数乘上逆变换系数等于\(\frac{1}{2 \pi}\) 即可。这是为了保证经过一次正变换和反变换之后,得到的信号与原信号幅值相同,与我们接下来的讨论关系不大。

1.理解变换公式

我们知道,根据欧拉公式,\(e^{-j w t}=\cos (w t)-j \sin (w t)\)。也就是说,傅里叶变换的本质就是:将原始信号乘上一组三角函数(正余弦),之后在整个时间域上积分。就这么简单!

y=sin(3t)图

我们来看一个信号:y = sin(3t),如下图:

img

很好的周期性质,且每个周期的积分值都是0。如果对这个函数在\((-\infty,+\infty)\)积分,那就是基本是0,因为 \((-\infty,+\infty)\)包含了无数个周期。

PS:虽然这个积分在高数上不可积,但是你应该明白这里我要表达的意思:因为良好的周期性,且每个周期积分值是0,那么最后在很长的一段时间区间上积分,得到的还是一个很小的数,近似为0。

我们来用一段较长的时间区间计算一下,\(\int_{0}^{50} \sin (3 t) d t=0.1002\),结果符合我们的预计。

y=sin(4t)sin(3t)图

现在,我们来将这个信号乘上一个sin(4t) ,则信号变为y1 = sin(3t)*sin(4t),如下图:

img

具有一个较短的小幅震动的周期和一个较长的主体周期,对吧?且每个主体周期的积分值都是0。同以上讨论,如果对这个函数在\((-\infty,+\infty)\)积分,基本还接近于0,因为 \((-\infty,+\infty)\)包含了无数个主体周期。

y=sin(3.1t)sin(3t)图

之后呢,我们来将这个信号乘上一个sin(3.1t) ,则信号变为y2 = sin(3t)*sin(3.1t),如下图:

img

同样是有一个较短的小幅震动的周期和一个较长的主体周期,对吧?可以判断,每个主体周期的积分值都是0(虽然(0,50)这个区间没有完整地展示主题周期)。那么,依然同以上讨论,如果对这个函数在\((-\infty,+\infty)\) 积分,基本还接近于0,因为\((-\infty,+\infty)\)包含了无数个主体周期。

我们来用之前的时间区间计算一下, \(\int_{0}^{50} \sin (3 t) * \sin (3.1 t) d t=-4.7731\)

咦?这一次怎么距离0这么远了呢?

原因就是:对于sin(3t)*sin(4t),它的主体周期较短,(0,50)是包含了好几个主体周期的,也就是说(0,50)在某种程度上是类似于 \((-\infty,+\infty)\) 的。但是,对于sin(3t)*sin(3.1t),它的主体周期很长,(0,50)连它的一个完整的主体周期都没有包含,那么(0,50)是不能类似于 \((-\infty,+\infty)\) 的,积分值自然就比较大。

我们此时可以这样小小总结一下,对于信号y = sin(3t),它的频率是3rad/s,(如果你喜欢用HZ,那就除以\(2\pi\) ,就是 \(\frac{3}{2\pi}\)HZ,这里使用rad/s,是为了与前面的傅里叶变换的公式中的w一致),而sin(4t)的频率是4rad/s,sin(3.1t)的频率是3.1rad/s。

如果在 \((-\infty,+\infty)\) 积分,那么y1 = sin(3t)sin(4t),y2 = sin(3t)sin(3.1t)的积分值都是0,也就是说,sin(4t)和sin(3.1t)在这里是没差别的。

但是!!!如果在一个有限区间内积分,由于sin(3.1t)的频率3.1rad/s,距离原信号y = sin(3t)的频率3rad/s更近,那么sin(3.1t)和sin(3t)的乘积,也就是y2 = sin(3t)*sin(3.1t)的积分的绝对值会更大,也就是会离0更远。这里已经显示出一定的频率选择性了。

最后,让我们请出我们今天的主角,将这个信号乘上一个自己同频率的sin(3t) ,则信号变为y3 = sin(3t)*sin(3t),如图:

img

Amazing!!!发现了什么?良好的周期性?还有呢?由于乘上了自己,任何时间的幅值都大于等于0了!不再满足周期内积分值为0这个条件了!那么,此时,我们对这个信号在 \((-\infty,+\infty)\) 积分,就会得到一个非常非常大的数字。这个很大很大的数字就告诉你,这个信号和你乘的信号是同频率的!这就是可以知道信号中具有哪些频率部分了,不是吗?

我们还是来用之前的时间区间计算一下, \(\int_{0}^{50} \sin (3 t) * \sin (3 t) d t=25.0833\)是不是比其他的积分值都大了好多?

好了,我们已经知道,▲.将一个信号乘上一个特定频率的sin函数,在 \((-\infty,+\infty)\) 上积分,可以将信号中与sin函数同频率的部分筛选出来。那么,原则上讲,只要乘上所有频率的sin函数,并积分,不就知道原始信号中的所有频率部分了吗?

但是这样做需要把所有频率乘进去,做无数次计算哈!算不出来的。所以,我们将所乘的sin函数的频率作为符号变量w,来进行积分,即:

\(X(w)=\int_{-\infty}^{+\infty} x(t) \sin (w t) d t\)

注意:这里的w只是一个符号变量,这样的话,就只需要做一次积分,可以计算了。

计算出来X(w)之后,想知道特定的频率w0对应的积分值,直接将w0带入X(w)就立马得到积分值。,如想知道是否含有w0=3rad/s分量,那么久计算X(3)看结果是否为0,这样就能知道原信号中是否含有这一频率的部分了。

好了,我们推导的这个式子,是不是与傅里叶变换的式子:

\(X(w)=\int_{-\infty}^{+\infty} x(t) e^{-j w t} d t=\int_{-\infty}^{+\infty} x(t)(\cos (w t)-j \sin (w t)) d t\)

很像了呢?

这就是傅里叶变换的原理!乘上带有符号变量的sin、cos函数,并积分,就知道原始信号中的所有频率部分啦!

2、傅里叶变换(FT)的正交性

傅里叶变换是一种变换。在变换中,我们将原始信号乘上的变化信号称为基函数

在傅里叶变换中,一系列不同频率的sin、cos等函数称为这个变换的基函数。至于为什么需要既使用sin,又使用cos,这涉及到一点点正交函数的概念。傅里叶变换中的不同频率的sin、cos等函数是正交函数,使用正交函数组成的基函数会带给变换一些方便。

什么是正交性?

向量正交:

我们都知道,向量 [公式] 的内积为[公式]正交的定义为内积为0,即\(a \cdot b=0\)。如\(a=(1,0)\) 表示x轴,$b=(0,1) $表示y轴,则 \(a \cdot b=0\)即意味着x轴与y轴正交。

假设有一个向量 \(v=\left(x_{0}, y_{0}\right)\)\(x_0\) 在x轴上定位v,\(y_0\) 在y轴上定位v。当x轴与y轴正交时,意味着x坐标和y坐标表示的信息是彼此独立的,两坐标可以完全定位v。

那么,当我们已知向量v,已知x轴a与y轴b,如果知道v的坐标呢?答案就是,投影/内积

将v向x轴a做投影/内积: [公式] ,可以得到[公式];将v向y轴b做投影/内积: [公式],可以得到[公式]

正交基

类似的,函数$ f_{1}(x), f_{2}(x) \(的内积定义如下:\)\int_{-\infty}^{+\infty} f_{1}(x) f_{2}(x) d x$。函数正交的定义同样为内积为0,即 \(\int_{-\infty}^{+\infty} f_{1}(x) f_{2}(x) d x=0\) 。PS:对于周期函数,定义中的积分区间为一个周期T。

我们来看,

\(\int_{0}^{T} \sin \left(w_{1} x\right) \cos \left(w_{2} x\right) d x=0 (1)\)

\(\int_{0}^{T} \sin \left(w_{1} x\right) \sin \left(w_{2} x\right) d x=0, w_{1} \neq w_{2}(2)\)

\(\int_{0}^{T} \cos \left(w_{1} x\right) \cos \left(w_{2} x\right) d x=0, w_{1} \neq w_{2}(3)\)

▲.因此,傅里叶变换中的不同频率的sin、cos等函数都是正交函数。我们将cos想象为一个轴,将sin想象为一个轴 ,这两个轴,就张了一个函数空间

我们已经知道,x(t)的傅里叶变换为:

\(X(w)=\int_{-\infty}^{+\infty} x(t)(\cos (w t)-j \sin (w t)) d t=\int_{-\infty}^{+\infty} x(t) \cos (w t) d t-j \int_{-\infty}^{+\infty} x(t) \sin (w t) d t\)

我们来看,实部 \(\int_{-\infty}^{+\infty} x(t) \cos (w t) d t\) 即相当于将x(t)与cos做内积,即向cos轴投影,得到的是在这个函数空间内的cos坐标,也就是与cos的相似度;虚部\(\int_{-\infty}^{+\infty} x(t) \sin (w t) d t\) 即相当于将x(t)与sin做内积,即向sin轴投影,得到的是在这个函数空间内的sin坐标,也就是与sin的相似度。cos坐标和sin坐标都确定后,该信号就确定了。

因此,傅里叶变换之后,实部是与cos的相似度,虚部是与sin的相似度,傅里叶变换也就是与cos的相似度和与sin的相似度的总和,也就表示了相应的频率信息。

相似度的理解:通过相似度,可以将跟本身频率相似的频率挑选出来

3、小例子

最后,我们对于y=sin(3t)做傅里叶变换(这里画图用的matlab的FFT,是FT的离散快速算法,不是积分出来的,但是原理是相同的,仅做展示使用),变换后的图像如下:

img

可以看到,与我们的预期相同,变换之后,在 \(\omega=3 rad/ s\) 处,出现了峰值,即表示该信号中包含了 的 \(\omega=3 rad/ s\) 频率成分。

这里需要说明三点:

1、这边的作图结果不是理想FT,如果是理想的FT,即在 \((-\infty,+\infty)\) 上积分,那么除了 \(\omega=3 rad/ s\) 的积分值将是 \(+\infty\) 之外,其他任意频率处的值都应该是0,得到的将是一个冲激函数。但是,这里我用的是离散傅里叶变换(详见我的下一篇文章),可以暂时理解为类似于我们讨论过的有限区间,当频率靠近3rad/s时(如之前提到的3.1rad/s的例子),也会积分出来一个较大的数值,所以这里不是一个冲击函数,而是一个山峰状的函数;

2、傅里叶变换之后是具有实部和虚部的,实部是与cos的相似度,虚部是与sin的相似度。我们要频率信息的时候,不管是与某一频率w的cos的相似还是与某一频率w的sin的相似,都是这一频率w嘛,不需要区分。因此,这里画图取了模,就不存在实部和虚部了。

3、在\(\omega=3 r a d / s\)处之所以没有出现我们讨论的很大很大很大的值,是因为画图之前对于变换之后的幅值统一除以了信号取样点的个数,统一压缩了一定的倍数。

总结

1.由于sin和cos是一对正交基,所以所有的信号都可以分解成正弦信号与余弦信号的线性叠加。

2.除了与乘以本身频率w的正余弦相同信号再在\((-\infty,+\infty)\)范围上积分的结果会是非0以外,乘以其他频率信号的结果都是0.通过这个特性可以将是否含有w频率信号鉴别出

2.傅里叶变化的基函数\(e^{jwt}\)由欧拉公式可分解为\(e^{-j w t}=\cos (w t)-j \sin (w t)\),所以根据公式:\(X(w)=\int_{-\infty}^{+\infty} x(t) \cos (w t) d t-j \int_{-\infty}^{+\infty} x(t) \sin (w t) d t\),前部分乘以cos(wt)能挑出cos成分中的信号频率,后部分乘以的sin(wt)能挑出sin成分中的信号频率。再将挑出的各信号频率线性叠加的结果就是分解信号总的含有各谐波分量的频率


傅里叶变换(FT)的缺点与短时傅里叶变换(STFT)

离散傅里叶变换(DFT)

在本文正式开始之前,我们需要明确一下实际信号进行的FT的一些特殊之处。实际采集的信号往往是这样的:

img

实际的信号往往具有两个特点:

1、离散性,就是采集数据不连续,很容易理解,采集信号肯定是一个一个数据采集的;

2、有限性,虽然理想的傅里叶变换是从 \((-\infty,+\infty)\) 进行积分,但是实际信号往往实在一个区间内(a,b)的。

所以,我们要用到离散傅里叶变换(DISCRETE FOURIER TRANSFORM,简称DFT),DFT与FT相比,就是多了两个特征:1、离散型,2、有限性

我们来一起试一试如何推导DFT公式。首先设采集了N个信号点,其时刻为 \(t_{0}, t_{1}, \ldots t_{N-1}\),对应时刻采集到的信号值为\(x\left(t_{0}\right), x\left(t_{1}\right), \ldots x\left(t_{N-1}\right)\) 。很自然的,原来信号连续,是积分,现在数据离散了,那就是把积分变成累加。于是我们得到: \(X(w)=\sum_{i=0}^{N-1} x\left(t_{i}\right) e^{-j w t_{i}}\)

这么一来,我们发现,原信号有N个数据点,DFT变换后的信号却变成连续的了,我们将之称为离散时间傅里叶变换(DISCRETE TIME FOURIER TRANSFORM,简称DTFT)

DTFT有两个缺点,第一, [公式] 且连续,需要进行无数次计算,计算机无法计算;第二,在进行计算的时候,我们需要已知:\(t_{0}, t_{1}, \ldots t_{N-1}\)\(x\left(t_{0}\right), x\left(t_{1}\right), \ldots x\left(t_{N-1}\right)\) ,但是调用过FFT函数的同学都知道,FFT只需要已知\(x\left(t_{0}\right), x\left(t_{1}\right), \ldots x\left(t_{N-1}\right)\)就可以进行。

这是怎么回事呢?

首先,我们使用相对采样时间 [公式] 代替真实采样时间 [公式]可以得到: [公式]

此时我们发现, \(X(w)\) 变成了以 \(2\pi\)为周期的函数,即 [公式]

那么,我们只需要计算 \(w \in(0,2 \pi)\) 区间的 \(X(w)\) ,就可以得到 \(w \in(-\infty,+\infty)\) 区间的\(X(w)\)了。也就是说,通过使用相对采样时间 [公式] 代替真实采样时间\(t_{0}, t_{1}, \ldots t_{N-1}\),我们将 [公式] 的范围从 \((-\infty,+\infty)\) 缩小到了(0,\(2\pi\)).

自然地,我们希望将(0,\(2\pi\))离散化称为N个点,这样计算机就可以计算了!则取\(w=\frac{2 \pi k}{N}, k=0,1, \ldots N-1\),则有: \(X(k)=\sum_{n=0}^{N-1} x(n) e^{-j \frac{2 \pi k}{N} n}, k=0,1 \ldots N-1\)

好了,这就是离散傅里叶变换DFT了!

接下来,来看DFT的两个性质:

  • 第一, [公式] ,即 [公式] 是所有元素的和,通常会比其他的元素大几个数量级。

  • 第二, [公式][公式][公式] ,即第二个元素和最后一个元素共轭,同理有[公式] 等等。

如下图所示,DFT之后的N 个元素中,第一个为均值;之后的 N-1个元素,只有一半元素是独立的。

img

需要说明,这里 [公式] 是一种相对频率,独立元素中,最小相对频域为1,最大相对频率为 [公式]要想把 [公式] 还原到真实的频率 [公式] ,只需要 [公式] ,将 [公式] 映射到 [公式] 即可[公式] 为采样频率。

PS:简单说一下,根据香农采样定理,当采样频率为 [公式] 时,能采到的最大信号频率为[公式]。因此,将相对频率 [公式] 通过公式 [公式] ,得到的 [公式] 就在区间[公式]内,也就是真实频率的区间了。

所以DFT公式为: [公式]

PS:DFT公式的形式很多,有的是从时域到频域,有的是从空间域到频域,但是本质都是一样的,抓住离散性有限性两个特点即可。离散性是指积分变成了累加,有限性是指积分/累加区间不是 \((-\infty,+\infty)\) 了,而是一个有限区间了。

傅里叶变换(FT)的缺点

应该说明,虽然本小节的题目是FT的缺点,但是FT和DFT的本质是相同的。由于信号都是有限长度的、离散的,所以接下来进行的都是DFT,不过在某些部分为了方便理解,还是写了FT的公式。在看本文的时候,你不需要刻意区分这两个概念。

我们现在来看两个信号,如下图:

\(y 1=\sin (5 t) *(0<t<25)+\sin (t) *(25<t<50)\)
\(y 2=\sin (t) *(0<t<25)+\sin (5 t) *(25<t<50)\)

img

img

这两个信号都是由sin(t)和sin(5t)组成的,y1是先出现了sin(5t),再出现了sin(t),y2是先出现了sin(t),再出现了sin(5t)。

我们对它们进行FT,看看他们包含怎么样的频率,如下图:

img y1,FT

img y2,FT

Amazing!发现了什么?变换后的结果是一模一样的,都在w=1rad/s和w=5rad/s出现了峰值!这就可以说明FT的缺点了——FT只能提供频域信息,而完全丢失了时域信息!!!

不管某一频率的信号出现的时间是早还是晚,FT都是将它一视同仁地乘上sin和cos(FT的变换基函数),然后在整个时间区间加和。因此,它不能提供某一频率信号出现的时间。

比如,对于上面两个信号,FT只能告诉我们,它们都有1rad/s和5rad/s的频率,而不能告诉我们1rad/s和5rad/s分别出现在哪个时间段。

所以,怎么办呢???

那就是把信号分成左右两半啊!左边进行一次FT,右边进行一次FT,很简单吧!好了,这就是短时傅里叶变换(STFT)的基本原理。

所以,接下来我们要正式开始步入——短时傅里叶变换(STFT),看看它是如何解决这个问题的。

短时傅里叶变换(STFT)

如上所述,我们将信号从中间截断,左边进行一次FT,右边进行一次FT,分别来看看。

imgy1左

imgy1右

imgy2左

imgy2右

可以看出,y1的左半部分是5rad/s,右半部分是1rad/s,y2恰好相反。这就说明,在y1中,(0, 25)的信号是5rad/s的频率,(25, 50)的信号是5rad/s的频率,y2恰好相反。这就是短时傅里叶变换的基本原理。

但是数学嘛,能用一个公式表达的,就别用一段话表达,截断、切开这些语句太不专业了。截断、切开的操作,更专业的讲叫作分窗,其实是可以通过数学上的处理变成DFT变换的基函数的一部分的。接下来我们来看一看。

首先,你可以想象一下,有一个窗子在这个信号上从左向右滑动,每次你都只能看到这个信号的一部分,所以我们把这个长度叫作窗长width。

现在我们来定义一个方窗函数\(y_{\text {window}}=1 *(-\text {width} / 2<t<\text {width} / 2)\) ,如下图,即是width = 10 的一个方窗函数:

img

定义了方窗函数之后,我们只需要对方窗函数进行平移,再与原信号作乘,就相当于原来的截断、切开的操作,因此这种操作更专业地叫作分窗。

那么,将方窗函数向右平移了 [公式](s可能是sliding的意思吧),再与原信号相乘,由于方窗函数除了中心的width部分是1外,其他部分都是0,这就相当于提取出了原信号在[公式]处,宽度为width的部分,这个信号分窗这个操作就可以写成: [公式]

如下两图所示,将 [公式][公式] 相乘,就相当于取出来了 [公式] 中的(20,30)中的一段。

那么,我们对原信号中被提取出来的这一部分进行FT,就可以写成: \(X\left(w, t_{s}\right)=\int_{-\infty}^{+\infty} y_{w i n d o w}\left(t-t_{s}\right) y(t) e^{-j w\left(t-t_{s}\right)} d t\)

PS:这里之所以 [公式] 要变成 [公式] ,是为了保证做FT的时候相乘的基函数具有统一性。

如此,变换后的[公式]代表原信号在[公式]处、宽度为width的部分的傅里叶变换,也就可以提取出来原信号在[公式]处、宽度为width的部分,包含各个频率部分的多少!带入不同的 [公式] ,也就是随着窗子的滑动,就可以知道不同的时间段内频率的成分。

我们采用width为10的方窗函数对\(y 3=\sin (20 t) *(0<t<25)+\sin (t) *(25<t<50)\)进行STFT,如下:

首先,方窗函数位于 [公式] 处,与原始信号相乘,选择出(0,10)的信号。

img

对选择出来的信号进行FT。可以看到,当t=5s时,选择的时间区间为(0,10),这一部分只包含了 [公式] 的频率成分。

img

之后,方窗函数向右移动,与原始信号相乘,选择出不同时间区间的信号,进行FT。这里选择t=25s进行展示。

img

img

可以看到,当t=25s时,选择的时间区间为(20,30),这一部分即包含了 [公式] 的频率成分,也包含了 [公式] 的频率成分。

重复以上过程,我们可以将方窗函数选择的不同时间区间的信号的FT的结果拼合起来,形成一张三维图。由此,我们即可知道,[公式] 的时间区间内,信号具有怎么样的频率成分。

img

通过width = 10的方窗的STFT结果,我们可以知道,对于信号:[公式] ,在(0,10)、(10,20)时间区间内,具有20rad/s的频率成分;在(20,30)时间区间内,具有1rad/s和20rad/s的频率成分;在(30,40)、(40,50)时间区间内,具有1rad/s的频率成分。

最后,进行三点重要的讨论。

第一点,变换之后的 [公式] 是一个三维函数,它有两个自变量, [公式] 和w。[公式] 指的是原信号在[公式]处,w上一篇文章我们已经讨论过了,就是频率。所以,STFT提取出来的信息就是:原信号在[公式]处、宽度为width的部分,包含的频率信息。

原则上讲,可以得到任一[公式]对应的频率成分,如下图。

img

但是 [公式] 是连续的,并不意味着你知道了每个时刻的频率成分,你知道的还只是 [公式]这一段区间内的频率信息。所以一般不需要计算所有的 [公式],每隔width计算一次即可。

你或许会想,我把width缩小一些,不就可以知道更精确的时间范围内的频率了吗?是的,你的猜想很对!但是,如此做也会带来一些频域分辨率的问题。这一点涉及到一些时域分辨率和频域分辨率的知识,我们下一篇文章会着重讲。

本质变化:

第二点,方窗函数是可以包含入变换基函数内部的,这组成了新的基函数,同时反映了STFT的本质。

我们来看, 如果定义 [公式] ,那么[公式]

那么,STFT的公式: [公式] 就可以写成: [公式]

我们在上一篇文章里说过,变换就是将原信号乘上一个基函数,再积分的过程,那么,SDFT的基函数就是 [公式]

Amazing!所以,STFT的本质是什么呢?

STFT的本质就是将FT的基函数 [公式] 乘上一个方窗函数,形成了一个新的基函数[公式]前面说的分窗、截断之类的都是表象,STFT的本质是基函数的改变!

那么,为什么STFT的基函数可以用于分窗,而FT的基函数不行呢?我们来看,我用正弦函数sin(5t)表示原来的基函数[公式] ,那么FT基函数和STFT基函数如下:

img

img

原因就是:FT的基函数是在时域无限延伸的,因此,无论怎么平移,都是任分布在整个时域的,起不到分窗的作用。而STFT的基函数只在时域一段不为0,在剩下的时域都是0,因此,STFT的基函数的平移,就相当于自动加了窗子啦!

紧支撑性:

这种只在时域一段不为0,在剩下的时域都是0的性质被称为“紧支撑性”(compactly supported),具有这种性质的函数,平移之后与一个信号相乘,就相当于分窗操作。这一点很重要,我们之后讲小波变换的基函数的时候还会讲。

第三点,我们前面对于分窗操作使用的函数一直称为“方窗函数”,这是一种最理想的窗函数。还有一些其他的窗函数,比如,汉宁窗、海明窗、高斯窗等。窗函数本质都是一个窗子而已,原理是一模一样的,上面所有的讨论也都成立,只是这些窗子会让信号稍稍变形一丢丢而已。你就想像方窗函数就是一面平面镜,其他的窗函数就是哈哈镜就行了。


短时傅里叶变换(STFT)的缺点与连续小波变换(CWT)


连续小波变换(CWT)的缺点与离散小波变换(DWT)

posted @ 2019-07-29 16:37  南邮果粒橙  阅读(1689)  评论(0编辑  收藏  举报