一种基于离散分数阶傅里叶变换相位信息的信号重建方法

摘要

本发明公开了一种仅仅通过离散分数阶傅里叶变换的相位信息来重建原始信号的方法，属于信号处理技术领域。本发明首先将信号重建问题转化为凸优化问题；然后，对原始信号进行离散分数阶傅里叶变换，并通过改变离散分数阶傅里叶变换的变换矩阵获得不同数目的相位信息；接着，将得到的相位信息进行存储或者传输；最后，利用块坐标下降法和内点法结合的幅度恢复算法，通过合适数目的相位信息将原始信号恢复出来，即重建原始信号。本发明方法利用相同数目下的相位信息包含的信息量大于幅度信息包含的信息量这一理论依据，实现了以较少数目的相位信息重建原始信号的目的。

说明书

技术领域

本发明涉及一种仅仅通过离散分数阶傅里叶变换域的相位信息来重建原始 信号的方法，属于信号处理技术领域。

背景技术

自然界中的场景都是三维的，然而照片和电脑屏幕显示的是伪三维的二维平 面场景。天文图像检测和地震信号检测，获取的只是这些图像的傅里叶振幅信息， 而观察不到希望得到的星云图像或地震反射序列信号。X射线晶体结构分析时， 分子结构图中只有结构因子的绝对值大小。电子显微镜或光学显微镜确定物体结 构时，相位信息也会明显丢失。然而相位信息中包含图像大量的信息，因此如何 利用已知的幅度信息来恢复丢失的相位信息，进而重建图像，成为需要解决的问 题。这个问题也被称为相位恢复问题，相应的解决方法我们称为相位恢复算法。 一些情况下，图像长期暴露在大气环境或被带有圆形孔径的散焦镜头照射等情 况，会使得期望获得的图像的傅里叶幅度被破坏，仅知道变换域的相位信息。因 而又引出了利用变换域相位信息恢复幅度信息，最终重建图像的问题。这个问题 我们称为幅度恢复问题，相应的解决方法我们称为幅度恢复算法。无论是相位恢 复问题，还是幅度恢复问题，都属于信号重建。通常变换域的相位和幅度相互独 立，仅从任何一项信息(相位或幅度)重建原始信号似乎是不可能的。然而，Hayes 和Oppenheim证明了在某些条件下，上述信号重建问题是可以解决的。

至今，国内外许多专家和学者已经提出了一些解决方法。相位恢复算法如 1968年Schiske提出的系列像重构算法、1972年Gerchberg和Saxton提出的GS 算法、1982年Fienup提出的Error Reduction算法、1994年Schewchuk提出的 Conjugate Gradient算法、1994年Yang和Gu提出的Y-G算法，然而这些算法都 是光学系统中由光波强度通过迭代得到相位信息，并重建图像的算法。尽管上述 算法在某些特定的情况下可以较好的重构图像，但并不能总是得到一致的鲁棒的 结果。幅度恢复算法包括1961年Srinivasan提出的直接法、1968年Goodman和 Knight提出的统计算法、1980年Hayes和Oppenheim提出的迭代算法和闭形算 法、1983年Levi和Stark提出的凸集投影算法和1992年Behar提出的已知部分 相位的信号重建算法。2008年Gang和Orchard提出的基于几何模型估计的信号 重建算法。2010年Loveimi等人通过最小二乘误差估计和重叠添加算法实现了 语音信号的重建。2013年，Boufounos用标准凸优化和贪婪算法实现了信号的重 建，理论和实验的结果都表明从相位信息精确重建原始信息是可能的。上述的这 些算法已经被用于声学光学全息技术领域、微电子技术领域、语音信号处理领域 和X射线晶体结构测量领域。然而，上述算法的先验条件都是针对原始信号， 如信号具有稀疏特性，信号Z变换的零极点分布特性等。

近年来，随着压缩感知、矩阵填充理论和凸优化技术的深入研究和发展，相 位恢复问题再次被提出。2008年Candes等人仅通过变换域的幅度信息提出了基 于矩阵填充理论的PhaseLift信号重建算法，该算法不仅将幅度的采集个数减少 为nlogn，而且在噪声存在的条件下具有鲁棒性。随后，2012年Waldspurger提 出了基于凸优化技术的PhaseCut算法，对音频信号处理的仿真结果表明该算法 的特性优于PhaseLift算法。PhaseLift算法和PhaseCut算法并没有添加原始信号 的先验信息，而是增加了变换矩阵A的维数。当变换矩阵的维数足够时，就可以 解决相位恢复问题。

对于先验信息添加在变换矩阵A的幅度恢复算法目前包括经典的 Gerchberg-Saxton算法和贪婪算法。

Gerchberg-Saxton算法：首先，找到一个合适的初始值，如X⁰＝diag(u)b。 然后，逐个去优化变换域的每个元素，即从而得 到较优的解X^N。最后，重建原始信号

贪婪算法：如权利要求书中所述将原始问题转化为其 中b是待优化的变换域幅度，分项写下问题：

$e_{2k}\left[n\right]=P{\left[V_{e_k}^T\right|0...0]}^T,e_{2k+1}\left[n\right]=P{[0...0|V_{o_k}^T]}^T,$

这个问题的解为：

$F^a[m,n]={\Sigma }_{k\in M}{e}_k\left[m\right]e^{-j\frac{\pi }{2}\mathrm{ka}}{e}_k\left[n\right],M=0\{0,...,p-2,(p-{\left(p\right)}_2)\}.$

优化b，首先随机一个初始值再利用式$V_e=\left(\begin{array}{ccc}{V}_{e_1}& ···& V_{e_k}\end{array}\right)$ k＝1,...,N逐个优化幅度值。最后得到较优的解

然而，上述的Gerchberg-Saxton算法和贪婪算法需要较多数目的相位信息才 可以较好的重建原始信号，而且这两种算法求解的稳定性较差。为了解决这些问 题，我们通过压缩感知和矩阵填充理论将幅度恢复问题重述为凸优化问题，利用 凸优化技术推导出幅度恢复算法，并最终实现较少数目相位信息下的稳定的信号 重建。

发明内容

发明目的：针对现有技术中存在的问题与不足，本发明提供一种基于离散分 数阶傅里叶变换相位信息的信号重建方法，能够实现从相位信息恢复幅度信息并 最终重建原始信号，解决已知任意信号离散分数阶变换域的相位信息来重建原始 信号的问题。

技术方案：一种基于离散分数阶傅里叶变换相位信息的信号重建方法，包括 将信号重建问题转化为一个凸优化问题，并利用“幅度恢复算法”解决该问题， 最终实现从离散分数阶傅里叶变换合适数目的相位信息到原始信号的重建，具体 包括以下步骤：

步骤A、大小为p的原始信号x，求解它的离散分数阶傅里叶变换；对于离 散分数阶傅里叶变换，采用特征值和特征向量分解方法，可知Hamiltonian矩阵 S为：

当p为偶数时，可以将任意向量分解为奇偶两部分的矩阵P为：

$P=\frac{1}{\sqrt{2}}\left(\begin{array}{cccc}\sqrt{2}& 0& 0& 0\\ 0& I_r& 0& J_r\\ 0& 0& 1& 0\\ 0& J_r& 0& -I_r\end{array}\right),$

当p为奇数时，矩阵P为

$P=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}\sqrt{2}& 0& 0\\ 0& I_r& J_r\\ 0& J_r& -I_r\end{array}\right).$

这里的I_r是r×r的单位矩阵(对角线上的值为1)，J_r是r×r的反单位矩阵(反 对角线上的值为1)。可以看到矩阵P满足P＝P^T＝P^-1。

步骤B、已知上述的矩阵S和P，计算PSP^T，且矩阵PSP^T是一个块对角矩 阵，于是分块得到${\mathrm{PSP}}^T=\left(\begin{array}{cc}\mathrm{Ev}& 0\\ 0& \mathrm{Od}\end{array}\right).$其中矩阵Ev和Od分别是矩阵PSP^T的偶特 征向量和奇特征向量。

步骤C、对矩阵上述两个矩阵Ev和Od做特征值分解，结果为Ev＝V_eΛ_eV_e^T和 Od＝V_oΛ_oV_o^T，

其中，(Ev的特征值e₁,...,e_k组成的矩阵)，

(Od的特征值o₁,...,o_k组成的矩阵)，

(Ev的特征向量组成的矩阵)，

(Od的特征向量组成的矩阵)。

再分别对特征值做降序排列，并对对应的特征向量排列得到： $e_{2k}\left[n\right]=P{\left[V_{e_k}^T\right|0...0]}^T,e_{2k+1}\left[n\right]=P{[0...0|V_{o_k}^T]}^T,$这里的e_2k[n]和e_2k+1[n]是离散 Hermite高斯矩阵，它是定义离散分数阶傅里叶变换的重要部分。

步骤D、计算变换矩阵F^a，

$F^a[m,n]={\Sigma }_{k\in M}{e}_k\left[m\right]e^{-j\frac{\pi }{2}\mathrm{ka}}{e}_k\left[n\right],M=\{0,...,p-2,(p-{\left(p\right)}_2)\}.$其中e_k是步骤C 中的离散Hermite高斯矩阵，是与k和a有关的指数函数，(p)₂≡pmod2， 即p除以2取余数。

最后原始信号x的离散分数阶傅里叶变换为X_a＝F^ax。令$E=P\left(\begin{array}{cc}{V}_e& 0\\ 0& V_o\end{array}\right),$则步骤D中的变换矩阵F^a可以简化为F^a＝Ediag(Λ)^aE^T，这里所需 要的计算复杂度为Ο(p³)。为了减少计算复杂度，将离散分数阶傅里叶变换写做 X_a＝F^ax＝E(Λ^a(E^Tx))，这时的计算复杂度仅为Ο(p³)。

步骤E、将原始信号重建问题转化为数学问题，如下所示：

$\left(\begin{array}{ccc}\mathrm{find}& & x\\ \mathrm{such}& \mathrm{that}& e^{j\angle \mathrm{Ax}}=u\end{array}\right)---\left(1\right)$

其中，x是大小为p的原始信号，是离散分数阶 傅里叶变换矩阵，是指维数为n×p的复数矩阵，∠Ax是Ax的相位角， 是变换域的相位信息。问题的目的就是找到一个合适的n 使原始信号得以重建。

步骤F、由于变换域的相位(维数为n的复数向量)和幅度(维数 为n的实数向量)的乘积diag(u)b是完整的变换域信号，那么如果恢复优化幅度 以至于满足Ax＝diag(u)b，那么这个优化后的幅度就是要求的最优解，从而精 确重建原始信号。于是，(1)可以转化为如下最小二乘问题；

我们知道(2)中的x可以近似为即(2)可以转化为；

进一步化简目标函数得：

其中，(·)^H是矩阵的共轭转置，是矩阵 的伪逆。

步骤G、若令(3)简化为；

再令表示正对称矩阵，于是(4)化简为；

$\left(\begin{array}{ccc}\min & \mathrm{Tr}\left(\mathrm{BM}\right)& \\ s.t.& B\ge 0,B\in S_n^{+},& \mathrm{rank}\left(B\right)=1\end{array}---\right)\left(5\right)$

(5)是一个非凸的秩约束问题，如果将秩约束放在后面考虑，那么可以近似 认为其是凸约束问题，如下所示：

minTr(BM),s.t.B≥0                   (6)

为了将约束条件添加到目标函数中，这里采用了对数障碍约束，μ是小于1 的约束系数，如下所示；

minTr(BM)-μlogdet(B),μ＞0         (7)

然后再利用幅度恢复算法解决上述凸优化问题，进而重建原始信号。

步骤H、首先，在算法的外层采用“块坐标下降算法”迭代，将矩阵B划分为 四块，$B=\left(\begin{array}{cc}P& y\\ y^T& b^2\end{array}\right),$$y=\left(\begin{array}{c}{b}_1{b}_n\\ .\\ .\\ .\\ b_{n-1}{b}_n\end{array}\right);b^2=b_n^2.$

由Powell提出的块矩阵行列式求解的结论 $\left(\begin{array}{cc}P& y\\ y^T& b^2\end{array}\right)\left(\begin{array}{cc}I& -P^{-1}y\\ 0& I\end{array}\right)=\left(\begin{array}{cc}P& 0\\ y^T& b^2-y^T{P}^{-1}y\end{array}\right),$得到矩阵B的行列式值为 det(B)＝det(P)det(b²-y^TP^-1y)。

步骤I、因为矩阵M和B都是Hermitian矩阵，所以可以把复数问题转化成 实数问题求解。对于两个Hermitian矩阵M和B来说，他们满足公式 2Tr(BM)＝Tr(Γ(B)Γ(M))，其中$\Gamma (·)=\left(\begin{array}{cc}\mathrm{Re}(·)& -\mathrm{Im}(·)\\ \mathrm{Im}(·)& \mathrm{Re}(·)\end{array}\right).$于是上述公式化简为： Tr(Γ(B)Γ(M))＝2Tr(BRe(M))，将其代入(7)得：

minTr(BRe(M))-μlogdet(B)                      (8)

步骤J、下面利用块坐标下降算法解决上述问题(8)，令Re(M)＝R，简化(8) 如下：

其中i^c∈{1,...,i-1,i+1,...,n},i＝{1,...,n}。再令问题(9)的目标函 数可以表示为：

$F({y^{\prime }}_{i^{c}i},b_i)=R_{i^{c}i}^T{y^{\prime }}_{i^{c}i}{b}_i+R_{\mathrm{ii}}{b}_i^2-\mu \log (1-{y^{\prime }}_{i^{c}i}^T{P}_{i^c{i}^c}^{-1}{y^{\prime }}_{i^{c}i})-\mu \log b_i^2.$

步骤K、为了得到(9)的最小值，在内层使用log-barrier内点法迭代解决上面 的凸优化问题。方便起见，令再分别对y'和b_i求 F(y',b_i)的偏导。

$\frac{\partial F(y^{\prime },b_i)}{{\partial y}^{\prime }}={\mathrm{Qb}}_i+\frac{{2\mathrm{\mu P}}^{-1}{y}^{\prime }}{(1-y^{\prime T}{P}^{-1}{y}^{\prime })}=I\left(y^{\prime }\right)---\left(10\right)$

$\frac{\partial F(y^{\prime },b_i)}{{\partial b}_i}{=Q}^T{y}^{\prime }+{2R}_{\mathrm{ii}}{b}_i-\frac{2\mu }{b_i}=J\left(b_i\right)---\left(11\right)$

然后，将式(10)和(11)进行二阶泰勒展开：

$I(y^{\prime }+{\mathrm{\Delta y}}^{\prime })=I\left(y^{\prime }\right)+\frac{2\mu {\left(P^{-1}\right)}^T}{(1-y^{\prime T}{P}^{-1}{y}^{\prime })}{\mathrm{\Delta y}}^{\prime }+\frac{{4\mathrm{\mu P}}^{-1}{y}^{\prime }{\left(P^{-1}{y}^{\prime }\right)}^T}{{(1-y^{\prime T}{P}^{-1}{y}^{\prime })}^2}{\mathrm{\Delta y}}^{\prime }---\left(12\right)$

$J(b_i+\Delta b_i)=J\left(b_i\right)+{2R}_{\mathrm{ii}}{\mathrm{\Delta b}}_i+\frac{2\mu }{b_i^2}{\mathrm{\Delta b}}_i---\left(13\right)$

令

G＝-μlog(1-y'^TP^-1y')

g_y'＝dGdy'＝2μP^-1y'(1-y'^TP^-1y')

H_y'＝d²Gdy'²＝2μ(P^-1)^T(1-y'^TP^-1y')+4μP^-1y'(P^-1y')^T(1-y'^TP^-1y')²

和$\left(\begin{array}{c}K=R_{\mathrm{nn}}{b}_n^2-\mu \log \left(b_n^2\right)\\ k_{b_n}={2R}_{\mathrm{nn}}{b}_n-2\mu /b_n\\ L_{b_n}={2R}_{\mathrm{nn}}+2\mu /b_n^2\end{array}\right),$于是(12)和(13)分别简化为$\left(\begin{array}{c}{\mathrm{Qb}}_i+g_{y^{\prime }}+H_{y^{\prime }}{\mathrm{\Delta y}}^{\prime }=0\\ Q^T{y}^{\prime }+k_{b_i}+L_{b_i}{\mathrm{\Delta b}}_i=0\end{array}\right),$求解得：

$\left(\begin{array}{c}{\mathrm{\Delta y}}^{\prime }=-{\left(H_{y^{\prime }}\right)}^{-1}({\mathrm{Qb}}_i+g_{y^{\prime }})\\ {\mathrm{\Delta b}}_i=-(Q^T{y}^{\prime }+k_{b_i})/L_{b_i}\end{array}\right)---\left(14\right)$

步骤L、利用公式(14)迭代变量y'和b_i。如果得到B的秩为1，那么B可分 解为B＝bb^T，b就是求得的最优解；如果B的秩大于1，那么就把矩阵B的特征 值最大值对应的特征向量作为近似的最优解。

步骤M、利用公式和最优解b，可以得到原始信号的重建结 果。把这个结果作为Gerchberg-Saxton算法的初始值，即

步骤N、逐个优化X⁰中的值，即直到获得精 确的信号。

有益效果：本发明涉及的利用变换域相位信息重建原始信号的方法，采用离 散分数阶傅里叶变换，在变换域获得合适数目的相位信息。相位作为已知的信息 被用于原始信号的重建，相比现有的信号重建算法，本发明考虑了另外一类先验 条件，即不再是对原始信号特性的先验，而是增加变换矩阵的维数。本发明将原 始的重建问题转化为凸优化问题，并利用经典的凸优化理论和算法解决问题，这 里使用了块坐标下降算法和内点法结合的幅度恢复算法。相比现有信号重建算 法，如：Gerchberg-Saxton算法和贪婪算法，幅度恢复算法可以使用较少数目的 相位信息去重建原始信号。

附图说明

图1为本发明从变换域相位到原始信号重建的系统结构框图；

图2为具体实施方式中所述对线性调频信号重建的均方重建误差；

图3为具体实施方式中所述对线性调频信号重建的结果(离散分数阶傅里叶 变换的阶数等于1)；

图4为具体实施方式中所述的一维随机实数信号的重建结果(a＝1)；

图5为具体实施方式中所述的一维随机信号重建的均方重建误差(a＝1)；

图6为具体实施方式中所述的二维图像重建的结果(a＝1)，(1)(2)(3)(4)是相 位数目分别为32×32,64×32,96×32,128×32的重建结果，(5)(6)(7)(8)是原始图 像分别减去重建图像(1)(2)(3)(4)的误差结果；

图7为具体实施方式中所述的一维随机复数信号的重建结果(a＝1)。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本 发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发 明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，本发明从离散分数阶傅里叶变换域相位到原始信号重建的系统 结构框图，其中文字部分表示整个系统的流程，公式部分表示每个流程的数学表 达或问题的数学表示。采用上述流程进行信号重建时，按照以下各步骤：

步骤1、数据输入装置输入信号或图像数据；

步骤2、信号或图像离散分数阶傅里叶变换装置按照以下方法对输入的数据 进行离散分数阶傅里叶变换，图3是当离散分数阶傅里叶变换的阶数等于1时对 线性调频信号重建的结果；

步骤201、将输入的信号或图像数据转化为标准的列向量x，对于二维的图 像数据，可以逐列进行重建。图6是二维图像重建的结果(离散分数阶傅里叶变 换的阶数等于1)，(1)(2)(3)(4)是相位数目分别为32×32,64×32,96×32,128×32 的重建结果，(5)(6)(7)(8)是原始图像分别减去重建图像(1)(2)(3)(4)的误差结果；

步骤202、将大小为p的列向量作为输入向量进行离散分数阶傅里叶变换， 得到一个列向量或矩阵。实施方式中采用四种不同大小的变换矩阵(p×p、 2p×p、3p×p、4p×p)，这样就会得到四个不同大小的列向量(p、2p、3p、 4p)，或四个不同大小的矩阵(p×p、2p×p、3p×p、4p×p)。离散分数阶傅 里叶变换具体按照以下公式，y＝Ax。式中，x＝[x₁,x₂,...,x_p]^T为输入向量； y＝[y₁,y₂,...,y_n]^T为输出向量，n∈{p,2p,3p,4p}；“T”代表矩阵转置；A为 离散分数阶傅里叶变换矩阵，A的大小分为四种情况，即p×p、2p×p、 3p×p、4p×p；

步骤203、提取四种变换矩阵下的离散分数阶傅里叶变换域的相位信息， 就是用输出向量y除以它的幅度信息，表达式为

步骤204、将步骤203得到的相位信息u作为已知的列向量数据，并记下分 数阶傅里叶变换域的系数矩阵A。

步骤3、通过问题转化装置，利用步骤2得到的相位信息u和变换矩阵A计 算凸优化问题需要的矩阵M，它的计算公式为其 中“H”代表矩阵或向量的共轭转置，代表矩阵的伪逆。

步骤4、通过幅度恢复算法装置将问题分解为多个部分去求解优化，这里采 用了块坐标下降算法和内点法。

步骤401、块坐标下降算法是将凸优化问题分解为多个部分的优化问题

i^c∈{1,...,i-1,i+1,...,n},i＝{1,...,n}，

其中

$b_{i^c}={[b_1,b_2,...,b_{i-1,}{b}_{i+1},...,b_n]}^T,$

此时原始凸优化问题就被分解为n个子问题；

步骤402、解决每个子问题，采用经典的log-barrier内点算法。对步骤401 中问题的目标函数分别对和b_i求偏导，

$\frac{\partial F(b_{i^c},b_i)}{{\partial b}_{i^c}}=\mathrm{Re}{\left(M\right)}_{i^{c}i}^T{b}_i+\frac{2{\mathrm{\mu b}}_{i^c{i}^c}^{-1}{b}_{i^c}}{(1-b_{i^c}^T{b}_{i^c{i}^c}^{-1}{b}_{i^c})}=I\left(b_{i^c}\right)$

$\frac{\partial F(b_{i^c},b_i)}{{\partial b}_i}=\mathrm{Re}{\left(M\right)}_{i^{c}i}^T{b}_{i^c}+2\mathrm{Re}{\left(M\right)}_{\mathrm{ii}}{b}_i-\frac{2\mu }{b_i}=J\left(b_i\right)$

然后，二阶泰勒展开上面两个式子：

$I(b_{i^c}+\Delta b_{i^c})=I\left(b_{i^c}\right)+\frac{{2\mu \left(b_{i^c{i}^c}^{-1}\right)}^T}{(1-b_{i^c}^T{b}_{i^c{i}^c}^{-1}{b}_{i^c})}\Delta b_{i^c}\frac{{4\mu P^{-1}{b}_{i^c}\left(b_{i^c{i}^c}^{-1}{b}_{i^c}\right)}^T}{{(1-b_{i^c}^T{b}_{i^c{i}^c}^{-1}{b}_{i^c})}^2}\Delta b_{i^c}$

$J(b_i+\Delta b_i)=J\left(b_i\right)+2\mathrm{Re}{\left(M\right)}_{\mathrm{ii}}\Delta b_i+\frac{2\mu }{b_i^2}\Delta b_i$

令

$G=-\mu \log (1-b_{i^c}^T{b}_{i^c{i}^c}^{-1}{b}_{i^c})$

$g_{b_{i^c}}=\mathrm{dG}/{\mathrm{db}}_{i^c}=2\mu b_{i^c{i}^c}^{-1}{b}_{i^{c}i}/(1-b_{i^c}^T{b}_{i^c{i}^c}^{-1}{b}_{i^c})$

$H_{b_{i^c}}=d^{2}G/{\mathrm{db}}_{i^c}^2=2\mu {\left(b_{i^c{i}^c}^{-1}\right)}^T/(1-b_{i^c}^T{b}_{i^c{i}^c}^{-1}{b}_{i^c})+4\mu P^{-1}{b}_{i^c}{\left(b_{i^c{i}^c}^{-1}{b}_{i^c}\right)}^T/{(1-b_{i^c}^T{b}_{i^c{i}^c}^{-1}{b}_{i^c})}^2$

和$\left(\begin{array}{c}K=\mathrm{Re}{\left(M\right)}_{\mathrm{ii}}{b}_i^2-\mu \log \left(b_i^2\right)\\ k_{b_i}=2\mathrm{Re}{\left(M\right)}_{\mathrm{ii}}{b}_i-2\mu /b_i\\ L_{b_i}=2\mathrm{Re}{\left(M\right)}_{\mathrm{ii}}+2\mu /b_i^2\end{array}\right)\text{,}$而二阶泰勒展开式分别简化为 $\left(\begin{array}{c}\mathrm{Re}{\left(M\right)}_{i^{c}i}{b}_i+g_{b_{i^c}+H_{b_{i^c}}{\mathrm{\Delta b}}_{i^c}=0}\\ \mathrm{Re}{\left(M\right)}_{i^{c}i}^T+b_{i^{c}i}+k_{b_i}+L_{b_i}{\mathrm{\Delta b}}_i=0\end{array}\right).$进一步求解得变量和b_i的变化量和Δb_i： $\left(\begin{array}{c}{\mathrm{\Delta b}}_{i^c}=-{\left(H_{b_{i^v}}\right)}^{-1}({\mathrm{Re}\left(M\right)}_{i^{c}i}{b}_i+g_{g_{i^c}})\\ {\mathrm{\Delta b}}_i=-\frac{{\mathrm{Re}\left(M\right)}_{i^{c}i}^T{b}_{i^c}+k_{b_i}}{L_{b_i}}\end{array}\right).$设置步长为s，最终得到变量每次迭代后的结 果为：

步骤403、对迭代最后的矩阵B分析，如果得到B的秩为1，那么B可分解 为B＝bb^T，b就是求得的最优解；如果B的秩大于1，那么就把矩阵B的特征值 最大值对应的特征向量作为近似的最优解。

这是因为：

步骤5、利用公式和最优解b，可以得到原始信号的恢复结 果。把这个结果作为Gerchberg-Saxton算法的初始值，即再逐个优化X⁰中的值，即最终获得较精确的解。

为了验证本发明方法的效果，进行了以下实验：

1、实验条件：

在一台64位操作系统的计算机上进行验证实验，该计算机配置为Intel(R) Core(TM)i5-2400处理器(3100兆赫兹)和4000兆字节随机存取存储器(RAM)， 编程语言用的是Matlab(7.11版本)。

2、实验方法：

实验采用幅度恢复算法框架(如附图1所示)，将图中离散分数阶傅里叶变 换矩阵的大小改变。实验采用的输入信号分别是线性调频信号LFM(大小为 p＝160)，图3是对线性调频信号重建的结果；Cameraman图像(p×p＝32×32)， (图6是二维图像的重建结果)，随机一维实数信号(大小为p＝32)，图4是一维随 机实数信号的重建结果；随机一维复数信号(大小为p＝32)，(图7是复数信号的 重建结果)。改变分数阶傅里叶变换矩阵的大小，然后实施信号的重建：

重建过程：首先，对一维数据直接进行分数阶傅里叶变换(不同大小矩阵， 矩阵A的大小分别为p×p、2p×p、3p×p、4p×p)，将二维数据看成多个一维 数据的组合并进行离散分数阶傅里叶变换(不同大小矩阵，矩阵A的大小分别为 p×p、2p×p、3p×p、4p×p)。然后，提取变换域的相位信息作为已知条件， 并根据已知的相位信息和矩阵A求解相关的矩阵M，再由该矩阵利用块坐标下 降算法和内点法逐步优化变换域幅度矩阵B，B＝bb^T。接着，利用 Gerchberg-Saxton算法修正结果，得到精确的重建信号。

3、实验结果的评价指标：

实验结果采用均方重建误差(Normalized Mean Squared Reconstruction Error， NMSR)和误差方差(Error Variance，EV)。

均方重建误差指的是原始信号同重建信号的差的平方和再除以原始信号本 身的平方和。

误差方差指的是均方重建误差平方和的平均值。

均方重建误差的公式表达为：

$\mathrm{NMSR}=\frac{{\left|\right|x-\hat{x}\left|\right|}_2^2}{{\left|\right|x\left|\right|}_2^2},-\log \left(\mathrm{NMSR}\right)=-\log \left(\frac{{\left|\right|x-\hat{x}\left|\right|}_2^2}{{\left|\right|x\left|\right|}_2^2}\right)$

误差方差的公式表达为：

$\mathrm{EV}=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{\left({\mathrm{NMSR}}_i\right)}^2$

其中x是原始信号，是重建信号，NMSR_i是第i次重建的均方误差，n是 重建次数。在同一信号和相同的迭代次数下，NMSR值越小，就代表重建的精确 度越高，效果越好。通常，采用NMSR对数的负数来显示，值越大，就代表重 建的精确度越高，效果越好。随着重建次数的增加，若误差方差基本不变，则表 明算法稳定。图2是线性调频信号重建的均方重建误差，图5是一维随机信号重 建的均方重建误差。

4、与现有技术的对比实验结果：

表1给出了一维随机实数信号分别采用Gerchberg-Saxton算法、贪 婪算法和幅度恢复算法进行信号重建的结果。测试结果同时给出了不同大小变换 矩阵下的均方重建误差、误差方差和运行时间。表2分别给出了一维随机实数信 号和一维随机复数信号在不同大小变换矩阵下的运行时间和均方重建误差。这里 仅考虑了当离散分数阶傅里叶变换的阶数等于1，即a＝1的情况，也就是说相位 信息是从原始信号进行标准离散傅里叶变换后获得的。

表1.不同算法的一维随机实数信号的信号重建结果比较

表2.一维随机实数信号和一维随机复数信号的重建结果比较

由表1和表2可以知道，随着变换矩阵的增大，贪婪算法、Gerchberg-Saxton 算法和幅度恢复算法都可以重建原始信号。同贪婪算法和Gerchberg-Saxton相比， 幅度恢复算法需要较少数目的相位信息就可以较好的重建原始信号。在相同的均 方重建误差条件下，幅度恢复算法需要至少64个相位信息，而另外两种算法需 要至少96个相位信息。从运行时间来看，相同大小的变换矩阵下，幅度恢复算 法比其他两种运行的时间长。尽管如此，可以发现当已知的相位信息数目同原始 信号数目相同的条件下，幅度恢复算法仍然可以重建信号。另外，贪婪算法和 Gerchberg-Saxton算法的误差方差有变化，幅度恢复算法基本不变，因此幅度恢 复算法较为稳定。对于复数信号，当采集的相位信息个数大于等于96时可以较 好的重建原始信号。