一种视频运动放大方法与流程

本发明属于图像处理方法，具体涉及一种视频运动放大方法，用于将微小的时空运动放大。

背景技术：

人眼具有有限的时空敏感性，很多真实的运动或者颜色变化不能被人眼所捕获，比如血液流经皮肤，皮肤颜色会发生轻微的变化，又比如由于心脏的周期性跳动引起的脉搏，然而人眼睛很难观察到皮肤颜色的变化和脉搏的跳动。采用视频运动放大方法，可以将这些不能引起视觉反应的运动放大，以便我们更好的观察各种现象。

Mallat和Meyer于二十世纪八十年代末提出了小波多分辨率分析(Multiresolution Analysis，MRA)，小波多分辨率分析是基于小波函数和尺度函数的缩放和平移，在Hilbert空间L²(R)的某个子空间建立基底，然后经过缩放和平移变换，将子空间的基底扩充到L²(R)中去。在小波多分辨率信号分解中，尺度函数φ_j,k和小波函数ψ_j,k共同构成了小波分解，其中j为缩放因子，k为平移因子，由尺度函数表征信号的低频部分，小波函数表征信号的细节部分。对比文件1和2，1989年7月发表，美国，S.G.Mallat,A theory for multiresolution signal decomposition:the wavelet representation,IEEE T PATTERN ANAL,11(1989)674-693；1989年9月发表，美国，S.G.Mallat,Multiresolution Approximations and Wavelet Orthonormal Bases of L²(R),T AM MATH SOC,315(1989)69-87。通过小波缩放能进行不同尺度的小波分解和重构，不同尺度的小波分解和重构构成了小波金字塔分解和重构。

现有的视频运动放大方法主要分为两种，第一种是线性欧拉视频运动放大方法，如对比文件3，2012年8月5日发表，美国，H.Wu,M.Rubinstein, E.Shih,J.Guttag,Fr,D.Durand,W.Freeman,Eulerian video magnification for revealing subtle changes in the world,ACM Trans.Graph.,31(2012)1-8。其采用拉普拉斯金字塔算法进行时空分解和重构，该方法存在两个缺陷，第一个缺陷是在视频的空间尺寸较小的时候，其只支持较小的放大因子，也就是说放大倍数有限，第二个缺陷是增加放大因子会显著增大噪声。第二种视频运动放大方法是基于复数可操作金字塔的视频运动放大方法，如对比文件4，2013年发表，美国，N.Wadhwa,M.Rubinstein,Fr,D.Durand,W.T.Freeman,Phase-based video motion processing,ACM Trans.Graph.,32(2013)1-10；该方法虽然解决了线性欧拉视频运动放大方法的两个缺陷，但是其引入了新的缺陷，第一个缺陷就是放大的运动发生失真，出现紊乱；第二个缺陷是对视频采集系统的振动较为敏感，也就是说相机的微小振动会在结果文件中体现出来。

针对欧拉视频运动放大存在对噪声比较敏感和基于相位的可操作金字塔算法运动放大在处理波传播过程中的紊乱，申请人致力于寻求一种能够解决上述问题的方法。由于小波在噪声处理方面存在优势，申请人将小波分解运用到视频运动放大中。Coiflet小波滤波器具有线性相位这一特性，在信号重构的时候可以极大的减少由于相位引起的信号失真，因而本方法采用的小波为Coiflet小波，由此申请人在前人工作的基础之上，提出了一种新的运动放大方法，称之为基于Coiflet小波的视频运动放大方法。该方法首先利用小波的多尺度分析的优点，采用小波金字塔算法对视频帧进行空间分解，将视频帧分解为不同尺度的空间频率带，然后为了提取特定的运动，采用一个与运动频率相适应的带通滤波器对空间频率带进行时域滤波，接着将滤波后的空间频带，由一个给定因子α放大，再加回到原始信号上，最后采用小波金字塔算法对视频重构，即可将微小不可见的波动放大为明显的波动。

为方便理解本发明，现对有关概念加以解释：

RGB色彩空间：一般视频帧的色彩空间都属于RGB色彩空间，RGB色彩空 间是通过对红(Red),绿(Green)和蓝(Blue)三个颜色通道的变化以及它们相互之间的叠加来得到各式各样的颜色，通常情况下，RGB色彩空间采用三个无符号整数(常用uint8，8位无符号整数)表示一种颜色，RGB各有256级亮度，用数字表示为从0、1、2...直到255，其亮度依次增加，通过各种不同的组合可以构成高达256³种颜色，因而RGB色彩空间几乎包括了人类视力所能感知的所有颜色，是目前运用最广的颜色系统之一。

YIQ色彩空间：YIQ色彩空间通常被北美的电视系统所采用，属于NTSC(National Television Standards Committee)系统。在YIQ色彩空间中，Y分量代表图像的亮度信息，I、Q两个分量则携带颜色信息，I分量代表从橙色到青色的颜色变化，而Q分量则代表从紫色到黄绿色的颜色变化。较其他颜色空间，YIQ颜色空间具有能将图像中的亮度分量分离提取出来的优点，而亮度分量可用于在自然条件下采集到的复杂背景下的运动目标的识别，因此有必要将视频由RGB色彩空间转换到YIQ色彩空间；

常见的小波函数有Haar，Daubechies，Coiflets，Symlets等；以Coiflets小波函数为例，该小波函数的缩写名为coif，根据支撑的长度不同，又分为coif1，coif2，coif3，coif4和coif5；所谓支撑长度，是指小波高通滤波器或低通滤波器是一个时间有限信号序列，在支撑长度区间内信号值非零，支撑区间外信号值为零，非零区间的长度为其支撑长度；

小波高通滤波器系数和小波低通滤波器系数：小波高通滤波器系数是一个向量，小波低通滤波器系数是另外一个向量，通过查阅相关文献可以得到对应于不同小波函数的小波高通滤波器系数和小波低通滤波器系数；如范延斌等，《小波理论算法与滤波器组》，160-192页，科学出版社，北京，2011年6月第一版。

当采用小波分解算法或小波重构算法时，对行、列进行低通滤波或高通滤 波，实质是将该行或列的元素与小波低通滤波器系数或小波高通滤波器系数求卷积，当小波低通滤波器系数或小波高通滤波器系数的长度小于该行或列的元素个数时，则需对短缺的位数补0。

计算机文件：任何一个计算机文件，都包含两个部分，第一部分为头文件部分，第二部分为数据部分。视频图像的头文件部分一般包括视频帧的宽度、高度、存储空间、格式类型、帧数率、数据压缩方式和说明部分等关键内容，视频图像的数据部分保存的是与头文件相关联的数据。计算机读取视频图像时，首先读取头文件以获取各种关键信息，尤其是数据文件的组成形式，比如视频帧的宽度、高度、帧数率等信息。在进行视频帧处理的时候，就是对数据的处理，如色彩空间变换、小波金字塔分解、小波金字塔重构和时域滤波。

技术实现要素：

本发明提供一种视频运动放大方法，属于线性欧拉视频运动放大方法，解决现有基于线性欧拉视频运动放大方法存在的增加放大因子会显著增大噪声的问题和基于复数可操作金字塔的视频运动放大方法运动失真等问题。

本发明所提供的一种视频运动放大方法，包括视频图像分解步骤、视频帧色彩空间转换步骤、N层金字塔分解步骤、时域带通滤波步骤、第二空间频率带组放大步骤、得到第四空间频率带组步骤、重构亮度矩阵步骤、视频帧色彩空间还原步骤、输出视频图像步骤，其特征在于：

(1)视频图像分解步骤：

录制运动主体在空间微小运动的RGB色彩空间视频图像，再根据时间先后将RGB色彩空间视频图像分解为一帧一帧的RGB色彩空间视频帧，读取所有的RGB色彩空间视频帧；

(2)视频帧色彩空间转换步骤：

将每一RGB色彩空间视频帧用R矩阵、G矩阵、B矩阵三个二维矩阵表示， R矩阵表示像素点红颜色强度，表示像素点绿色强度，B矩阵表示像素点蓝色强度；

通过下式对R矩阵、G矩阵和B矩阵进行矩阵运算，得到YIQ色彩空间视频帧的Y矩阵、I矩阵和Q矩阵三个二维矩阵：

Y＝0.299R+0.587G+0.114B，

I＝0.596R-0.275G-0.321B，

Q＝0.212R-0.523G+0.311B；

其中，Y矩阵表示YIQ色彩空间视频帧各像素的亮度值，I矩阵表示YIQ色彩空间视频帧各像素从橙色到青色的颜色强度，Q矩阵表示表示YIQ色彩空间视频帧各像素从紫色到黄绿色的颜色强度；

(3)N层金字塔分解步骤：

对各YIQ色彩空间视频帧中的Y矩阵，采用小波分解算法将其分解为四个矩阵：第1近似系数矩阵cA₁、第1水平细节矩阵cH₁、第1垂直细节矩阵cV₁和第1对角细节矩阵cD₁；此过程为第一层金字塔分解，cH₁、cV₁和cD₁构成金字塔第一层；

再采用小波分解算法对第1近似系数矩阵cA₁进行分解，将其分解为第2近似系数矩阵cA₂、第2水平细节矩阵cH₂、第2垂直细节矩阵cV₂和第2对角细节矩阵cD₂；此过程为第二层金字塔分解，cH₂、cV₂和cD₂构成金字塔第二层；

如此依次进行N次分解，最后得到第N近似系数矩阵cA_N、第N水平细节矩阵cH_N、第N垂直细节矩阵cV_N和第N对角细节矩阵cD_N；cA_N、cH_N、cV_N和cD_N构成金字塔第N层；N≥3，图片尺寸越大，N值越大，反之则N值越小；

对所有YIQ色彩空间视频帧进行N次金字塔分解后，由各帧YIQ色彩空间视频帧的第n水平细节矩阵cH_n组成尺度为n的水平细节空间频率带，由各帧YIQ色彩空间视频帧的第n垂直细节矩阵cV_n组成尺度为n的垂直细节空间频率带，由各帧YIQ色彩空间视频帧的第n对角细节矩阵cD_n组成尺度为n的对角细节空间频率带；其中n称为尺度，表示不同金字塔层，n＝1、2、…、 N；由各帧YIQ色彩空间视频帧的第N近似系数矩阵cA_N组成近似系数空间频率带；上述各空间频率带构成第一空间频率带组；

(4)时域带通滤波步骤：

采用第一低通IIR(Infinite Impulse Response，无限脉冲响应)数字滤波器对构成第一空间频率带组的各空间频率带分别进行时域滤波，滤波后得到的各空间频率带构成第一临时空间频率带组；

采用第二低通IIR数字滤波器对构成第一空间频率带组的各空间频率带分别进行时域滤波，滤波后得到的各空间频率带构成第二临时空间频率带组；

将第一临时空间频率带组中各空间频率带与第二临时空间频率带组中对应的各空间频率带相减，所得到的各空间频率带共同构成第二空间频率带组；

第一低通IIR数字滤波器和第二低通IIR数字滤波器的截止频率分别为f₁和f₂，0<f₂<f₁<f_s/2，f_s为视频图像的录制帧速率；

(5)第二空间频率带组放大步骤：

分别计算金字塔每一层的最大放大倍数α(n)_max，分别比较是否α(n)_max≤α，是则将α(n)_max作为相应层实际放大倍数，否则将α作为相应层实际放大倍数；其中α为设定放大倍数，5<α<30；

将构成第二空间频率带组的各空间频率带按照其所在金字塔的层数分别乘以相应层实际放大倍数，所得到的各空间频率带共同构成第三空间频率带组；

(6)得到第四空间频率带组步骤：

将构成第三空间频率带组的各空间频率带分别与所述第一空间频率带组中相应的空间频率带进行矩阵加法运算，所得到的各空间频率带共同构成第四空间频率带组；所述矩阵加法运算为两组对应尺度、对应位置的矩阵运算；

(7)重构亮度矩阵Y步骤：

采用小波重构算法将第四空间频率带组重构YIQ色彩空间视频帧的Y⁴矩阵，具体过程如下：

对第四空间频率带组中不同金字塔层的每一层进行逐层重构，从金字塔第N层开始，采用小波重构算法，将构成第四空间频率带组中各帧YIQ色彩空间视频帧的第N近似系数矩阵cA_N⁴与第N水平细节矩阵cH_N⁴、第N垂直细节矩阵cV_N⁴和第N对角细节矩阵cD_N⁴重构第N-1近似系数矩阵cA_N-1⁴；

再进行金字塔第N-1层重构，对各帧YIQ色彩空间视频帧，采用小波重构算法，将上一层重构的第N-1近似系数矩阵cA_N-1⁴与第N-1水平细节矩阵cH_N-1⁴、第N-1垂直细节矩阵cV_N-1⁴、第N-1对角细节矩阵cD_N-1⁴重构第N-2近似系数矩阵cA_N-2⁴；

如此依次进行，直到得到重构的第1近似系数矩阵cA₁⁴；然后采用小波重构算法进行金字塔第1层重构，对各帧YIQ色彩空间视频帧，采用小波重构算法，将第1近似系数矩阵cA₁⁴与第1水平细节矩阵cH₁⁴、第1垂直细节矩阵cV₁⁴、第1对角细节矩阵cD₁⁴重构为Y⁴矩阵；

(8)视频帧色彩空间还原步骤：

通过下式对步骤(7)重构的Y矩阵、步骤(2)中的I矩阵和Q矩阵进行矩阵运算，得到新的R矩阵、G矩阵和B矩阵：

R＝1.000Y⁴+0.956I+0.621Q，

G＝1.000Y⁴-0.272I-0.647Q，

B＝1.000Y⁴-1.106I+1.703Q；

新的R矩阵、G矩阵和B矩阵合成新RGB色彩空间视频帧；

依次对每一新YIQ色彩空间视频帧进行上述转换，就将YIQ色彩空间视频帧还原为RGB色彩空间视频帧；

(9)输出视频图像步骤：

构造一个视频图像文件，包括头文件部分和数据部分，首先设置视频图像文件存储位置，然后按照输入视频图像的头文件格式，构造输出视频图像的头文件，并将其写入硬盘指定位置，接着将每一帧视频帧数据硬盘指定位置，由此完成RGB色彩空间视频图像的构造。

所述的视频运动放大方法，其特征在于：

所述步骤(3)中，对Y矩阵，采用小波分解算法将其分解为四个矩阵包括下述子步骤：

(3.1)对Y矩阵的每一行进行低通滤波，然后进行列下采样，接着对每一列进行低通滤波，最后进行行下采样得到第1近似系数cA₁矩阵；

(3.2)对Y矩阵的每一行进行低通滤波，然后进行列下采样，接着对每一列进行高通滤波，最后进行行下采样得到第1水平细节cH₁矩阵；

(3.3)对Y矩阵的每一行进行高通滤波，然后进行列下采样，接着对每一列进行低通滤波，最后进行行下采样得到第1垂直细节cV₁矩阵；

(3.4)对Y矩阵的每一行进行高通滤波，然后进行列下采样，接着对每一列进行高通滤波，最后进行行下采样得到第1对角细节cD₁矩阵；

采用小波分解算法对第n近似系数矩阵cA_n-1进行分解，子步骤与上述相同，n＝1、2、…、N；

所述列下采样就是保留偶数列，丢弃奇数列；所述行下采样就是保留偶数行，丢弃奇数行；

所述对行进行低通滤波就是将该行元素与小波低通滤波器系数求卷积，得到一行新元素值，代替原有行；

所述对行进行高通滤波就是将该行元素与小波高通滤波器系数求卷积，得到一行新元素值，代替原有行；

所述对列进行低通滤波就是将该列元素与小波低通滤波器系数求卷积，得到一列新元素值，代替原有列；

所述对列进行高通滤波就是将该列元素与小波高通滤波器系数求卷积，得到一列新元素值，代替原有列；

所述小波低通滤波器系数和小波高通滤波器系数通过查阅相关文献可以得到。

所述的视频运动放大方法，其特征在于：

所述步骤(4)中，采用第一低通IIR数字滤波器对构成第一空间频率带组的各空间频率带分别进行时域滤波包括下述子步骤：

(4.1)计算低通IIR滤波器的滤波器系数r₁：

r₁＝2πf₁/f_s，

式中f₁为第一低通IIR数字滤波器的截止频率，f_s为视频图像的录制帧速率；

(4.2)计算第一低通IIR数字滤波器的输出Y(m)：

Y(m)＝(1-r₁)×Y(m-1)+r₁×X(m)，

式中，X(m)为滤波器的输入，m为视频帧序号，m＝1、2、…、K；其中K为视频帧的总帧数，当m为1，X(1)是已知的，Y(0)为：

采用第二低通IIR数字滤波器对构成第一空间频率带组的各空间频率带分别进行时域滤波，其过程与上述过程完全相同，区别仅在于将f₁换成f₂。

所述的视频运动放大方法，其特征在于：

所述步骤(5)中，计算金字塔每一层的最大放大倍数α(n)_max包括下述子步骤：

(5.1)计算金字塔每一层的空间波长λ(n)：

式中，W_n和H_n分别为金字塔第n层的宽与高，单位为像素；

(5.2)计算位移函数δ(t)：

式中，λ_c为空间临界波长，取值为16～20像素；

(5.3)计算金字塔每一层的最大放大倍数α(n)_max：

所述的视频运动放大方法，其特征在于：

所述步骤(7)中，采用小波重构算法，将构成第四空间频率带组中各帧YIQ色彩空间视频帧的第N近似系数矩阵cA_N⁴与第N水平细节矩阵cH_N⁴、第N垂直细节矩阵cV_N⁴和第N对角细节矩阵cD_N⁴重构第N-1近似系数矩阵cA_N-1⁴，包括下述子步骤：

(7.1)对矩阵cA_N⁴进行行上采样，然后对行上采样后矩阵的每一列进行低通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行低通滤波得到A_N1矩阵；

(7.2)对矩阵cH_N⁴进行行上采样，然后对行上采样后矩阵的每一列进行高通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行低通滤波得到A_N2矩阵；

(7.3)对cV_N⁴矩阵进行行上采样，然后对行上采样后矩阵的每一列进行低通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行高通滤波得到A_N3矩阵；

(7.4)对cD_N⁴矩阵进行行上采样，然后对行上采样后矩阵的每一列进行高通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行高通滤波得到A_N4矩阵；

(7.5)将上述A_N1矩阵、A_N2矩阵、A_N3矩阵和A_N4矩阵进行矩阵加法运算得到第N-1近似系数矩阵cA_N-1⁴；

所谓矩阵行上采样，就是将矩阵的行数增加一倍，在原始矩阵任意两行之间插入一行零向量，最后对新构矩阵的最后一行补零；所谓矩阵列上采样，就是将矩阵的列数增加一倍，在原始矩阵任意两列之间插入一列零向量，最 后对新构矩阵的最后一列补零；

所述对行进行低通滤波就是将该行元素与小波低通滤波器系数求卷积，得到一行新元素值，代替原有行；

所述对行进行高通滤波就是将该行元素与小波高通滤波器系数求卷积，得到一行新元素值，代替原有行；

所述对列进行低通滤波就是将该列元素与小波低通滤波器系数求卷积，得到一列新元素值，代替原有列；

所述对列进行高通滤波就是将该列元素与小波高通滤波器系数求卷积，得到一列新元素值，代替原有列；

所述小波低通滤波器系数和小波高通滤波器系数通过查阅相关文献可以得到。

视频图像的相位梯度与视频图像中的运动有直接的关系，现有基于相位的复数可操作金字塔运动放大方法的基本思想是将图像的相位梯度进行放大，达到将图像中的运动放大的目的，由于图像的相位梯度发生畸变，导致放大的运动也会发生畸变。本发明并没有涉及图像的相位梯度，因而不存在运动发生畸变。

现有的线性欧拉视频运动放大方法存在增加放大因子会显著增大噪声的问题，本发明在步骤(4)中，第一低通IIR数字滤波器的差分方程为：Y(m)＝(1-r₁)×Y(m-1)+r₁×X(m)，实际操作中，现有线性欧拉视频运动放大方法将滤波器的初始输入X(1)和Y(0)进行了舍弃，直接从X(2)和Y(1)进行时域滤波，其中令Y(1)＝X(1)，具体计算过程为：

Y(1)＝X(1)

Y(2)＝(1-r₁)*Y(1)+r₁*X(2)

Y(3)＝(1-r₁)*Y(2)+r₁*X(3)

Y(m)＝(1-r₁)*Y(m-1)+r₁*X(m)

Y(K)＝(1-r₁)*Y(K-1)+r₁*X(K)

本发明步骤(4)的具体计算过程为：

Y(1)＝(1-r₁)*Y(0)+r₁*X(1)

Y(2)＝(1-r₁)*Y(1)+r₁*X(2)

Y(m)＝(1-r₁)*Y(m-1)+r₁*X(m)

Y(K)＝(1-r₁)*Y(K-1)+r₁*X(K)

通过对现有线性欧拉视频运动放大方法的输出视频进行噪声分析，发现其噪声呈现初始小，然后快速增加，接着缓慢衰减的趋势。采用本发明的试验结果显示，输出视频中的噪声很平稳，基本保持不变，另外从噪声的大小来看，本发明极大的减少了输出视频中的噪声。

综上所述，本发明解决了现有基于线性欧拉视频运动放大方法存在的增加放大因子会显著增大噪声的问题和现有基于复数可操作金字塔的视频运动放大方法运动失真等问题。

附图说明

图1为本发明的流程示意图；

图2为原始视频帧序列；

图3为小波金字塔分解示意图；

图4(A)为原始图片；

图4(B)为三层小波金字塔结构示意图；

图5为第一空间频率带组示意图；

图6为小波分解示意图，

图7(A)为第一低通IIR数字滤波器的频响特性示意图；

图7(B)为第二低通IIR数字滤波器的频响特性示意图；

图7(C)为使用第一、第二低通IIR数字滤波器两个不同截止频率的低通滤波器所构造的带通滤波器的频响特性示意图；

图8为第四空间频率带组；

图9为小波重构示意图；

图10为重构的视频序列。

具体实施方式

以下结合附图和实施例，对本发明进一步说明。

如图1所示，本发明包括视频图像分解步骤、视频帧色彩空间转换步骤、N层金字塔分解步骤、时域带通滤波步骤、第二空间频率带组放大步骤、得到第四空间频率带组步骤、重构亮度矩阵步骤、视频帧色彩空间还原步骤、输出视频图像步骤。

本发明的一个实施例，包括下述步骤：

(1)视频图像分解步骤：

录制运动主体在空间微小运动的RGB色彩空间视频图像，再根据时间先后将RGB色彩空间视频图像分解为一帧一帧的RGB色彩空间视频帧，读取所有的RGB色彩空间视频帧；如图2所示，申请人采用单反相机录制了两颗银杏树在空中微小运动的视频，视频的帧速率f_s为25帧每秒，视频大小为512×512像素。

(2)视频帧色彩空间转换步骤：

将每一RGB色彩空间视频帧用R矩阵、G矩阵、B矩阵三个二维矩阵表示， R矩阵表示像素点红颜色强度，表示像素点绿色强度，B矩阵表示像素点蓝色强度；

通过下式对R矩阵、G矩阵和B矩阵进行矩阵运算，得到YIQ色彩空间视频帧的Y矩阵、I矩阵和Q矩阵三个二维矩阵：

Y＝0.299R+0.587G+0.114B，

I＝0.596R-0.275G-0.321B，

Q＝0.212R-0.523G+0.311B；

其中，Y矩阵表示YIQ色彩空间视频帧各像素的亮度值，I矩阵表示YIQ色彩空间视频帧各像素从橙色到青色的颜色强度，Q矩阵表示表示YIQ色彩空间视频帧各像素从紫色到黄绿色的颜色强度；

(3)N层金字塔分解步骤：

对各YIQ色彩空间视频帧中的Y矩阵，采用小波分解算法将其分解为四个矩阵：第1近似系数矩阵cA₁、第1水平细节矩阵cH₁、第1垂直细节矩阵cV₁和第1对角细节矩阵cD₁；此过程为第一层金字塔分解，cH₁、cV₁和cD₁构成金字塔第一层；

再采用小波分解算法对第1近似系数矩阵cA₁进行分解，将其分解为第2近似系数矩阵cA₂、第2水平细节矩阵cH₂、第2垂直细节矩阵cV₂和第2对角细节矩阵cD₂；此过程为第二层金字塔分解，cH₂、cV₂和cD₂构成金字塔第二层；

如此依次进行N次分解，最后得到第N近似系数矩阵cA_N、第N水平细节矩阵cH_N、第N垂直细节矩阵cV_N和第N对角细节矩阵cD_N；cA_N、cH_N、cV_N和cD_N构成金字塔第N层；N≥3，图片尺寸越大，N值越大，反之则N值越小；

对所有YIQ色彩空间视频帧进行N次金字塔分解后，由各帧YIQ色彩空间视频帧的第n水平细节矩阵cH_n组成尺度为n的水平细节空间频率带，由各帧YIQ色彩空间视频帧的第n垂直细节矩阵cV_n组成尺度为n的垂直细节空间频率带，由各帧YIQ色彩空间视频帧的第n对角细节矩阵cD_n组成尺度为n的对角细节空间频率带；其中n称为尺度，表示不同金字塔层，n＝1、2、…、 N；由各帧YIQ色彩空间视频帧的第N近似系数矩阵cA_N组成近似系数空间频率带；上述各空间频率带构成第一空间频率带组；

图3表示采用小波金字塔算法进行小波金字塔分解图像数据的存储形式，该图是对一幅尺寸为512×512像素的图片进行小波金字塔分解，C表示分解后矩阵的个数，S表示各矩阵的尺寸信息。图4(B)表示小波分解后的三层小波金字塔结构示意图，首先采用小波分解算法将图片分解为cA₁，cH₁，cV₁和cD₁，图中左上角为cA₁，将cA₁再分解为第二层金字塔cA₂，cH₂，cV₂，cD₂,cA₁被cA₂，cH₂，cV₂，cD₂取代，将cA₂再分解为第三层金字塔cA₃，cH₃，cV₃，cD₃，cA₂被cA₃，cH₃，cV₃，cD₃取代，于是就构成三层金字塔的结构；图5表示了第一空间频率带组，此图结构布局类似于图4(B)，为了简化，仅仅进行了1层金字塔分解，金字塔的总层数为1；

如图6所示，采用小波分解算法将Y矩阵其分解为四个矩阵包括下述子步骤：

(3.1)对Y矩阵的每一行进行低通滤波，然后进行列下采样，接着对每一列进行低通滤波，最后进行行下采样得到第1近似系数cA₁矩阵；

(3.2)对Y矩阵的每一行进行低通滤波，然后进行列下采样，接着对每一列进行高通滤波，最后进行行下采样得到第1水平细节cH₁矩阵；

(3.3)对Y矩阵的每一行进行高通滤波，然后进行列下采样，接着对每一列进行低通滤波，最后进行行下采样得到第1垂直细节cV₁矩阵；

(3.4)对Y矩阵的每一行进行高通滤波，然后进行列下采样，接着对每一列进行高通滤波，最后进行行下采样得到第1对角细节cD₁矩阵；

采用小波分解算法对第n近似系数矩阵cA_n-1进行分解，子步骤与上述相同，n＝1、2、…、N；

所述列下采样就是保留偶数列，丢弃奇数列；所述行下采样就是保留偶数行，丢弃奇数行；

所述对行进行低通滤波就是将该行元素与小波低通滤波器系数求卷积， 得到一行新元素值，代替原有行；

所述对行进行高通滤波就是将该行元素与小波高通滤波器系数求卷积，得到一行新元素值，代替原有行；

所述对列进行低通滤波就是将该列元素与小波低通滤波器系数求卷积，得到一列新元素值，代替原有列；

所述对列进行高通滤波就是将该列元素与小波高通滤波器系数求卷积，得到一列新元素值，代替原有列；

本实施例中，所采用的小波函数为coif5，通过查阅参考文件1和参考文件2，得到低通滤波器和高通滤波器，如表1所列；

参考文件1：中国，轩建平等，“高阶Coiflet小波系构造与应用”，振动工程学报，(2008)521-529；参考文件2：禹如杰，“基于零矩尺度函数的波传播问题求解方法研究”，华中科技大学硕士论文，1-52，2013年发表；

(4)时域带通滤波步骤：

采用第一低通IIR数字滤波器(Infinite Impulse Response，无限脉冲响应数字滤波器)对构成第一空间频率带组的各空间频率带分别进行时域滤波，滤波后得到的各空间频率带构成第一临时空间频率带组；

采用第二低通IIR数字滤波器对构成第一空间频率带组的各空间频率带分别进行时域滤波，滤波后得到的各空间频率带构成第二临时空间频率带组；

将第一临时空间频率带组中各空间频率带与第二临时空间频率带组中对应的各空间频率带相减，所得到的各空间频率带共同构成第二空间频率带组；

本实施例中，第一低通IIR数字滤波器和第二低通IIR数字滤波器的截止频率f₁和f₂分别为0.4Hz和0.2Hz，图7(A)表示第一低通IIR数字滤波器的频响特性、图7(B)表示第二低通IIR数字滤波器的频响特性、图7(C)表示由第一和第二低通IIR数字滤波器构成的带通滤波器的频响特性。

(5)第二空间频率带组放大步骤：

分别计算金字塔每一层的最大放大倍数α(n)_max，分别比较是否α(n)_max ≤α，是则将α(n)_max作为相应层实际放大倍数，否则将α作为相应层实际放大倍数；其中α为设定放大倍数，5<α<30；

将构成第二空间频率带组的各空间频率带按照其所在金字塔的层数分别乘以相应层实际放大倍数，所得到的各空间频率带共同构成第三空间频率带组；

本实施例中，计算金字塔每一层的最大放大倍数α(n)_max包括下述子步骤：

(5.1)计算金字塔每一层的空间波长λ(n)：

式中，W_n和H_n分别为金字塔第n层的宽与高，单位为像素；以第1层和第2层为例，第1层宽W₁与高H₁分别为256和256，第二层的宽W₂与高H₂分别为128和128，参考图3。

(5.2)计算位移函数δ(t)：

式中，λ_c为空间临界波长，取值为16像素，设定放大倍数α取为15；

(5.3)计算金字塔每一层的最大放大倍数α(n)_max：

以第1层和第2层为例，计算最大放大倍数的具体过程为：首先利用已知第一层的宽W₁与高H₁分别为256和256，第二层的宽W₂与高H₂分别为128和128，将其代入计算公式，计算每一层的空间波长，得到λ(1)＝361.984和λ(2)＝180.992，然后根据设定的空间临界波长λ_c＝16和设定放大倍数α＝15计算得到δ(t)＝0.125，最后将这三个计算值代入，求得第一层和第二层的最大放大倍数分别为：α(1)_max＝119.661，α(2)_max＝59.331。

(6)得到第四空间频率带组步骤：

将构成第三空间频率带组的各空间频率带分别与所述第一空间频率带组中相应的空间频率带进行矩阵加法运算，所得到的各空间频率带共同构成第 四空间频率带组；所述矩阵加法运算为两组对应尺度、对应位置的矩阵运算；

如图8所示，该图就是叠加后的空间频率带。此叠加过程是线性欧拉的核心近似部分，在线性欧拉视频运动放大中，其核心思想可以表述为：

假设一个一维图像中任意一个像素点的强度I(x)随着时间t的变化为I(x,t)，将初始值写为I(x,0)＝h(x)，那么任意时刻该像素点的强度可以表示为I(x,t)＝h(x+δ(t))，当物体的运动为小运动的时候，可以通过泰勒近似展开如式1所示：

在公式1中，h(x+δ(t))为图像原始像素点强度，表征着第一空间频率带，为表征物体的运动像素点强度变化，为了获取感兴趣的运动，可以用一个带通滤波器对第一空间频率带进行时域带通滤波。通过这种时域处理得到第二空间频率带，第二空间频率带可以用表示，将第二空间频率带的强度值乘以一个放大因子α，于是就得到了第三空间频率带，用表示第三空间频率带，再将其加到原始第一空间频率带中，用下式表示该相加过程：

公式2右边还可以近似表示为：

通过回加，就得到第四空间频率带，通过对比公式1的左边和公式3的右边，可以发现表征随着时间变化的运动被放大了(1+α)倍。

(7)重构亮度矩阵步骤：

采用小波重构算法将第四空间频率带组重构YIQ色彩空间视频帧的Y⁴矩阵，具体过程如下：

对第四空间频率带组中不同金字塔层的每一层进行逐层重构，从金字塔第N层开始，采用小波重构算法，将构成第四空间频率带组中各帧YIQ色彩空间视频帧的第N近似系数矩阵cA_N⁴与第N水平细节矩阵cH_N⁴、第N垂直细节 矩阵cV_N⁴和第N对角细节矩阵cD_N⁴重构第N-1近似系数矩阵cA_N-1⁴；

再进行金字塔第N-1层重构，对各帧YIQ色彩空间视频帧，采用小波重构算法，将上一层重构的第N-1近似系数矩阵cA_N-1⁴与第N-1水平细节矩阵cH_N-1⁴、第N-1垂直细节矩阵cV_N-1⁴、第N-1对角细节矩阵cD_N-1⁴重构第N-2近似系数矩阵cA_N-2⁴；

如此依次进行，直到得到重构的第1近似系数矩阵cA₁⁴；然后采用小波重构算法进行金字塔第1层重构，对各帧YIQ色彩空间视频帧，采用小波重构算法，将第1近似系数矩阵cA₁⁴与第1水平细节矩阵cH₁⁴、第1垂直细节矩阵cV₁⁴、第1对角细节矩阵cD₁⁴重构为Y⁴矩阵；

如图9所示，本实施例中，采用小波重构算法，将构成第四空间频率带组中各帧YIQ色彩空间视频帧的第N近似系数矩阵cA_N⁴与第N水平细节矩阵cH_N⁴、第N垂直细节矩阵cV_N⁴和第N对角细节矩阵cD_N⁴重构第N-1近似系数矩阵cA_N-1⁴，包括下述子步骤：

(7.1)对矩阵cA_N⁴进行行上采样，然后对行上采样后矩阵的每一列进行低通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行低通滤波得到A_N1矩阵；

(7.2)对矩阵cH_N⁴进行行上采样，然后对行上采样后矩阵的每一列进行高通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行低通滤波得到A_N2矩阵；

(7.3)对cV_N⁴矩阵进行行上采样，然后对行上采样后矩阵的每一列进行低通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行高通滤波得到A_N3矩阵；

(7.4)对cD_N⁴矩阵进行行上采样，然后对行上采样后矩阵的每一列进行高通滤波，接着对滤波后矩阵进行列上采样，再接着对列上采样得到矩阵的每一行进行高通滤波得到A_N4矩阵；

(7.5)将上述A_N1矩阵、A_N2矩阵、A_N3矩阵和A_N4矩阵进行矩阵加法运算得 到第N-1近似系数矩阵cA_N-1⁴；

所谓矩阵行上采样，就是将矩阵的行数增加一倍，在原始矩阵任意两行之间插入一行零向量，最后对新构矩阵的最后一行补零；所谓矩阵列上采样，就是将矩阵的列数增加一倍，在原始矩阵任意两列之间插入一列零向量，最后对新构矩阵的最后一列补零；

所述对行进行低通滤波就是将该行元素与小波低通滤波器系数求卷积，得到一行新元素值，代替原有行；

所述对行进行高通滤波就是将该行元素与小波高通滤波器系数求卷积，得到一行新元素值，代替原有行；

所述对列进行低通滤波就是将该列元素与小波低通滤波器系数求卷积，得到一列新元素值，代替原有列；

所述对列进行高通滤波就是将该列元素与小波高通滤波器系数求卷积，得到一列新元素值，代替原有列；

本实施例中，所采用的小波函数为coif5，通过查阅前述参考文件1和参考文件2，得到小波低通滤波器系数和小波高通滤波器系数如表1所列；

(8)视频帧色彩空间还原步骤：

通过下式对步骤(7)重构的Y⁴矩阵、步骤(2)中的I矩阵和Q矩阵进行矩阵运算，得到新的R矩阵、G矩阵和B矩阵：

R＝1.000Y⁴+0.956I+0.621Q，

G＝1.000Y⁴-0.272I-0.647Q，

B＝1.000Y⁴-1.106I+1.703Q；

新的R矩阵、G矩阵和B矩阵合成新RGB色彩空间视频帧；

依次对每一新YIQ色彩空间视频帧进行上述转换，就将YIQ色彩空间视频帧还原为RGB色彩空间视频帧；

(9)输出视频图像步骤：

构造一个视频图像文件，包括头文件部分和数据部分，首先设置视频图 像文件存储位置，然后按照输入视频图像的头文件格式，构造输出视频图像的头文件，并将其写入硬盘指定位置，接着将每一帧视频帧数据硬盘指定位置，由此完成RGB色彩空间视频图像的构造。所构造的RGB色彩空间视频图像如图10所示，将其与图2对比，可以很明显的看到视频中的运动被放大了。

表1 coif5低通滤波器系数和高通滤波器系数