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

實(shí)驗(yàn)三：頻次域位場辦理與變換實(shí)驗(yàn)報告實(shí)驗(yàn)三：頻次域位場辦理與變換實(shí)驗(yàn)報告/實(shí)驗(yàn)三：頻次域位場辦理與變換實(shí)驗(yàn)報告《重磁資料辦理與解說》實(shí)驗(yàn)三頻次域位場辦理和變換實(shí)驗(yàn)專業(yè)名稱：地球物理學(xué)學(xué)生姓名：學(xué)生學(xué)號：指導(dǎo)老師：王萬銀、紀(jì)新林、紀(jì)曉琳、邱世燦提交日期：2016-12-2基本源理1.1位場的方程由場論知識可知，位場方程分為兩大類：有源的Possion方程2U0，2以及無源的Laplace方程U0。Laplace方程的第一邊值問題U|sf1平常為Dirichlet問題，第二邊值問題U|f2平常稱為Nueman問題。若P點(diǎn)在S平面內(nèi)稱為內(nèi)部問題，反之ns稱為外面問題。由獨(dú)一性定理可知，Dirichlet的內(nèi)部和外面問題的解是獨(dú)一的，而Nueman內(nèi)部問題的解不是獨(dú)一的，有一常數(shù)差，但其外面問題解是獨(dú)一的。外面問題的解的獨(dú)一性的原由：Ur0;U0。nr無源地區(qū)位場能夠表示為：W(p)1GWWGds（1-1）4nnW(x,y,z)zW(,,)32dd2x2y2z2W(x,y,)*h(x,y,z)（1-2）1.2二維傅里葉變換及卷積性質(zhì)：（）傅里葉變換：G(u,v)Fg(x,y)i2(uxvy)dxdy（1-3）1g(x,y)e1i2(uxvy)（1-4）g(x,y)FG(u,v)G(u,v)edudv（2）卷積性質(zhì)：F[g(x,y)*p(x,y)]G(u,v)*P(u,v)（1-5）F1G(u,v)*P(u,v)g(x,y)*p(x,y)（1-6）1.3頻次域位場延拓原理當(dāng)已知實(shí)測平面的異樣時，換算場源之外的異樣稱之為延拓，分為向上延拓和向下延拓。半空間狄利克萊問題解析解：FW(x,y,z)FW(x,y,)Fh(x,y,z)FW(x,y,)e2(u*uv*v)(z)（1-5）2(u*uv*v)(z)此中：e稱為延拓因子，ζ為計算面Z坐標(biāo),Z軸向下為正方向，F(xiàn)W(x,y,)為計算面頻次域位場，F(xiàn)W(x,y,z)為延拓面的頻次域位場。輸入/輸出數(shù)據(jù)格式設(shè)計2.1輸入數(shù)據(jù)格式設(shè)計觀察面位場數(shù)據(jù)保留在filename_obser文件中，為.grd格式。計算延拓偏差時的精準(zhǔn)場值文件保留在filename_obser2中，為.grd格式。2.2輸出數(shù)據(jù)格式設(shè)計實(shí)質(zhì)計算獲得的數(shù)據(jù)保留在filename_output文件中，為.grd格式。2.3參數(shù)文件數(shù)據(jù)格式設(shè)計將所要讀取的參數(shù)保留在一個文件中，該文件名變量為cmdfile，字符串變量，長度不超出80，全路徑名。在該文件中保留的參數(shù)以下：filename_obser:觀察面位場數(shù)據(jù)文件名變量filename_output：延拓后位場數(shù)據(jù)文件名變量filename_obser2:延拓后正確位場數(shù)據(jù)文件名變量factor_m：擴(kuò)邊因子distance:延拓的高度（m/z軸向下為正方向）整體設(shè)計此次程序采納IPO構(gòu)造設(shè)計，第一經(jīng)過讀取cmd文件，獲得有關(guān)輸入?yún)?shù)：觀察面位場數(shù)據(jù)文件名變量、延拓后位場數(shù)據(jù)文件名變量、延拓后正確位場數(shù)據(jù)文件名變量、擴(kuò)邊因子、延拓的高度（m/z軸向下為正方向）；此后確立確立擴(kuò)邊網(wǎng)格的大小，擴(kuò)邊數(shù)據(jù)點(diǎn)號地點(diǎn)；再從觀察面位場數(shù)據(jù)文件中讀取數(shù)據(jù)。下一步，進(jìn)行二維余弦擴(kuò)邊，將擴(kuò)完邊的數(shù)據(jù)進(jìn)行迅速二維傅里葉變換，變換到頻次域；接下來計算延拓因子而且將擴(kuò)完邊的數(shù)據(jù)進(jìn)行迅速二維傅里葉變換后在頻次域與延拓因子相乘；最后進(jìn)行迅速二維傅里葉反變換而且去除擴(kuò)邊部分后輸出。整體設(shè)計見表1。輸入?yún)?shù)：filename_obser:觀察面位場數(shù)據(jù)文件名變量、filename_output：延拓后位場數(shù)據(jù)文件名變量、filename_obser2:延拓后正確位場數(shù)據(jù)文件名變量、factor_m：擴(kuò)邊因子、distance:延拓的高度（m/z軸向下為正方向）確立擴(kuò)邊網(wǎng)格的大小m*n(m,n均為2的冪次方)從觀察面位場數(shù)據(jù)文件中讀取數(shù)據(jù)二維擴(kuò)邊子程序并為信號虛部賦值擴(kuò)邊后的數(shù)據(jù)進(jìn)行二維迅速傅里葉正變換計算延拓因子將擴(kuò)完邊并進(jìn)行迅速二維傅里葉正變換后的復(fù)數(shù)數(shù)組在頻次域與延拓因子相乘對上步結(jié)果進(jìn)行迅速二維傅里葉反變換去除擴(kuò)邊部分后輸出結(jié)果圖1整體設(shè)計N-S圖測試結(jié)果4.1測試參數(shù)（1）向上延拓原始場值數(shù)據(jù)保留在’’a20_mag.中g(shù)rd，”向上延拓3.3m，延拓后理論數(shù)據(jù)保存在“a53_mag.grd中”。網(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，擴(kuò)邊因子3.3延拓的高度（m/z軸向下為正方向）（2）向下延拓原始場值數(shù)據(jù)保留在’’a53_mag.中g(shù)rd，”向下延拓3.3m，延拓后理論數(shù)據(jù)保留在“a20_mag.grd中”。網(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，擴(kuò)邊因子-3.3，延拓的高度（m/z軸向下為正方向）4.2測試結(jié)果25201510512m/0Y-53

2000180016001400120010008006004002000-10-200-400-15-600-800-20-25-25-20-15-10-50510152025X/km圖1z=-2.0km平面理論磁異樣圖25201590085080075010700650600155055004502400m350k0300Y250200-51503100500-10-50-100-150-200-15-250-300-350-20-25-25-20-15-10-50510152025X/km圖2z=-5.3km平面理論磁異樣圖2520151051m2kY-53-10-15-20-25

850800750700650600550500450400350300250200150100500-50-100-150-200-250-300-350-25-20-15-10-50510152025圖3-5.3km平面延拓磁異樣圖（均方偏差152.1404）252015200018001016001400512m/0Y-53

120010008006004002000-10-200-400-15-600-800-20-25-25-20-15-10-50510152025X/km圖4-2.0km平面延拓磁異樣圖（均方偏差2488.719）比較圖2，圖3，能夠看出，向上延拓的結(jié)果與理論磁異樣基本一致，均方偏差僅為

152.1404；比較圖

1和圖

3能夠看出，向上延拓對地下

3個淺部異樣體的的磁異樣進(jìn)行了壓制；而比較圖1，圖滑，與理論磁異樣平面圖比較有較大偏差，為

4，向下延拓的磁異樣平面圖不太光2488.719；比較圖2和圖4能夠看出，向下延拓突出了地下三個淺部異樣體的磁異樣特點(diǎn)。結(jié)論及建議由測試結(jié)果可知，向上延拓能夠?qū)\部異樣進(jìn)行壓制，且偏差較??；向下延拓可突出淺部異樣，且延拓結(jié)果偏差較大。經(jīng)此次試驗(yàn)，對位場延拓基本步驟有了進(jìn)一步認(rèn)識，達(dá)成程序編寫且測試正確。附錄：程序源代碼!****************************************************************************************************!程序功能：實(shí)現(xiàn)頻次域位場延拓cmd文件參數(shù):cmdfile:寄存有關(guān)參數(shù)的文件名變量filename_obser:觀察面位場數(shù)據(jù)文件!filename_output：延拓后位場數(shù)據(jù)文件filename_obser2:延拓后正確位場數(shù)據(jù)文件!factor_m：擴(kuò)邊因子distance:延拓的高度（m/z軸向下為正方向）.grd文件參數(shù)：!N_point,N_line：點(diǎn)數(shù)(x方向)、線數(shù)(y方向)!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值Ur:初始觀察面場值擴(kuò)邊參數(shù)：!m1,m2：x方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號地點(diǎn)!1,m：x方向擴(kuò)邊后數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號地點(diǎn)!n1,n2：y方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號地點(diǎn)!1,n3：y方向擴(kuò)邊后數(shù)據(jù)起點(diǎn)和終點(diǎn)點(diǎn)號地點(diǎn)延拓參數(shù)：Ur:初始觀察面信號的實(shí)部!Ui：初始觀察面信號的虛部!Fr：延拓觀察面的信號的實(shí)部!Fi：延拓觀察面的信號的虛部!U：延拓后正確場值!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ù)說明：!filename_obser：寄存初始觀察面的數(shù)據(jù)文件!filename_output：寄存延拓觀察面的數(shù)據(jù)文件!filename_obser2：寄存參照觀察面的數(shù)據(jù)文件distance:延拓的高度（m/z軸向下為正方向）!factor_m：擴(kuò)邊因子!****************************************************************************subroutineread_cmd(cmdfile,filename_obser,filename_output,factor_m,filename_obser2,distance)implicitnonecharacter*(*)cmdfile,filename_obser,filename_output,filename_obser2integerfactor_mrealdistance!real*8eigvalcharacter*80stropen(10,file=cmdfile,status='old')read(10,*)filename_obserread(10,*)filename_obser2read(10,*)filename_outputread(10,*)factor_mread(10,*)distanceclose(10)endsubroutineread_cmd!***************************************************************************!子程序：read_grd功能：從原始觀察.grd文件中讀取有關(guān)參數(shù)輸入?yún)?shù)說明：!filename_obser：輸入文件名輸出參數(shù)說明：!N_point,N_line：點(diǎn)數(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功能：確立擴(kuò)邊數(shù)據(jù)點(diǎn)號地點(diǎn)輸入?yún)?shù)說明：!factor_m:擴(kuò)邊比率因子（>1.0）!a,b：點(diǎn)數(shù)、線數(shù)輸出參數(shù)說明：!

m1,m2:x

方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)點(diǎn)號!

m:擴(kuò)邊后數(shù)據(jù)終點(diǎn)地點(diǎn)點(diǎn)號

(起點(diǎn)地點(diǎn)點(diǎn)號為

1)!

n1,n2:y

方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)點(diǎn)號!n:擴(kuò)邊后數(shù)據(jù)終點(diǎn)地點(diǎn)點(diǎn)號(起點(diǎn)地點(diǎn)點(diǎn)號為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-1write(*,*)m,apausemtemp=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ù)輸入?yún)?shù)說明：!filename_obser：輸入文件名!m1,m2:x方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)點(diǎn)號!m:擴(kuò)邊后數(shù)據(jù)終點(diǎn)地點(diǎn)點(diǎn)號(起點(diǎn)地點(diǎn)點(diǎn)號為1)!n1,n2:y方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)點(diǎn)號!n:擴(kuò)邊后數(shù)據(jù)終點(diǎn)地點(diǎn)點(diǎn)號(起點(diǎn)地點(diǎn)點(diǎn)號為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功能：二維擴(kuò)邊子程序并為信號虛部賦值輸入?yún)?shù)說明：!m1,m2:x方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)點(diǎn)號!m:擴(kuò)邊后數(shù)據(jù)終點(diǎn)地點(diǎn)點(diǎn)號(起點(diǎn)地點(diǎn)點(diǎn)號為1)!n1,n2:y方向?qū)嵸|(zhì)數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)點(diǎn)號!n:擴(kuò)邊后數(shù)據(jù)終點(diǎn)地點(diǎn)點(diǎn)號(起點(diǎn)地點(diǎn)點(diǎn)號為1)Ur:初始觀察面信號的實(shí)部!Ui：初始觀察面信號的虛部輸出參數(shù)說明：Ur:初始觀察面信號的實(shí)部!Ui：初始觀察面信號的虛部!*************************************************************************subroutineexpand_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))doi=1,m,1doj=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功能：一維擴(kuò)邊子程序輸入?yún)?shù)說明：!n0,n3：擴(kuò)邊后數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)!n1,n2：實(shí)質(zhì)數(shù)據(jù)起點(diǎn)地點(diǎn)和終點(diǎn)地點(diǎn)!feild(i),(i=n1,n1+1,...,n2)：實(shí)質(zhì)數(shù)據(jù)輸出參數(shù)說明：!field(i),(i=n0,...,n1-1)：擴(kuò)邊后左側(cè)的數(shù)據(jù)!field(i),(i=n2+1,...,n3)：擴(kuò)邊后右側(cè)的數(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變換輸入?yún)?shù)說明：!m0,m3：x方向的起點(diǎn)和終點(diǎn)!n0,n3：y方向的起點(diǎn)和終點(diǎn)!field：輸入信號(需要賦值給Freal，實(shí)部)!m,n：x,y方向擴(kuò)邊后數(shù)據(jù)終點(diǎn)點(diǎn)號地點(diǎn)（初步點(diǎn)號為1）!NF:正、反變換標(biāo)記量(1:反變換;2:正變換)輸出參數(shù)說明：!Freal：信號的實(shí)部!Fimage：信號的虛部(關(guān)于實(shí)信號而言，賦值為0)對應(yīng)頻次散布說明：!

數(shù)據(jù)

Freal(m,n)

和

Fimage(m,n)

對應(yīng)的頻次散布地點(diǎn)為

:!

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.1)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ù)實(shí)部!

Ximage(n):

輸入數(shù)據(jù)虛部!

N:點(diǎn)數(shù)（N

必然為

2的冪次方）!

NF:

正、反變換標(biāo)記量

(1:反變換

;2:正變換

)輸出參數(shù)說明：!Xreal(n):變換后頻譜實(shí)部!Ximage(n):變換后頻譜虛部對應(yīng)頻次散布說明：!

數(shù)據(jù)

Xreal(n)

和

Ximage(n)

對應(yīng)的頻次散布地點(diǎn)為

:!

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/2nu1=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.1)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=ittj1=j2enddoENDFUNCTIONIBITR!***************************************************************************!子程序：cal_dxdy功能：計算x,y方向的點(diǎn)距輸入?yún)?shù)說明：!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!N_point,N_line：點(diǎn)數(shù)(x方向)、線數(shù)(y方向)輸出參數(shù)說明：!dx,dy：x,y方向的點(diǎn)距!***************************************************************************subroutinecal_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ù)說明：!dx:x方向點(diǎn)距!dy:y方向線距!m:點(diǎn)數(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方向擴(kuò)邊后數(shù)據(jù)終點(diǎn)點(diǎn)號地點(diǎn)（初步點(diǎn)號為1）!W：徑向圓頻次輸出參數(shù)說明：!Y：延拓因子!*******************************************************************subroutinecalculate_Y(m,n,W,y,distance)implicitnoneintegerm,nrealpi,distancerealy(1:m,1:n),W(1:m,1:n)reali,jpi=3.1415926doi=1,m,1doj=1,n,1Y(I,J)=1*exp(-1*w(i,j)*distance)enddoenddoendsubroutinecalculate_Y!***************************************************************************!子程序：calculate_FmulY功能：將延拓因子與原信號做乘積輸入?yún)?shù)說明：!m,n：x,y方向擴(kuò)邊后數(shù)據(jù)終點(diǎn)點(diǎn)號地點(diǎn)（初步點(diǎn)號為1）!Y：延拓因子!Ur：初始信號的實(shí)部!Ui：初始信號的虛部輸出參數(shù)說明：!Fr：延拓信號的實(shí)部!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!******************************************************************************