Spitz補間
6507 ワード
参考文献(リンク)Spitz S.Seismic trace interpolation in the F-X domain[J]. Geophysics, 1991, 56(6):785-794.
Spitz補間は最も古典的な補間法であり,前のブログで述べた線形同相軸を持つ地震データが周波数‐空間領域で線形回帰できる理論を利用した.(回顧:線形回帰フィルタの長さは同相軸(異なる傾き)の数に1を加え、実際の応用ではより長い自己回帰フィルタが採用され、アルゴリズムをより安定させることができる.自己回帰の利点は、傾き情報を知る必要がなく、地震データだけを自分に使うことです)
Spitz補間の構想は,まず低周波数無擬似周波数部分を用いて自己回帰係数(回帰)を計算し,チャネル暗号化後の隣接チャネル間の位相シフトを元の半分にすることで,暗号化後のデータの自己回帰係数を導出し,次いで自己回帰関係を用いて既知のデータから未知のデータ(予測)を導出することである.回顧と予測はいずれも最小二乗によって実現されるので,アルゴリズムは比較的安定である.欠点は:1、線形同相軸にしか使用できない;2.ルール・グリッドでのみ使用できます.
以下の記号はSpitz論文と一致している.
L L個の同相軸とN N個のNチャネルを含む地震データは周波数領域で以下のように表すことができる.
g k ( f ) = ∑ j = 1 L a j ( f ) z j k − 1 ( f ) , k = 1 , . . . , N g_k(f)=\sum_{j=1}^La_j(f)z_j^{k-1}(f),k=1,...,N gk(f)=j=1∑Laj(f)zjk−1(f),k=1,...,N
ここでg k(f)g_k(f)gk(f)はk番目のkチャネルデータを表し、a j(f)a_j(f)aj(f)は第j j j個の同相軸に対応するサブ波スペクトル(第1パスデータ)、z j(f)=e x p(2πi f p j)z_j(f)=exp(2pi i f p_j)zj(f)=exp(2πifpj)は、第j j jの同相軸隣接チャネル間の位相シフトpj p_に対応するj pj(傾きに対応).前のブログで紹介したCanalesメソッドによれば,前−後の単一ステップ予測フィルタリングは,g k(f)=Σj=1 L P j(f)g k−j(f),k=L+1と書くことができる.N g_k(f)=\sum_{j=1}^L P_j(f)g_{k-j}(f),k=L+1,...,N gk(f)=j=1∑LPj(f)gk−j(f),k=L+1,...,N g k ∗ ( f ) = ∑ j = 1 L P j ( f ) g k + j ∗ ( f ) , k = 1 , . . . , N − L g_k^*(f)=\sum_{j=1}^L P_j(f)g_{k+j}^*(f),k=1,...,N-L gk∗(f)=j=1∑LPj(f)gk+j∗(f),k=1,...,N−LのうちP j(f)P_j(f)Pj(f)は予測フィルタの係数である.(前方-後方フィルタリングの共役関係は付録から得ることができる)
Spitz法は従来の地震路の2つの中間に新しい地震路を挿入するので,補間後の地震データはg k′(f)=Σj=1 L aj(f)z′j k−1(f),k=1,.,を満たす.2 N − 1 g_k'(f)=\sum_{j=1}^La_j(f){z'}_j^{k-1}(f),k=1,...,2N-1 gk′(f)=j=1∑Laj(f)z′jk−1(f),k=1,...,2 N−1は、g k′(f)=Σj=1 LPj′(f)g′k−j(f)、k=L+1の順方向から逆方向への単一ステップ予測フィルタ形式を再び書くことができる.2 N − 1 g'_k(f)=\sum_{j=1}^L P'_j(f){g'}_{k-j}(f),k=L+1,...,2N-1 gk′(f)=j=1∑LPj′(f)g′k−j(f),k=L+1,...,2N−1 g k ′ ( f ) = ∑ j = 1 L P ′ j ∗ ( f ) g ′ k + j ( f ) , k = 1 , . . . , 2 N − L − 1 g'_k(f)=\sum_{j=1}^L {P'}^*_j(f){g'}_{k+j}(f),k=1,...,2N-L-1 gk′(f)=j=1∑LP′j∗(f)g′k+j(f),k=1,...,2N−L−1
上記の方程式はg j’g'_についてj gj′の線形方程式群、係数はPj′P′_j Pj′は既知の値と未知の値をそれぞれ1つに配置し、整理結果:A(P′(f))[g 2′(f)g 4′(f)..g 2 N−2′(f)]T=T=B(P′(f))[g 1′(f)g 3′(f)..g 2 N−1′(f)]T A(P′(f)[g′_2(f)]T(g′_2(f)text}g′_4(f)text{}…text{}g′_{{2 N−2 N−2}(f)^T=B(P′(f)[g′_1(f)]text{xt{xt{xt}...text{}g′}g′{g′′{2 N−g′{2 N−2 N−}g'_3(f)text{}…text{}g'{2 N-1}(f)]^T A(P’(f)[g 2’(f)g 4’(f)…g 2 N−2’(f)]T=B(P′(f))[g1′(f) g3′(f) ... g2N−1′(f)]T
P’P’P’P’が既知であれば、上記式は2(2 N−L−1)2(2 N−L−1)2(2 N−L−1)個の方程式(前方後方フィルタ等式の個数)とN−1 N−1 N−1個の未知数(文献では2 N−1 2 N−1 2 N−1と表記されているが、書き間違えたのか、すべてのパスを含むのか)とを含み、最小二乗法(既知から未知へのフィルタに等価)で解くことができる.
( A + A ) − 1 A + B (A^+A)^{-1}A^+B (A+A)−1A+B
ここで、A+A^+A+は、A A A Aの複素共役転置である.
実際にはP’P’P’は未知であるが,P P Pで導くことができる.注意出力(既知の点を含む)と入力断面との位相シフト関係は、z j’(f)=z j(f/2)z'j(f)=z_j(f/2)zj’(f)=zj(f/2)(注:zj’(f)=e x p(2πi fpj/2)=zj(f/2)z'j(f)=exp(2\pi i f p_j/2)=z_j(f/2)zj’(f)=exp(2πifpj/2)=zj(f/2)であり、補間後の隣接チャネル間の位相シフトは元の半分となる.
P P Pとz z zの関係から得られる:
P j ′ ( f ) = P j ( f/2 ) P'_j(f)=P_j(f/2) Pj′(f)=Pj(f/2)
擬似周波数のない低周波数部分のみが利用されていることがわかる.
以上の導出はいくつかの詳細から与えられず、続く.
Spitzは別の方法でCanales文の導出を行った.付導出過程は、地震データx T(f)=[x 1(f),.,x N(f)]x^T(f)=[x_1(f),...,x_N(f)]xT(f)=[x 1(f),...,xN(f)]に対して、x(f)=g(f)+n(f)x(f)=g(f)+n(f)x(f)=g(f)+n(f)x(f)=g(f)+n(f)とする.定義z変換S j T(f)=[1,z j(f),.,z j N−1(f)]S_j^T(f)=[1,z_j(f),…,z_j^{N−1}(f)]SjT(f)=[1,zj(f),…,zjN−1(f)],ここで,z j(f)=e x p(2πi f p i)z_j(f)=exp(2pi ifp_i)zj(f)=exp(2πifpi)の場合:g(f)=Σj=1 L aj(f)S j(f)g(f)=sum_{j=1}^La_j(f)S_j(f)g(f)=j=1ΣL aj(f)Sj(f)ここでa j(f)a_j(f)aj(f)は、j番目jの同方向軸に対応するサブ波のフーリエ変換である.分けて、g 1=Σj=1 L a j⋅1 g 2=Σj=1 L a j⋅z j⋯g N=Σj=1 L a j⋅z j−1begin{aligned}g_と書く1&=&\sum_{j=1}^La_j\cdot 1\\g_2&=&\sum_{j=1}^La_j\cdot z_j\\&\cdots&\\g_N&=&\sum_{j=1}^La_j\cdot z_j^{N−1}end{aligned}g 1 g 2 gN=⋯=j=1ΣL aj⋅1 j=1ΣL aj⋅zj=1ΣL aj⋅zj=1ΣL aj⋅zjN−1行列として書く形式は、(3)[g 1 g 2⋮gn]=[1⋯1 z 1 z 2⋯z L⋮z 1 N−1 z 2 N−1⋯zL−1 N−1]× [ a 1 a 2 ⋮ a n ]\left[\begin{matrix} g_1\\g_2\\\vdots\\g_n\end{matrix}\right]=\left[\begin{matrix} 1 & 1&\cdots & 1\\z_1 & z_2&\cdots & z_L\\\vdots\\z_1^{N-1} & z_2^{N-1}&\cdots & z_L^{N-1}\end{matrix}\right]\times\left[\begin{matrix} a_1\\a_2\\\vdots\\a_n\end{matrix}\right]\tag{3} ⎣⎢⎢⎢⎡g1g2⋮gn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡1z1⋮z1N−11z2z2N−1⋯⋯⋯1zLzLN−1⎦⎥⎥⎥⎤×⎣⎢⎢⎡a 1 a 2⋮an⎥⎥⎤(3)は、マトリクスをS S S、すなわちS=[S 1,⋯ ,S L]S=[S_1,cdots,S_L]S=[S 1,⋯,SL]とする.(3)を行列形式と記す:x=S a+n x=Sa+n x=Sa+nがS S S Sとして形成された行列をファンデモン行列,S j S_S_j Sj間は線形無関係であり、S S SのランクはL L L(N>=L N>=L N>=L N>=L)である.S S S Sのいずれの行も、他のL L行の線形結合として表すことができる.例えば、L＋1 L＋1 L＋1行目は、前L L行の線形の組み合わせとして表すことができる.(4)(z 1 L,⋯ ,z L)=(P L,⋯ ,P 1)と書く× [ 1 1 ⋯ 1 z 1 z 2 ⋯ z L ⋮ z 1 L − 1 z 2 L − 1 ⋯ z L L − 1 ] (z_1^L,\cdots,z_L^L)=(P_L,\cdots,P_1)\times\left[\begin{matrix} 1 & 1&\cdots & 1\\z_1 & z_2&\cdots & z_L\\\vdots\\z_1^{L-1} & z_2^{L-1}&\cdots & z_L^{L-1}\end{matrix}\right]\tag{4} (z1L,⋯,zLL)=(PL,⋯,P1)×⎣⎢⎢⎡1 z 1⋮z 1 L−1 z 2 z 2 L−1⋯⋯⋯⋯⋯⋯1 zL zLL−1⎦⎥⎤(4)上式右側を左側に移動し、z j z_について整理するj zjの方程式は,任意のz j z_に対してj zj形式は,いずれも,Z L=z L−P 1 z L−1−⋯−P L=0 Z_である.L=z^L-P_1z^{L-1}-\cdots-P_L=0 ZL=zL−P 1 zL−1−⋯−PL=0はまたz j z_j zjはそれぞれ異なり、上記方程式にはL L個の解(零点)があるので、z j z_j zjは上記方程式のKK個の零点であり、Z L=(z−z 1)(z−z 2)⋯(z−z L)Z_と書くことができるL=(z−z_1)(z−z_2)cdots(z−z_L)ZL=(z−z 1)(z−z 2)⋯(z−zL)上記2式と比較してP 1=z 1+⋯+zLP 2=−(z 1 z 2+⋯+zL−1 zL)⋯(−1)L+1 z 1 z 1 z 2⋯+zL+1 zL)P L=(−1)L+1 z 1 z 2⋯zLbegin{aligned}P_1&=&z_1+\cdots+z_L\\P_2&=&-(z_1z_2+\cdots+z_{L-1}z_L)\\&\cdots&\\P_L&=&(-1)^{L+1}z_1z_2\cdots z_Lend{aligned}P 1 P 2 PL=⋯=z 1+⋯+zL−(z 1 z 2+⋯+zL−1 zL)(−1)L+1 z 1 z 2⋯+zL上式は,我々がこれまで推さなかった一般解である.
コードの説明:
概要
Spitz補間は最も古典的な補間法であり,前のブログで述べた線形同相軸を持つ地震データが周波数‐空間領域で線形回帰できる理論を利用した.(回顧:線形回帰フィルタの長さは同相軸(異なる傾き)の数に1を加え、実際の応用ではより長い自己回帰フィルタが採用され、アルゴリズムをより安定させることができる.自己回帰の利点は、傾き情報を知る必要がなく、地震データだけを自分に使うことです)
Spitz補間の構想は,まず低周波数無擬似周波数部分を用いて自己回帰係数(回帰)を計算し,チャネル暗号化後の隣接チャネル間の位相シフトを元の半分にすることで,暗号化後のデータの自己回帰係数を導出し,次いで自己回帰関係を用いて既知のデータから未知のデータ(予測)を導出することである.回顧と予測はいずれも最小二乗によって実現されるので,アルゴリズムは比較的安定である.欠点は:1、線形同相軸にしか使用できない;2.ルール・グリッドでのみ使用できます.
理論
以下の記号はSpitz論文と一致している.
L L個の同相軸とN N個のNチャネルを含む地震データは周波数領域で以下のように表すことができる.
g k ( f ) = ∑ j = 1 L a j ( f ) z j k − 1 ( f ) , k = 1 , . . . , N g_k(f)=\sum_{j=1}^La_j(f)z_j^{k-1}(f),k=1,...,N gk(f)=j=1∑Laj(f)zjk−1(f),k=1,...,N
ここでg k(f)g_k(f)gk(f)はk番目のkチャネルデータを表し、a j(f)a_j(f)aj(f)は第j j j個の同相軸に対応するサブ波スペクトル(第1パスデータ)、z j(f)=e x p(2πi f p j)z_j(f)=exp(2pi i f p_j)zj(f)=exp(2πifpj)は、第j j jの同相軸隣接チャネル間の位相シフトpj p_に対応するj pj(傾きに対応).前のブログで紹介したCanalesメソッドによれば,前−後の単一ステップ予測フィルタリングは,g k(f)=Σj=1 L P j(f)g k−j(f),k=L+1と書くことができる.N g_k(f)=\sum_{j=1}^L P_j(f)g_{k-j}(f),k=L+1,...,N gk(f)=j=1∑LPj(f)gk−j(f),k=L+1,...,N g k ∗ ( f ) = ∑ j = 1 L P j ( f ) g k + j ∗ ( f ) , k = 1 , . . . , N − L g_k^*(f)=\sum_{j=1}^L P_j(f)g_{k+j}^*(f),k=1,...,N-L gk∗(f)=j=1∑LPj(f)gk+j∗(f),k=1,...,N−LのうちP j(f)P_j(f)Pj(f)は予測フィルタの係数である.(前方-後方フィルタリングの共役関係は付録から得ることができる)
Spitz法は従来の地震路の2つの中間に新しい地震路を挿入するので,補間後の地震データはg k′(f)=Σj=1 L aj(f)z′j k−1(f),k=1,.,を満たす.2 N − 1 g_k'(f)=\sum_{j=1}^La_j(f){z'}_j^{k-1}(f),k=1,...,2N-1 gk′(f)=j=1∑Laj(f)z′jk−1(f),k=1,...,2 N−1は、g k′(f)=Σj=1 LPj′(f)g′k−j(f)、k=L+1の順方向から逆方向への単一ステップ予測フィルタ形式を再び書くことができる.2 N − 1 g'_k(f)=\sum_{j=1}^L P'_j(f){g'}_{k-j}(f),k=L+1,...,2N-1 gk′(f)=j=1∑LPj′(f)g′k−j(f),k=L+1,...,2N−1 g k ′ ( f ) = ∑ j = 1 L P ′ j ∗ ( f ) g ′ k + j ( f ) , k = 1 , . . . , 2 N − L − 1 g'_k(f)=\sum_{j=1}^L {P'}^*_j(f){g'}_{k+j}(f),k=1,...,2N-L-1 gk′(f)=j=1∑LP′j∗(f)g′k+j(f),k=1,...,2N−L−1
上記の方程式はg j’g'_についてj gj′の線形方程式群、係数はPj′P′_j Pj′は既知の値と未知の値をそれぞれ1つに配置し、整理結果:A(P′(f))[g 2′(f)g 4′(f)..g 2 N−2′(f)]T=T=B(P′(f))[g 1′(f)g 3′(f)..g 2 N−1′(f)]T A(P′(f)[g′_2(f)]T(g′_2(f)text}g′_4(f)text{}…text{}g′_{{2 N−2 N−2}(f)^T=B(P′(f)[g′_1(f)]text{xt{xt{xt}...text{}g′}g′{g′′{2 N−g′{2 N−2 N−}g'_3(f)text{}…text{}g'{2 N-1}(f)]^T A(P’(f)[g 2’(f)g 4’(f)…g 2 N−2’(f)]T=B(P′(f))[g1′(f) g3′(f) ... g2N−1′(f)]T
P’P’P’P’が既知であれば、上記式は2(2 N−L−1)2(2 N−L−1)2(2 N−L−1)個の方程式(前方後方フィルタ等式の個数)とN−1 N−1 N−1個の未知数(文献では2 N−1 2 N−1 2 N−1と表記されているが、書き間違えたのか、すべてのパスを含むのか)とを含み、最小二乗法(既知から未知へのフィルタに等価)で解くことができる.
( A + A ) − 1 A + B (A^+A)^{-1}A^+B (A+A)−1A+B
ここで、A+A^+A+は、A A A Aの複素共役転置である.
実際にはP’P’P’は未知であるが,P P Pで導くことができる.注意出力(既知の点を含む)と入力断面との位相シフト関係は、z j’(f)=z j(f/2)z'j(f)=z_j(f/2)zj’(f)=zj(f/2)(注:zj’(f)=e x p(2πi fpj/2)=zj(f/2)z'j(f)=exp(2\pi i f p_j/2)=z_j(f/2)zj’(f)=exp(2πifpj/2)=zj(f/2)であり、補間後の隣接チャネル間の位相シフトは元の半分となる.
P P Pとz z zの関係から得られる:
P j ′ ( f ) = P j ( f/2 ) P'_j(f)=P_j(f/2) Pj′(f)=Pj(f/2)
擬似周波数のない低周波数部分のみが利用されていることがわかる.
以上の導出はいくつかの詳細から与えられず、続く.
付録:自己回帰のもう一つの導出
Spitzは別の方法でCanales文の導出を行った.付導出過程は、地震データx T(f)=[x 1(f),.,x N(f)]x^T(f)=[x_1(f),...,x_N(f)]xT(f)=[x 1(f),...,xN(f)]に対して、x(f)=g(f)+n(f)x(f)=g(f)+n(f)x(f)=g(f)+n(f)x(f)=g(f)+n(f)とする.定義z変換S j T(f)=[1,z j(f),.,z j N−1(f)]S_j^T(f)=[1,z_j(f),…,z_j^{N−1}(f)]SjT(f)=[1,zj(f),…,zjN−1(f)],ここで,z j(f)=e x p(2πi f p i)z_j(f)=exp(2pi ifp_i)zj(f)=exp(2πifpi)の場合:g(f)=Σj=1 L aj(f)S j(f)g(f)=sum_{j=1}^La_j(f)S_j(f)g(f)=j=1ΣL aj(f)Sj(f)ここでa j(f)a_j(f)aj(f)は、j番目jの同方向軸に対応するサブ波のフーリエ変換である.分けて、g 1=Σj=1 L a j⋅1 g 2=Σj=1 L a j⋅z j⋯g N=Σj=1 L a j⋅z j−1begin{aligned}g_と書く1&=&\sum_{j=1}^La_j\cdot 1\\g_2&=&\sum_{j=1}^La_j\cdot z_j\\&\cdots&\\g_N&=&\sum_{j=1}^La_j\cdot z_j^{N−1}end{aligned}g 1 g 2 gN=⋯=j=1ΣL aj⋅1 j=1ΣL aj⋅zj=1ΣL aj⋅zj=1ΣL aj⋅zjN−1行列として書く形式は、(3)[g 1 g 2⋮gn]=[1⋯1 z 1 z 2⋯z L⋮z 1 N−1 z 2 N−1⋯zL−1 N−1]× [ a 1 a 2 ⋮ a n ]\left[\begin{matrix} g_1\\g_2\\\vdots\\g_n\end{matrix}\right]=\left[\begin{matrix} 1 & 1&\cdots & 1\\z_1 & z_2&\cdots & z_L\\\vdots\\z_1^{N-1} & z_2^{N-1}&\cdots & z_L^{N-1}\end{matrix}\right]\times\left[\begin{matrix} a_1\\a_2\\\vdots\\a_n\end{matrix}\right]\tag{3} ⎣⎢⎢⎢⎡g1g2⋮gn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡1z1⋮z1N−11z2z2N−1⋯⋯⋯1zLzLN−1⎦⎥⎥⎥⎤×⎣⎢⎢⎡a 1 a 2⋮an⎥⎥⎤(3)は、マトリクスをS S S、すなわちS=[S 1,⋯ ,S L]S=[S_1,cdots,S_L]S=[S 1,⋯,SL]とする.(3)を行列形式と記す:x=S a+n x=Sa+n x=Sa+nがS S S Sとして形成された行列をファンデモン行列,S j S_S_j Sj間は線形無関係であり、S S SのランクはL L L(N>=L N>=L N>=L N>=L)である.S S S Sのいずれの行も、他のL L行の線形結合として表すことができる.例えば、L＋1 L＋1 L＋1行目は、前L L行の線形の組み合わせとして表すことができる.(4)(z 1 L,⋯ ,z L)=(P L,⋯ ,P 1)と書く× [ 1 1 ⋯ 1 z 1 z 2 ⋯ z L ⋮ z 1 L − 1 z 2 L − 1 ⋯ z L L − 1 ] (z_1^L,\cdots,z_L^L)=(P_L,\cdots,P_1)\times\left[\begin{matrix} 1 & 1&\cdots & 1\\z_1 & z_2&\cdots & z_L\\\vdots\\z_1^{L-1} & z_2^{L-1}&\cdots & z_L^{L-1}\end{matrix}\right]\tag{4} (z1L,⋯,zLL)=(PL,⋯,P1)×⎣⎢⎢⎡1 z 1⋮z 1 L−1 z 2 z 2 L−1⋯⋯⋯⋯⋯⋯1 zL zLL−1⎦⎥⎤(4)上式右側を左側に移動し、z j z_について整理するj zjの方程式は,任意のz j z_に対してj zj形式は,いずれも,Z L=z L−P 1 z L−1−⋯−P L=0 Z_である.L=z^L-P_1z^{L-1}-\cdots-P_L=0 ZL=zL−P 1 zL−1−⋯−PL=0はまたz j z_j zjはそれぞれ異なり、上記方程式にはL L個の解(零点)があるので、z j z_j zjは上記方程式のKK個の零点であり、Z L=(z−z 1)(z−z 2)⋯(z−z L)Z_と書くことができるL=(z−z_1)(z−z_2)cdots(z−z_L)ZL=(z−z 1)(z−z 2)⋯(z−zL)上記2式と比較してP 1=z 1+⋯+zLP 2=−(z 1 z 2+⋯+zL−1 zL)⋯(−1)L+1 z 1 z 1 z 2⋯+zL+1 zL)P L=(−1)L+1 z 1 z 2⋯zLbegin{aligned}P_1&=&z_1+\cdots+z_L\\P_2&=&-(z_1z_2+\cdots+z_{L-1}z_L)\\&\cdots&\\P_L&=&(-1)^{L+1}z_1z_2\cdots z_Lend{aligned}P 1 P 2 PL=⋯=z 1+⋯+zL−(z 1 z 2+⋯+zL−1 zL)(−1)L+1 z 1 z 2⋯+zL上式は,我々がこれまで推さなかった一般解である.
コードの説明:
コアコード: PF = prediction_filter(x1,npf,pre1); #
[y] = interpolate_freq(x2,PF,pre2); #
自己回帰コード: for j=1:ns-np
C(j,:)=VEC(j+np:-1:j);
end
A = [C(:,2:np+1);conj(fliplr(C(:,1:np)))]; # （ ）
B = [C(:,1);conj(C(:,np+1))];#
R = A'*A;
g = A'*B;
mu = (pre/100.)*trace(R)/np;
PF = (R+mu*eye(np))\g; #
補間コード: TMPF1 = [fliplr(PF.'),-1]; #
W1 = convmtx(TMPF1,ny-np);
TMPF2 = conj(fliplr(TMPF1)); #
W2 = convmtx(TMPF2,ny-np);
WT = [W1;W2];
A = WT(:,2:2:ny); #
B = -1*WT(:,1:2:ny)*x.'; # ， ，
R= A'*A;
g = A'*B;
mu = (pre/100.)*trace(R)/(nx-1);
y1 = (R+mu*eye(nx-1))\g; #
y = zeros(1,ny);
y(1:2:ny)=x;
y(2:2:ny)=y1.';
PF = prediction_filter(x1,npf,pre1); #
[y] = interpolate_freq(x2,PF,pre2); #
for j=1:ns-np
C(j,:)=VEC(j+np:-1:j);
end
A = [C(:,2:np+1);conj(fliplr(C(:,1:np)))]; # （ ）
B = [C(:,1);conj(C(:,np+1))];#
R = A'*A;
g = A'*B;
mu = (pre/100.)*trace(R)/np;
PF = (R+mu*eye(np))\g; #
TMPF1 = [fliplr(PF.'),-1]; #
W1 = convmtx(TMPF1,ny-np);
TMPF2 = conj(fliplr(TMPF1)); #
W2 = convmtx(TMPF2,ny-np);
WT = [W1;W2];
A = WT(:,2:2:ny); #
B = -1*WT(:,1:2:ny)*x.'; # ， ，
R= A'*A;
g = A'*B;
mu = (pre/100.)*trace(R)/(nx-1);
y1 = (R+mu*eye(nx-1))\g; #
y = zeros(1,ny);
y(1:2:ny)=x;
y(2:2:ny)=y1.';