實驗三：頻率域位場處理與轉(zhuǎn)換實驗報告

〈〈重磁資料處理與解釋》實驗三頻率域位場處理和轉(zhuǎn)換實驗專業(yè)名稱：地球物理學(xué)學(xué)生姓名：學(xué)生學(xué)號：指導(dǎo)老師：王萬銀、紀(jì)新林、紀(jì)曉琳、邱世燦提交日期：2016-12-2

1基本原理1.1位場的方程2由場論知識可知，位場方程分為兩大類：有源的Possion方程U0,以及無源的Laplace方程U0。Laplace方程的第一邊值問題U|sf1通常為Dirichlet問題，第二邊值問題一|f通常稱為Nueman同題。若P點在S平面內(nèi)稱為內(nèi)部同題，反之n稱為外部問題。由唯一性定理可知，Dirichlet的內(nèi)部和外部問題的解是唯一的,而Nueman內(nèi)部問題的解不是唯一的，有一常數(shù)差，但其外部問題解是唯一的外部問題的解的唯一性的原因:0;外部問題的解的唯一性的原因:0;Un1W(p)4(1-1)W(xy,z)W(,,)x21W(p)4(1-1)W(xy,z)W(,,)x2y2d232zWGWdsn(1-2)W(x,y,)*h(x,y,z)1.2二維傅里葉變換及卷積性質(zhì)：(1-2)(1)傅里葉變換：G(u,v)Fg(x,y)i2(1)傅里葉變換：G(u,v)Fg(x,y)i2(uxvy)g(x,y)edxdy(1-3)1g(x,y)FG(u,v)i2(uxvy)dudvG(u,v)e(1-4)(2)卷積性質(zhì):iFG(u,v)*P(u,v)g(x,y)*p(x,y)(1-6)1.3頻率域位場延拓原理當(dāng)已知實測平面的異常時，換算場源以外的異常稱之為延拓，分為向上延拓和向下延拓。半空間狄利克萊問題解析解：2J(u*uv*v)(z)z、FW(x,y,z)FW(x,y,)Fh(x,y,z)FW(x,y,)e()(1-5)其中：e?”刃稱為延拓因子，z為計算面Z坐標(biāo),Z軸向下為正方向，F(xiàn)W(x,y,)為計算面頻率域位場，F(xiàn)W(x,y,z)為延拓面的頻率域位場。2輸入/輸出數(shù)據(jù)格式設(shè)計2.1輸入數(shù)據(jù)格式設(shè)計觀測面位場數(shù)據(jù)保存在filename_obser文件中，為.grd格式。計算延拓誤差時的精確場值文件保存在filename_obser2中，為.grd格式。2.2輸出數(shù)據(jù)格式設(shè)計實際計算得到的數(shù)據(jù)保存在filename_output文件中，為.grd格式。2.3參數(shù)文件數(shù)據(jù)格式設(shè)計將所要讀取的參數(shù)保存在一個文件中，該文件名變量為cmdfile,字符串變量，長度不超過80,全路徑名。在該文件中保存的參數(shù)如下：filename_obse:觀測面位場數(shù)據(jù)文件名變量filename_output:延拓后位場數(shù)據(jù)文件名變量filename_obser2延拓后準(zhǔn)確位場數(shù)據(jù)文件名變量factor_m:擴邊因子distanee延拓的高度（m/z軸向下為正方向）3總體設(shè)計此次程序采用IPO結(jié)構(gòu)設(shè)計，首先通過讀取emd文件，得到相關(guān)輸入?yún)?shù)：觀測面位場數(shù)據(jù)文件名變量、延拓后位場數(shù)據(jù)文件名變量、延拓后準(zhǔn)確位場數(shù)據(jù)文件名變量、擴邊因子、延拓的高度（m/z軸向下為正方向）；然后確定確定擴邊網(wǎng)格的大小，擴邊數(shù)據(jù)點號位置；再從觀測面位場數(shù)據(jù)文件中讀取數(shù)據(jù)。下一步，進行二維余弦擴邊，將擴完邊的數(shù)據(jù)進行快速二維傅里葉變換，轉(zhuǎn)換到頻率域；接下來計算延拓因子并且將擴完邊的數(shù)據(jù)進行快速二維傅里葉變換后在頻率域與延拓因子相乘；最后進行快速二維傅里葉反變換并且去除擴邊部分后輸出?？傮w設(shè)計見表1。圖1總體設(shè)計N-S圖4.1測試參數(shù)(1)向上延拓原始場值數(shù)據(jù)保存在’’a20_mag.中d向上延拓3.3m,延拓后理論數(shù)據(jù)保存在“a53_mag.gr中。網(wǎng)格大小為27*27(m)Xmin=-26m,Xmax=26m,Ymin=26m,Ymax=26m.有關(guān)參數(shù)保存在filename.cmd文件中，如下:A20_mag.grd,存放初始觀測面的數(shù)據(jù)文件A53_mag.grd,存放延拓后理論數(shù)據(jù)文件Output.grd,存放延拓觀測面的數(shù)據(jù)文件1,擴邊因子3.3延拓的高度(m/z軸向下為正方向)(2)向下延拓原始場值數(shù)據(jù)保存在'a53_mag.grd向下延拓3.3m，延拓后理論數(shù)據(jù)保存在“a20_mag.gr(中o網(wǎng)格大小為27*27(m)Xmin=-26m,Xmax=26m,Ymin=26m,Ymax=26m.有關(guān)參數(shù)保存在filename.cmd文件中，如下:A53_mag.grd,存放初始觀測面的數(shù)據(jù)文件A20_mag.grd,存放延拓后理論數(shù)據(jù)文件Output.grd,存放延拓觀測面的數(shù)據(jù)文件1,擴邊因子-3.3,延拓的高度(m/z軸向下為正方向)oomoomoomoooomoomooo，6，2wo—o)-!2佗1也158切4AO如說趾-825201510-r1r~20-15-10-50圖1z=-2.0km平1111111L-10-15--20-25-I-25I-20I-153X)_8aaqo0q5o05o0o5o0o5o0o50o5o0o5o0o5o0035DC02dcaaJ9-8'7"7、6-6.5如'4.3、3-2'21150I-10I-5r0X/kmI10I15~i-20~i-25圖2z=-5.3km平面理論磁異常圖252015-10-2o8850800750700650600550500450400350-5-10-15300250200150100500-20--2550100150200250300350-25-20-15-10-5510152025252015r10「-5.-10--20圖3-5.3km平面延拓磁異常圖（均方誤差152.1404）』■Lil-25--25-20-15-10-505X/km10152025圖4-2.0km平面延拓磁異常圖（均方誤差2488.719）對比圖2,圖3,可以看出，向上延拓的結(jié)果與理論磁異?；疽恢拢秸`差僅為152.1404;對比圖1和圖3可以看此向上延拓對地下3個淺部異常體的的磁異常進行了壓制；而對比圖1,圖4,向下延拓的磁異常平面圖不太光滑，與理論磁異常平面圖相比有較大誤差，為2488.719;對比圖2和圖4可以看出，向下延拓突出了地下三個淺部異常體的磁異常特征。5結(jié)論及建議由測試結(jié)果可知，向上延拓可以對淺部異常進行壓制，且誤差較??；向下延拓可突出淺部異常，且延拓結(jié)果誤差較大。經(jīng)此次試驗，對位場延拓基本步驟有了進一步了解，完成程序編寫且測試正確。附錄：程序源代碼!****************************************************************************************************!程序功能：實現(xiàn)頻率域位場延拓!cmd文件參數(shù)：!cmdfile:存放有關(guān)參數(shù)的文件名變量!filename_obser:觀測面位場數(shù)據(jù)文件!filename_output:延拓后位場數(shù)據(jù)文件!filename_obser2:延拓后準(zhǔn)確位場數(shù)據(jù)文件!factor_m:擴邊因子!distance:延拓的高度（m/z軸向下為正方向）!.grd文件參數(shù)：!N_point,N_line:點數(shù)（x方向）、線數(shù)（y方向）!x_min,x_max:x的最小值和最大值!y_min,y_max:y的最小值和最大值!Ur:初始觀測面場值!擴邊參數(shù)：!m1,m2:x方向?qū)嶋H數(shù)據(jù)起點和終點點號位置!1,m:x方向擴邊后數(shù)據(jù)起點和終點點號位置!n1,n2:y方向?qū)嶋H數(shù)據(jù)起點和終點點號位置!1,n3:y方向擴邊后數(shù)據(jù)起點和終點點號位置!延拓參數(shù):Ur:初始觀測面信號的實部!Ui:初始觀測面信號的虛部!Fr:延拓觀測面的信號的實部!Fi:延拓觀測面的信號的虛部!U:延拓后準(zhǔn)確場值!W(m,n):徑向圓頻率!Y(m,n):延拓因子!誤差參數(shù)：!error:延拓后場值與真實場值的均方誤差!****************************************************************************************************programprolongationparameter(eigval=3.701411*1e5)character*(80)cmdfile,filename_obser,filename_output,filename_obser2real,allocatable::Ur(:,:),Ui(:,:),y(:,:),Fr(:,:),Fi(:,:),U(:,:),W(:,:)integerN_point,N_lineintegerm,n,m1,m2,n1,n2integerfactor_mrealxmin,xmax,ymin,ymax,dx,dyrealdistance,error!real*8eigvalcmdfile='filename.cmd'callread_cmd(cmdfile,filename_obser,filename_output,factor_m,filename_obser2,distance)callread_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax)callcalculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(Ur(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,1:n),Fi(1:m,1:n),U(1:m,1:n),W(1:m,1:n))callinput_grd(Ur,filename_obser,m1,m2,n1,n2,m,n)callinput_grd(U,filename_obser2,m1,m2,n1,n2,m,n)callexpand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)callFFT2(Ur,Ui,m,n,2)CALLcal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)callWAVE2D(m,n,dx,dy,W)callcalculate_Y(m,n,W,y,distance)callcalculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)callFFT2(Fr,Fi,m,n,1)calloutput_grd(filename_output,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax,Ymin,Ymax)callcalculate_error(Fr,U,m1,m2,n1,n2,m,n,error)deallocate(Ur,Ui,y,Fr,Fi,u,w)endprogram****************************************************************************子程序:read_cmd功能：讀取參數(shù)文件輸入?yún)?shù)說明:

cmdfile:參數(shù)文件名輸出參數(shù)說明:!filenameobser:存放初始觀測面的數(shù)據(jù)文件:存放延拓觀測面的數(shù)據(jù)文件!filename_output!filenameobser2:存放參照觀測面的數(shù)據(jù)文件distance:延拓的高度（m/z軸向下為正方向）factorm:擴邊因子!****************************************************************************subroutineread_cmd(cmdfile,filename_obser,filename_output,factor_m,filename_obser2,distance)implicitnonecharacter*(*)cmdfile,filename_obser,filename_output,filename_obser2integerfactormrealdistance!real*8eigval!filenameobser:存放初始觀測面的數(shù)據(jù)文件:存放延拓觀測面的數(shù)據(jù)文件!filename_output!filenameobser2:存放參照觀測面的數(shù)據(jù)文件read(10,*)filename_obserread(10,*)filename_obser2read(10,*)filename_outputread(10,*)factor_mread(10,*)distanceclose(10)endsubroutinereadcmdI***************************************************************************!子程序：read_grd!功能：從原始觀測.grd文件中讀取相關(guān)參數(shù)filename_obser:輸入文件名!輸出參數(shù)說明：!N_point,N_line:點數(shù)、線數(shù)!x_min,x_max:x的最小值和最大值y_min,y_max:y的最小值和最大值subroutineread_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax)implicitnonecharacter*(*)filename_obserintegerN_point,N_linerealXmin,Xmax,Ymin,Ymaxopen(10,file=filename_obser,status='old')Read(10,*)Read(10,*)N_line,N_pointRead(10,*)Xmin,XmaxRead(10,*)Ymin,YmaxClose(10)endsubroutineread_grd!**************************************************************************!子程序：calculate_mn!功能：確定擴邊數(shù)據(jù)點號位置!輸入?yún)?shù)說明:***************************************************************************!factor_m:擴邊比例因子（＞1.0）!a,b:點數(shù)、線數(shù)!輸出參數(shù)說明：m1,m2:x方向?qū)嶋H數(shù)據(jù)起點位置和終點位置點號!m:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!n1,n2:y方向?qū)嶋H數(shù)據(jù)起點位置和終點位置點號!n:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!**************************************************************************subroutinecalculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicitnoneintegera,b,m,n,m1,m2,n1,n2integermtemp,mu,nuintegerfactor_mmtemp=aDOWHILE((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2EnddoIF(mtemp.eq.1)THENm=a*2ELSEmu=int(log(float(a))/0.693147+factor_m)m=2**muENDIFm1=1+(m-a)/2m2=m1+a-1pausemtemp=bDOWHILE((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2EnddoIF(mtemp.eq.1)THENn=b*2ELSEnu=int(log(float(b))/0.693147+factor_m)n=2**nuENDIFn1=1+(n-b)/2n2=n1+b-1write(*,*)m1,m2,n1,n2,m,npauseendsubroutinecalculate_mn!*************************************************************************!子程序：INPUT_GRD!功能：讀取grd文件中的數(shù)據(jù)!filename_obser:輸入文件名!m1,m2:x方向?qū)嶋H數(shù)據(jù)起點位置和終點位置點號!m:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!n1,n2:y方向?qū)嶋H數(shù)據(jù)起點位置和終點位置點號!n:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!輸出參數(shù)說明：!A:存放輸出數(shù)據(jù)的二維數(shù)組名!*************************************************************************SUBROUTINEINPUT_GRD(A,filename_obser,m1,m2,n1,n2,m,n)character*(*)filename_obserintegerm1,m2,n1,n2,m,nrealxmin,xmax,ymin,ymaxrealA(1:m,1:n)reali,j,kdoj=1,n,1doi=1,m,1A(i,j)=3.701411*1e10enddoenddoOpen(20,file=filename_obser,status='old')read(20,*)read(20,*)read(20,*)xmin,xmaxread(20,*)ymin,ymaxread(20,*)read(20,*)((A(i,j),i=m1,m2),j=n1,n2)Close(20)ENDSUBROUTINEINPUT_GRD!*************************************************************************!子程序：expand_cos_2D!功能：二維擴邊子程序并為信號虛部賦值!輸入?yún)?shù)說明：!m1,m2:x方向?qū)嶋H數(shù)據(jù)起點位置和終點位置點號!m:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)!n1,n2:y方向?qū)嶋H數(shù)據(jù)起點位置和終點位置點號!n:擴邊后數(shù)據(jù)終點位置點號(起點位置點號為1)Ur:初始觀測面信號的實部!Ui:初始觀測面信號的虛部!輸出參數(shù)說明：!Ur:初始觀測面信號的實部Ui:初始觀測面信號的虛部I*************************************************************************expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicitnoneintegerm,n,m1,m2,n1,n2realUr(1:m,1:n),Ui(1:m,1:n)real,allocatable::u(:),r(:)integerj,i,kallocate(u(1:m))doj=n1,n2,1doi=1,m,1u(i)=Ur(i,j)enddocallexpand_cos_1d(1,m1,m2,m,u(1))doi=1,m,1Ur(i,j)=u(i)enddoenddodeallocate(u)allocate(r(1:n))doj=1,n,1r(j)=Ur(i,j)enddocallexpand_cos_1d(1,n1,n2,n,r(1))doj=1,n,1Ur(i,j)=r(j)enddoenddodeallocate(r)doi=1,mdoj=1,nUi(i,j)=0enddoenddoendsubroutineexpand_cos_2D!**************************************************************************!子程序：expand_cos_1d!功能：一維擴邊子程序!輸入?yún)?shù)說明：n0,n3:擴邊后數(shù)據(jù)起點位置和終點位置n1,n2:實際數(shù)據(jù)起點位置和終點位置!feild(i),(i=n1,n1+1,...,n2):實際數(shù)據(jù)!輸出參數(shù)說明：!field(i),(i=n0,...,n1-1):擴邊后左邊的數(shù)據(jù)!field(i),(i=n2+1,...,n3):擴邊后右邊的數(shù)據(jù)!**************************************************************************Subroutineexpand_cos_1d(n0,n1,n2,n3,Field)RealField(n0:n3)pi=3.141592654Field(n0)=(Field(n1)+Field(n2))/2.0Field(n3)=Field(n0)i1=n0i2=n1DOi=i1,i2-1,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1))Enddoi1=n2i2=n3DOi=i1+1,i2,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1))

EnddoEndsubroutineexpand_cos_1d***********************************************************************!功能：FFT2!功能：復(fù)數(shù)組2-D快速Fourier變換EnddoEndsubroutineexpand_cos_1d!m0,m3:x方向的起點和終點!n0,n3:y方向的起點和終點!field:輸入信號(需要賦值給Freal,實部)TOC\o"1-5"\h\z!m,n:x,y方向擴邊后數(shù)據(jù)終點點號位置(起始點號為1)NF:正、反變換標(biāo)志量(1:反變換；2:正變換)輸出參數(shù)說明：!Freal：信號的實部!Fimage:信號的虛部(對于實信號而言，賦值為0)!對應(yīng)頻率分布說明：!數(shù)據(jù)Freal(m,n)和Fimage(m,n)對應(yīng)的頻率分布位置為!m方向：0,1,,m/2-1,m/2,-(m/2-1),……,-1!n方向：0,1,,n/2-1,n/2,-(n/2-1),..…,-1!************************************************************************SUBROUTINEFFT2(Freal,Fimage,m,n,nf)implicitnoneINTEGERm,n,nfREALFreal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE::Treal(:),Timage(:)integernmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0)STOPDOi=1,m,1IF(n.ne.1)THENdoj=1,n,1Treal(j)=Freal(i,j)Timage(j)=Fimage(i,j)enddocallFFT(Treal,Timage,n,nf)doj=1,n,1Freal(i,j)=Treal(j)Fimage(i,j)=Timage(j)enddoENDIFENDDODOj=1,n,1IF(m.ne.l)THENdoi=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)enddocallFFT(Treal,Timage,m,nf)doi=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)enddoENDIFENDDOdeallocate(Treal,Timage,STAT=ierr)ENDSUBROUTINEFFT2!******************************************************************!子程序：FFT!功能：復(fù)數(shù)組1-D快速Fourier變換!輸入?yún)?shù)說明：!Xreal(n):輸入數(shù)據(jù)實部!Ximage(n):輸入數(shù)據(jù)虛部!N:點數(shù)(N必須為2的窯次方)!NF:正、反變換標(biāo)志量(1:反變換；2:正變換)!輸出參數(shù)說明：!Xreal(n):變換后頻譜實部!Ximage(n):變換后頻譜虛部!對應(yīng)頻率分布說明：!數(shù)據(jù)Xreal(n)和Ximage(n)對應(yīng)的頻率分布位置為!0,1,.......,n/2-1,n/2,-(n/2-1),,-1*****************************************************************SUBROUTINEFFT(Xreal,Ximage,n,nf)implicitnoneINTEGERn,nfREALXreal(1:n),Ximage(1:n)integernu,n2,nu1,k,k1,k1n2,l,i,ibitrrealf,p,arg,c,s,treal,timagenu=int(log(float(n))/0.693147+0.001)n2=n/2nu仁nu-1f=float((-1)**nf)k=0DOl=1,nu,1DOwhile(k.lt.n)doi=1,n2,1p=ibitr(k/2**nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=Xreal(k1n2)*c+Ximage(k1n2)*stimage=Ximage(k1n2)*c-Xreal(k1n2)*sXreal(k1n2)=Xreal(k1)-trealXimage(k1n2)=Ximage(k1)-timageXreal(k1)=Xreal(k1)+trealXimage(k1)=Ximage(k1)+timagek=k+1enddok=k+n2ENDDOk=0nu1=nu1-1n2=n2/2ENDDODOk=1,n,1i=ibitr(k-1,nu)+1if(i.gt.k)thentreal=Xreal(k)timage=Ximage(k)Xreal(k)=Xreal(i)Ximage(k)=Ximage(i)Xreal(i)=trealXimage(i)=timageendifENDDOIF(nf.ne.l)THENdoi=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)enddoENDIFENDSUBROUTINEFFTFUNCTIONIBITR(J,NU)implicitnoneintegeribitr,j,nuintegerj1,itt,i,j2j1=jitt=0doi=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=itt

j1=j2enddoENDFUNCTIONIBITR!***************************************************************************!子程序：cal_dxdy!功能：計算x,y方向的點距!輸入?yún)?shù)說明：!x_min,x_max:x的最小值和最大值!y_min,y_max:y的最小值和最大值!N_point,N_line:點數(shù)（x方向）、線數(shù)（y方向）!輸出參數(shù)說明：subroutine!dx,dy:x,y方向的點距**************************************************************************cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicitnonerealxmin,xmax,ymin,ymaxintegerN_POINT,N_LINErealdx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEendsubroutinecal_dxdy*****************************************************************!子程序：WAVE2D!功能：計算2D徑向圓頻率W!輸入?yún)?shù)說明：subroutine!dx:x方向點距!dy:y方向線距!m:點數(shù)(M必須為2的窯次方)!n:線數(shù)(N必須為2的窯次方)!輸出參數(shù)說明：!W(m,n):徑向圓頻率!******************************************************************SUBROUTINEWAVE2D(m,n,dx,dy,W)implicitnoneINTEGERm,nREALdx,dyREALW(1:m,1:n),U(1:m),V(1:n)realpi,delx,delyintegermidm,midn,i,j,xx,yypi=3.1415926midm=m/2+1midn=n/2+1delx=pi*2.0/float(m)/dxdely=pi*2.0/float(n)/dydoj=1,n,1yy=jif(yy.gt.midn)yy=yy-nv(j)=dely*(yy-1)enddodoi=1,m,1xx=iif(xx.gt.midm)xx=xx-mu(i)=delx*(xx-1)enddodoj=1,n,1doi=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j))enddoenddoENDSUBROUTINEWAVE2D!*******************************************************************!子程序：calculate_Y!功能：計算延拓因子Y!輸入?yún)?shù)說明：distance:延拓的高度（m/z軸向下為正方向）m,n:x,y方向擴邊后數(shù)據(jù)終點點號位置（起始點號為1）!W:徑向圓頻率!輸出參數(shù)說明：!Y:延拓因子!*******************************************************************subroutinecalculate_Y(m,n,W,y,distance)implicitnoneintegerm,nrealpi,distancerealy(1:m,1:n),W(1:m,1:n)reali,jdoj=1,n,1Y(I,J)=1*exp(-1*w(i,j)*distance)enddoenddoendsubroutinecalculate_Y!***************************************************************************!子程序：calculate_FmulY!功能：將延拓因子與原信號做乘積!輸入?yún)?shù)說明：!m,n:x,y方向擴邊后數(shù)據(jù)終點點號位置(起始點號為1)!Y:延拓因子!Ur:初始信號的實部!Ui:初始信號的虛部!輸出參數(shù)說明：!Fr:延拓信號的實部!Fi:延拓信號的虛部***************************************************************************subroutinecalculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)implicitnoneintegerm,nrealUr(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,1:n),Fi(1:m,1:n)reali,jdoi=1,mdoj=1,nFr(i,j)=Ur(i,j)*y(i,j)Fi(i,j)=Ui(i,j)*y(i,j)enddoenddoendsubroutinecalculate_FmulY!*************************************************************************************