Algorithmic determination of wave cycle parts
位相回復とは、 位相問題 に対する解を アルゴリズム的に 求めるプロセスです。 振幅 、位相の 複素スペクトルが与えられた場合、以下のようになります 。 F ( k ) {\displaystyle F(k)} | F ( k ) | {\displaystyle |F(k)|} ψ ( k ) {\displaystyle \psi (k)}
F ( k ) = | F ( k ) | e i ψ ( k ) = ∫ − ∞ ∞ f ( x ) e − 2 π i k ⋅ x d x {\displaystyle F(k)=|F(k)|e^{i\psi (k)}=\int _{-\infty }^{\infty }f(x)\ e^{-2\pi ik\cdot x}\,dx} ここで、 x は M 次元空間座標、 kは M 次元 空間周波数 座標である 。位相回復とは、測定された振幅に対して一連の制約を満たす位相を見つけることである。位相回復の重要な応用分野には、 X線結晶構造解析 、 透過型電子顕微鏡 、 コヒーレント回折イメージング などがあり、 これらについては [1]、 位相回復問題の1次元および2次元の両方のケースに対する一意性定理(位相のない1次元 逆散乱問題 を含む)は、Klibanovと彼の協力者によって証明された(参考文献を参照)。 M = 2 {\displaystyle M=2}
ここでは1次元離散フーリエ変換 (DFT)位相回復問題を考える 。複素信号のDFTは 次のように表される。 f [ n ] {\displaystyle f[n]}
F [ k ] = ∑ n = 0 N − 1 f [ n ] e − j 2 π k n N , = | F [ k ] | ⋅ e j ψ [ k ] k = 0 , 1 , … , N − 1 {\displaystyle F[k]=\sum _{n=0}^{N-1}f[n]e^{-j2\pi {\frac {kn}{N}},}=|F[k]|\cdot e^{j\psi [k]}\quad k=0,1,\ldots ,N-1} 、
のオーバーサンプリングDFTは 次のように与えられる。 x {\displaystyle x}
F [ k ] = ∑ n = 0 N − 1 f [ n ] e − j 2 π k n M , k = 0 , 1 , … , M − 1 {\displaystyle F[k]=\sum _{n=0}^{N-1}f[n]e^{-j2\pi {\frac {kn}{M}},}\quad k=0,1,\ldots ,M-1} 、
どこ 。 M > N {\displaystyle M>N}
DFT演算子は全単射であるため、これは位相を復元することと等価です 。信号の復元は、フーリエ振幅ではなく自己相関系列から行うのが一般的です。つまり、ゼロ パディング後の ベクトルを と表記します 。 の自己相関系列 は次のように定義されます。 ψ [ k ] {\displaystyle \psi [k]} f ^ {\displaystyle {\hat {f}}} f {\displaystyle f} N − 1 {\displaystyle N-1} f ^ {\displaystyle {\hat {f}}}
g [ m ] = ∑ i = max { 1 , m + 1 } N f ^ i f ^ i − m ¯ , m = − ( N − 1 ) , … , N − 1 {\displaystyle g[m]=\sum _{i=\max\{1,m+1\}}^{N}{\hat {f}}_{i}{\overline {{\hat {f}}_{i-m}}},\quad m=-(N-1),\ldots ,N-1} 、
の DFT は と 表記され 、 を満たします 。 g [ m ] {\displaystyle g[m]} G [ k ] {\displaystyle G[k]} G [ k ] = | F [ k ] | 2 {\displaystyle G[k]=|F[k]|^{2}}
方法
エラー削減アルゴリズム 誤差低減は、 Gerchberg–Saxtonアルゴリズム の一般化です。4段階のプロセスを反復することで、 の測定値から を解きます。 回の反復における手順は以下のとおりです。 f ( x ) {\displaystyle f(x)} | F ( u ) | {\displaystyle |F(u)|} k {\displaystyle k}
ステップ(1): 、、 は それぞれ 、、、 の推定値です 。最初のステップでは、 の フーリエ変換を 計算します。 G k ( u ) {\displaystyle G_{k}(u)} ϕ k {\displaystyle \phi _{k}} g k ( x ) {\displaystyle g_{k}(x)} F ( u ) {\displaystyle F(u)} ψ {\displaystyle \psi } f ( x ) {\displaystyle f(x)} g k ( x ) {\displaystyle g_{k}(x)}
G k ( u ) = | G k ( u ) | e i ϕ k ( u ) = F ( g k ( x ) ) . {\displaystyle G_{k}(u)=|G_{k}(u)|e^{i\phi _{k}(u)}={\mathcal {F}}(g_{k}(x)).} ステップ(2):信号方程式 [ 説明が必要 ] を介して回折パターンから計算されたの実験値が に代入され 、フーリエ変換の推定値が得られます。 | F ( u ) | {\displaystyle |F(u)|} | G k ( u ) | {\displaystyle |G_{k}(u)|}
G k ′ ( u ) = | F ( u ) | e i ϕ k ( u ) , {\displaystyle G'_{k}(u)=|F(u)|e^{i\phi _{k}(u)},} ここで、' は後で破棄される中間結果を示します。
ステップ(3):フーリエ変換の推定値 は逆フーリエ変換される。 G k ′ ( u ) {\displaystyle G'_{k}(u)}
g k ′ ( x ) = | g k ′ ( x ) | e i θ k ′ ( x ) = F − 1 ( G k ′ ( u ) ) . {\displaystyle g'_{k}(x)=|g'_{k}(x)|e^{i\theta '_{k}(x)}={\mathcal {F}}^{-1}(G'_{k}(u)).} ステップ(4): オブジェクトの新しい推定値、 がオブジェクトの制約を満たすように変更する必要があります [ 説明が必要 ] 。 したがって、は次のように区分的に定義されます。 g k ′ ( x ) {\displaystyle g'_{k}(x)} g k + 1 ( x ) {\displaystyle g_{k+1}(x)} g k + 1 ( x ) {\displaystyle g_{k+1}(x)}
g k + 1 ( x ) ≡ { g k ′ ( x ) , x ∉ γ , 0 , x ∈ γ , {\displaystyle g_{k+1}(x)\equiv {\begin{cases}g'_{k}(x),&x\notin \gamma ,\\0,&x\in \gamma ,\end{cases}}} ここで、は オブジェクト制約を満たさない 領域です。新たな推定値 が得られ、4段階のプロセスが繰り返されます。 γ {\displaystyle \gamma } g k ′ ( x ) {\displaystyle g'_{k}(x)} g k + 1 ( x ) {\displaystyle g_{k+1}(x)}
このプロセスは、フーリエ制約と物体制約の両方が満たされるまで続けられます。理論的には、このプロセスは常に 収束に至ります [1]。 しかし 、満足のいく画像を生成するために必要な反復回数が非常に多いため(通常2000回以上)、誤差低減アルゴリズム単体では実用には適していません。
ハイブリッド入出力アルゴリズムは、誤差低減アルゴリズムの修正版です。最初の3つの段階は同一です。ただし、 は の 推定値としてではなく、 の推定値である 出力関数 に対応する入力関数として機能します 。 [1] 4番目のステップでは、関数が オブジェクトの制約に違反する場合、 の値は ゼロに近づきますが、最適にはゼロにはなりません。ハイブリッド入出力アルゴリズムの主な利点は、関数に以前の反復に関する フィードバック情報 が含まれているため 、停滞の可能性が低減されることです。ハイブリッド入出力アルゴリズムは、誤差低減アルゴリズムよりも大幅に速く解に収束することが示されている。その収束率は、ステップサイズ最適化アルゴリズムによってさらに向上させることができます。 [2] g k ( x ) {\displaystyle g_{k}(x)} f ( x ) {\displaystyle f(x)} g k ′ ( x ) {\displaystyle g'_{k}(x)} f ( x ) {\displaystyle f(x)} g k ′ ( x ) {\displaystyle g'_{k}(x)} g k + 1 ( x ) {\displaystyle g_{k+1}(x)} g k ( x ) {\displaystyle g_{k}(x)}
g k + 1 ( x ) ≡ { g k ′ ( x ) , x ∉ γ , g k ( x ) − β g k ′ ( x ) , x ∈ γ . {\displaystyle g_{k+1}(x)\equiv {\begin{cases}g'_{k}(x),&x\notin \gamma ,\\g_{k}(x)-\beta {g'_{k}(x)},&x\in \gamma .\end{cases}}} これは 0から1までの値を取るフィードバックパラメータです。ほとんどのアプリケーションでは、 最適な結果が得られます。{ Scientific Reports volume 8, Article number: 6436 (2018)} β {\displaystyle \beta } β ≈ 0.9 {\displaystyle \beta \approx 0.9}
シュリンクラップ 2次元位相回復問題では、 とその共役が同一のフーリエ係数を持つという解の退化が存在する 。 これ は「画像ツイン化」につながり、位相回復アルゴリズムが停滞し、物体とその共役の両方の特徴を持つ画像が生成される 。 [ 3 ] シュリンクラップ法は、物体振幅の現在の推定値にローパスフィルタリング( ガウス 分布との畳み込みによる)を施し、閾値を適用することで、サポートの推定値を定期的に更新し、画像の曖昧性を低減する。 [4] f ( x ) {\displaystyle f(x)} f ∗ ( − x ) {\displaystyle f^{*}(-x)}
位相回復は不適切問題です。基となる信号を一意に識別するには、 Gerchberg–Saxtonアルゴリズム のような事前情報を追加する手法に加えて、短時間フーリエ変換(STFT)のような振幅のみの測定を追加する方法もあります。
以下に紹介する方法は主にJaganathanら [5] の研究に基づいています 。
からサンプリングされた 離散信号が与えられた場合、長さ W のウィンドウを使用して の STFT を計算します 。 これは と表記されます。 x = ( f [ 0 ] , f [ 1 ] , . . . , f [ N − 1 ] ) T {\displaystyle \mathbf {x} =(f[0],f[1],...,f[N-1])^{T}} f ( x ) {\displaystyle f(x)} w = ( w [ 0 ] , w [ 1 ] , . . . , w [ W − 1 ] ) T {\displaystyle \mathbf {w} =(w[0],w[1],...,w[W-1])^{T}} f {\displaystyle \mathrm {f} } Y {\displaystyle \mathbf {Y} }
Y [ m , r ] = ∑ n = 0 N − 1 f [ n ] w [ r L − n ] e − i 2 π m n N {\displaystyle Y[m,r]=\sum _{n=0}^{N-1}{f[n]w[rL-n]e^{-i2\pi {\frac {mn}{N}}}}}
および の 場合 、パラメータ は 隣接する短時間セクション間の時間的な間隔を示し、パラメータ は 考慮される短時間セクションの数を示します。 0 ≤ m ≤ N − 1 {\displaystyle 0\leq m\leq N-1} 0 ≤ r ≤ R − 1 {\displaystyle 0\leq r\leq R-1} L {\displaystyle L} R = ⌈ N + W − 1 L ⌉ {\displaystyle R=\left\lceil {\frac {N+W-1}{L}}\right\rceil }
STFTのもう一つの解釈(スライディングウィンドウ解釈と呼ばれる)は、離散フーリエ変換(DFT)を用いて行うことができます。ここでは、 ウィンドウ要素をシフトおよび反転したウィンドウ要素とします 。すると、 w r [ n ] = w [ r L − n ] {\displaystyle w_{r}[n]=w[rL-n]} w {\displaystyle \mathbf {w} }
Y = [ Y 0 , Y 1 , . . . , Y R − 1 ] {\displaystyle \mathbf {Y} =[\mathbf {Y} _{0},\mathbf {Y} _{1},...,\mathbf {Y} _{R-1}]} 、 どこ 。 Y r = x ∘ w r {\displaystyle \mathbf {Y} _{r}=\mathbf {x} \circ \mathbf {w} _{r}}
問題の定義 を の STFT の振幅の 2 乗に対応する測定値 と し 、 を対角要素を持つ 対角行列 とすると 、STFT 位相回復は次のように表すことができます。 Z w [ m , r ] = | Y [ m , r ] | 2 {\displaystyle {Z}_{w}[m,r]=|Y[m,r]|^{2}} N × R {\displaystyle N\times R} x {\displaystyle \mathbf {x} } W r {\displaystyle \mathbf {W} _{r}} N × N {\displaystyle N\times N} ( w r [ 0 ] , w r [ 1 ] , … , w r [ N − 1 ] ) . {\displaystyle \left(w_{r}[0],w_{r}[1],\ldots ,w_{r}[N-1]\right).}
および に対して となる値 を見つけます 。ここでは 点逆 DFT 行列 の 番目の列 です 。 x {\displaystyle \mathbf {x} } Z w [ m , r ] = | ⟨ f m , W r x ⟩ | 2 {\displaystyle Z_{w}[m,r]=\left|\left\langle \mathbf {f} _{m},\mathbf {W} _{r}\mathbf {x} \right\rangle \right|^{2}} 0 ≤ m ≤ N − 1 {\displaystyle 0\leq m\leq N-1} 0 ≤ r ≤ R − 1 {\displaystyle 0\leq r\leq R-1} f m {\displaystyle \mathbf {f} _{m}} m {\displaystyle m} N {\displaystyle N}
直感的には、 とともに増大する計算複雑性により、 この方法は実用的ではないと考えられます。しかし実際には、 を 満たす任意の パラメータについて、に対応する測定値のみを考慮すれば十分です 。 N {\displaystyle N} 0 ≤ m ≤ M {\displaystyle 0\leq m\leq M} M {\displaystyle M} 2 W ≤ M ≤ N {\displaystyle 2W\leq M\leq N}
より具体的には、信号とウィンドウの両方が 消失して いない場合、つまり、 すべての に対して であり 、 すべての に対して 、次の要件が満たされていれば、信号は その STFT 振幅から一意に識別できます。 x [ n ] ≠ 0 {\displaystyle x[n]\neq 0} 0 ≤ n ≤ N − 1 {\displaystyle 0\leq n\leq N-1} w [ n ] ≠ 0 {\displaystyle w[n]\neq 0} 0 ≤ {\displaystyle 0\leq } n ≤ W − 1 {\displaystyle n\leq W-1} x {\displaystyle \mathbf {x} }
L < W ≤ N / 2 {\displaystyle L<W\leq N/2} 、 2 W ≤ M ≤ N {\displaystyle 2W\leq M\leq N} 。 証明はジャガナサンの研究 [5] に記載されており、 STFT位相回復を次の最小二乗問題として再定式化している。
min x ∑ r = 0 R − 1 ∑ m = 0 N − 1 ( Z w [ m , r ] − | ⟨ f m , W r x ⟩ | 2 ) 2 {\displaystyle \min _{\mathbf {x} }\sum _{r=0}^{R-1}\sum _{m=0}^{N-1}\left(Z_{w}[m,r]-\left|\left\langle \mathbf {f} _{m},\mathbf {W} _{r}\mathbf {x} \right\rangle \right|^{2}\right)^{2}} 。
このアルゴリズムは、理論的な回復保証はありませんが、隣接する短時間セクション間に大きな重複がある場合、経験的にグローバル最小値に収束することができます。
半正定値緩和ベースのアルゴリズム 回復保証を確立する方法の一つは、問題を半正定値計画(SDP)として定式化することです。これは、変換を用いて問題を高次元空間に埋め込み 、ランク1制約を緩和して凸計画を得る方法です。再定式化された問題は以下のように記述されます。 X = x x ∗ {\displaystyle \mathbf {X} =\mathbf {x} \mathbf {x} ^{\ast }}
次 を解いて 求めます 。 X ^ {\displaystyle \mathbf {\hat {X}} } m i n i m i z e t r a c e ( X ) s u b j e c t t o Z [ m , r ] = t r a c e ( W r ∗ f m f m ∗ W r X ) X ⪰ 0 {\displaystyle {\begin{aligned}&\mathrm {minimize} ~~~\mathrm {trace} (\mathbf {X} )\\[6pt]&\mathrm {subject~to} ~~Z[m,r]=\mathrm {trace} (\mathbf {W} _{r}^{\ast }\mathbf {f} _{m}\mathbf {f} _{m}^{\ast }\mathbf {W} _{r}\mathbf {X} )\\[0pt]&~~~~~~~~~~~~~~~~~~~\mathbf {X} \succeq 0\end{aligned}}} 1 ≤ m ≤ M {\displaystyle 1\leq m\leq M} 0 ≤ r ≤ R − 1 {\displaystyle 0\leq r\leq R-1}
が見つかったら、 ランク 1 の近似値を使用して 信号を回復できます。 X ^ {\displaystyle \mathbf {\hat {X}} } x {\displaystyle \mathbf {x} }
アプリケーション 位相回復は、コヒーレント回折イメージング (CDI)の重要な要素です。CDIでは、対象物から散乱された 回折 パターンの強度を測定します。次に、位相回復アルゴリズムを用いて回折パターンの位相を取得し、対象物の画像を構築します。このように、位相回復により、 光学レンズを 必要とせずに回折パターンを画像に変換することができます 。
位相回復アルゴリズムを用いることで、複雑な光学系とその収差を特徴づけることが可能です。 [6] 例えば、位相回復は ハッブル宇宙望遠鏡 の欠陥のある光学系の診断と修復に使用されました。 [7] [8]
位相回復の他の応用としては、 X線結晶構造解析 [9] や 透過型電子顕微鏡法 などがある。
参照
参考文献 ^ abc Fienup, JR (1982-08-01). 「位相回復アルゴリズム:比較」 . 応用光学 . 21 (15): 2758–69 . Bibcode :1982ApOpt..21.2758F. doi :10.1364/AO.21.002758. ISSN 0003-6935. PMID 20396114. ^ Marchesini, S. (2007年1月25日). 「招待論文:位相回復のための反復射影アルゴリズムの統一的評価」 Review of Scientific Instruments . 78 (1): 011301–011301–10. arXiv : physics/0603201 . Bibcode :2007RScI...78a1301M. doi :10.1063/1.2403783. ISSN 0034-6748. PMID 17503899. S2CID 7462041. ^ Fienup, JR; Wackerman, CC (1986年11月1日). 「位相回復停滞の問題と解決策」 . Journal of the Optical Society of America A. 3 ( 11): 1897. Bibcode :1986JOSAA...3.1897F. doi :10.1364/JOSAA.3.001897. ISSN 1084-7529. ^ Marchesini, S.; He, H.; Chapman, HN; Hau-Riege, SP; Noy, A.; Howells, MR; Weierstall, U.; Spence, JCH (2003-10-28). 「回折パターンのみからのX線画像再構成」. Physical Review B. 68 ( 14) 140101. arXiv : physics/0306174 . Bibcode :2003PhRvB..68n0101M. doi :10.1103/PhysRevB.68.140101. ISSN 0163-1829. S2CID 14224319. ^ ab Jaganathan, Kishore; Eldar, Yonina C.; Hassibi, Babak (2016年6月). 「STFT位相検索:一意性保証と回復アルゴリズム」. IEEE Journal of Selected Topics in Signal Processing . 10 (4): 770– 781. arXiv : 1508.02820 . Bibcode :2016ISTSP..10..770J. doi : 10.1109/JSTSP.2016.2549507 . ISSN 1941-0484. ^ Fienup, JR (1993-04-01). 「複雑な光学系のための位相回復アルゴリズム」 . 応用光学 . 32 (10): 1737– 1746. Bibcode :1993ApOpt..32.1737F. doi :10.1364/AO.32.001737. ISSN 2155-3165. PMID 20820307. ^ 「一人称:科学者の発見が宇宙に焦点を当てる」 www.wbur.org 2022年4月 2022年 5月30日 閲覧 。 ロバート・ゴンサルベス教授へのインタビュー。 ^ Krist, JE; Burrows, CJ (1995-08-01). 「ハッブル宇宙望遠鏡の修理前と修理後の画像の位相回復解析」. 応用光学 . 34 (22): 4951–64 . Bibcode :1995ApOpt..34.4951K. doi :10.1364/AO.34.004951. PMID 21052338. ^ Millane, Rick P.; Arnal, Romain D. (2015). 「高分子結晶学的位相問題の一意性」 . Acta Crystallographica Section A: Foundations and Advances . 71 (6): 592– 598. doi :10.1107/S2053273315015387. PMID 26522408. Klibanov, MV (1985). 「コンパクトに支えられた関数のフーリエ変換の絶対値からの決定の一意性について」 『ソビエト数学 - Doklady 』 32 : 668–670 . Klibanov, MV (1987). 「コンパクトな台を持つ関数のフーリエ変換の絶対値からの決定と逆散乱問題」. 微分方程式 . 22 : 1232–1240 . Klibanov, MV (1987). 「逆散乱問題とフーリエ変換の絶対値からの関数の復元」. Siberian Math. J. 27 ( 5): 708– 719. doi :10.1007/bf00969199. S2CID 120840929. Klibanov, MV (1989). 「連続力学モデルにおけるX線回折による結晶格子の歪みの決定の一意性」. 微分方程式 . 25 : 520–527 . Klibanov, MV & Sacks, PE (1992). 「位相のない逆散乱と光学における位相問題」. J. Math. Phys . 33 (11): 2813– 3821. Bibcode :1992JMP....33.3813K. doi :10.1063/1.529990. Klibanov, MV; Sacks, PE (1994). 「逆散乱の位相問題におけるポテンシャルの部分知識の利用」 J. Comput. Phys . 112 (2): 273– 281. Bibcode :1994JCoPh.112..273K. doi :10.1006/jcph.1994.1099. Klibanov, MV; Sacks, PE; Tikhonravov, AV (1995). 「位相回復問題」. 逆問題 . 11 (1): 1– 28. Bibcode :1995InvPr..11....1K. doi :10.1088/0266-5611/11/1/001. S2CID 250916850. Klibanov, MV (2006). 「フーリエ変換の絶対値からの2次元関数の復元について」 J. Math. Anal. Appl . 323 (2): 818– 843. doi : 10.1016/j.jmaa.2005.10.079 .