频率域位场处理与转换(12.25-)
实验二延拓——精选推荐
实验二频率域位场处理和转换实验******专业：勘查技术与工程学号：**********指导老师：王万银，纪新林目录一、基本原理 (1)1.1空间域延拓原理 (1)1.2波数域延拓原理 (1)二、输入/输出数据格式设计 (2)1、输入数据 (2)2、输出数据 (2)三、总体设计 (2)1、程序流程图 (2)2、参数说明 (3)四、测试结果 (3)五、结论及建议 (7)附录 (8)一、基本原理1.1空间域延拓原理延拓分为：向上延拓和向下延拓。
向上延拓是只将实测平面上的数据换到平面之上的平面上的计算。
对于二度体（令z 向下为正）：U(x,z)=-zπ U (ξ,0)(ξ−x)2+z 2d ξ +∞−∞（z<0） （1）其中U （ξ，0）为剖面上各点的实测值。
Z 为延拓的高度，即转换后的平面和观测平面的距离，由于z 坐标向下为正，因此z<0。
空间域的延拓实际是积分的计算。
向下延拓是由实测磁场向磁源方向延拓。
其计算公式和（1）相同。
1.2波数域延拓原理设场源位于z=H 平面以下（H>0），则磁场为在z=H 平面以上对x 、y 、z 连续的函数，若z=0观测平面的磁场T （x ，y ，0）为已知，则由外部的狄利克莱问题：T （x ，y ，z ）=−z2π T （ξ，η，0）[(x −ξ)2+(y −η)2+z 2]3/2+∞−∞+∞−∞d ξd η （2）对其进行变量x ，y 进行傅里叶变换成S T u ，v ，z = T （x ，y ，z ）e −2πi(ux +vy )dxdy ∞−∞∞−∞ （3） 因此：S T u ，v ，z = S T u ，v ，0 e 2π(u2+v 2)1/2z（4）在波数域中，延拓变成了对观测数据的傅里叶变换乘以延拓因子。
二、输入/输出数据格式设计1、输入数据（1）控制变量名：’cmd.txt’（2）观测面高度之差：height=3.3m或height=-3.3m（3）低高度网格化数据：从文件“A20_mag.grd”输入（4）高高度网格化数据：从文件“A53_mag.grd”输入（5）特征值：eigval=1.701411∗1038（6）输出文件名：out.grd2、输出数据通过修改参数，运行两次程序，向上延拓结果及向下延拓结果均输入“output.grd”先做向上延拓，将结果利用surfer画图后，再做向下延拓三、总体设计1、程序流程图图3.1 程序流程框图2、参数说明mpoint,line 点数，线数height 观测面高度差xmin,xmax,ymax,ymin x方向最小值、最大值，y方向最大值、最小值eigval 特征值m0,m1,m2,m3 x方向扩边点位n0,n1,n2,n3 y方向扩边点位Field 场值Field_real,Field_image 傅里叶变换中场值的实部，虚部Fconti_real Fconti_image 延拓因子实部、虚部四、测试结果1、原场值图：图4.1 低高度场值图图4.2 高高度场值2图4.3 向上延拓结果图4.4 向下延拓结果分析：将底高度向上延拓3.3m 后结果与高高度原场值图形近似，但是场值变小了，且异常部位相对不明显了，将高高度向下延拓3.3m 后，异常部位没变，但是场值绝对值变大了近5倍，且在边缘部位有几个异常点。
频率域变换
数字图像处理本章包含的主要内容傅立叶变换 卷积和卷积定理 频率域低通滤波 频率域高通滤波2问题1：傅立叶变换•空间域/灰度•频率域/幅值与频率4• 傅立叶变换的预备知识¾ 点源和狄拉克函数一幅图像可以看成由无穷多像素组成，每个像素可以看成 一个点源， 点源可以用狄拉克函数δ表示：⎧∞ δ ( x, y ) = ⎨ ⎩0∞x = 0, y = 0 其他ε满足−∞∫ ∫ δ ( x, y ) dxdy = ∫ ∫ ε δ ( x, y ) dxdy = 1−ε为任意小的正数5• 狄拉克函数δ具性有的性质9 δ函数为偶函数δ ( − x, − y ) = δ ( x, y )∞ ∞9位移性 或f ( x, y ) =−∞ −∞∫∫f (α , β )δ ( x − α , y − β ) d α d βf ( x, y ) = f ( x, y ) ∗ δ ( x, y )9 9可分性 筛选性δ ( x, y ) = δ ( x)δ ( y )f (α , β ) =∞ ∞ −∞ −∞ ∞ ∞∫∫f ( x, y )δ ( x − α , y − β )dxdy当α=β=0时f (0, 0) =−∞ −∞∫∫f ( x, y )δ ( x, y )dxdy6¾ 二维线性不变系统满足线性和齐次性条件的系统称为二维线性系统 9 9 线性T[ f1(x, y) + f2 (x, y)] = T[ f1(x, y)] +T[ f2 (x, y)]齐次性T [ af ( x, y )] = aT [ f ( x, y )]9二维线性系统T[a1 f1(x, y) +a2 f2(x, y)] = aT 1 [ f1(x, y)] + a2T[ f2 (x, y)]7 点扩散函数 当输入为单位脉冲δ(x,y),系统的输出为脉冲响应， 用h(x,y)表示，在图像处理中用作对点源的响应， 称为点扩散函数 位移不变系统 当输入的脉冲响应延迟了α，β单位，即δ(x- α,yβ)，如果输出为h(x- α,y- β)，称此系统为位移不 变系统 对于位移不变系统，系统的输出仅和输入函数的性 态有关，和作用的起点无关8• 二维线性位移不变系统输入，输出关系图f (x, y) g(x,y)=f (x, y)* h ( x, y ) T∞ ∞g ( x , y ) = T [ f ( x , y )] = T−∞ −∞ ∞ ∞∫∫f (α , β )δ ( x − α , y − β ) d α d β线性 位移不变= =−∞ −∞ ∞ ∞∫∫ ∫∫f (α , β )T [δ ( x − α , y − β )]d α d β f (α , β ) h ( x − α , y − β ) d α d β−∞ −∞线性位移不变系统的输出等于系统的输入和系统脉冲响 应（点扩散函数）的卷积9¾傅立叶积分的定义 • 调谐信号：ejω t= cos( ω t ) + j sin( ω t )其中j2=-1 变换因子• 傅立叶积分：H ( f ) =∫∞ −∞h ( t )e− j 2 π ftdt其中t代表时间，f代表频率10振幅图形傅立叶变换为：),00v v u u −−)]//00N vy M ux +把图像进行傅立叶变换后，往往要把中心移到∑∑−−=11),(1),(N N y x f y x f 8、平均值问题2：卷积和卷积定理¾g（x－α）32¾相关的几何解释g(x1+α)g(2x1+α)g(3x1+α)g(+α)g(4x1+α) g(-x1+α)g(5x1+α)37389相关描述的是两个函数图像的形似程度，当完全相同时，相关函数就会出现一个相关峰值。
信号时域频域及其转换
信号分析方法概述：通用的基础理论是信号分析的两种方法：1 是将信号描述成时间的函数 2 是将信号描述成频率的函数。
也有用时域和频率联合起来表示信号的方法。
时域、频域两种分析方法提供了不同的角度，它们提供的信息都是一样，只是在不同的时候分析起来哪个方便就用哪个。
思考：原则上时域中只有一个信号波（时域的频率实际上是开关器件转动速度或时钟循环次数，时域中只有周期的概念），而对应频域（纯数学概念）则有多个频率分量。
人们很容易认识到自己生活在时域与空间域之中（加起来构成了三维空间），所以比较好理解时域的波形（其参数有：符号周期、时钟频率、幅值、相位）、空间域的多径信号也比较好理解。
但数学告诉我们，自己生活在N维空间之中，频域就是其中一维。
时域的信号在频域中会被对应到多个频率中，频域的每个信号有自己的频率、幅值、相位、周期（它们取值不同，可以表示不同的符号，所以频域中每个信号的频率范围就构成了一个传输信道。
时域中波形变换速度越快（上升时间越短），对应频域的频率点越丰富。
所以：OFDM中，IFFT把频域转时域的原因是：IFFT的输入是多个频率抽样点（即各子信道的符号），而IFFT之后只有一个波形，其中即OFDM符号，只有一个周期。
时域时域是真实世界，是惟一实际存在的域。
因为我们的经历都是在时域中发展和验证的，已经习惯于事件按时间的先后顺序地发生。
而评估数字产品的性能时，通常在时域中进行分析，因为产品的性能最终就是在时域中测量的。
时钟波形的两个重要参数是时钟周期和上升时间。
时钟周期就是时钟循环重复一次的时间间隔，通产用ns度量。
时钟频率Fclock，即1秒钟内时钟循环的次数，是时钟周期Tclock的倒数。
Fclock=1/Tclock上升时间与信号从低电平跳变到高电平所经历的时间有关，通常有两种定义。
一种是10-90上升时间，指信号从终值的10%跳变到90%所经历的时间。
这通常是一种默认的表达方式，可以从波形的时域图上直接读出。
数字图像处理数字图像处理第二章(第二讲)空域变换、频率域变换
2.1 引言 2.2 空域变换 2.3 频率域变换 2.4 离散余弦变换 2.5 KL变换 2.6 其他正交变换
第二章 常用的数学变换
2.2 空域变换——2.2.2. 遥感影像几何校正 2.2.2.2 几何校正类型
商用遥感数据(如：SPOT-Image， Digital Globe， Space Image)都已消 除了大多数系统误差 。两种常用的几何校正方法：
第二章 常用的数学变换
2.2 空域变换——2.2.2. 几何变换
第二章 常用的数学变换
2.2 空域变换——2.2.2. 几何变换
第二章 常用的数学变换
2.2 空域变换——2.2.2. 几何变换 灰度插值——高阶插值
如果简化计算，仅取原点周围有限范围函数：
第二章 常用的数学变换
2.2 空域变换——2.2.2. 几何变换
2.2.2.2 几何校正类型 ➢ 从影像到影像的配准
从影像到影像的配准 是平移和旋转过程的结合，通过 两幅影像中的同名点进行匹配，使同名地物出现在配准后 的 影像的相同位置。若不需要使每个像元都具有特定的地 图投影坐标(x， y)，就可以使用这种纠正方法。例如，我们 可能希望用光标查看两幅不同时相影像是否发生变化。
国家级精品资源共享课
2.2.2.2 几何校正类型
➢ 影像校正/配准的混合方法
影像校正和配准所用的基本原理是相同的。所不同的是：从影像到 地图的校正中，参考的是有标准地图投影的地图；而从影像到影像的配 准中，参考的是另一景影像。 如果采用已校正过的影像(而不是传统地 图)为参考，那么得到的所有配准影像都会带有原参考影像中包含的几 何误差。 因此，高精度地球科学遥感研究中，应采用从影像到地图的 校正。 然而，对两个或多个遥感数据进行精确的变化检测时，选择从 影像到地图的校正和从影像到影像的配准相结合的混合校正法就显得十 分有用。
频域和时域的转换公式
频域和时域的转换公式
（）
近年来，互联网及相关科技的迅猛发展推动了社会全面的发展，频域和时域的
转换也被广泛的应用到各个领域当中去。
频域和时域的转换即频率域转换到时域转换或者时域转换到频率域转换，它是将信号中影响信号特性的因素从一个域中提取到另一个域进行分析的过程。
频率域转换到时域转换是将频率域（作用于不同频率的、随时间变化的信号）
中的特征映射到时域的信号的一个过程，其常用的转换公式为傅里叶变换
（Fourier Transform）、拉普拉斯变换（Laplace Transform）和索尔兹变换（Z Transform）。
傅里叶变换是一种几乎可以算出任何频率域函数对应的时域函数形
式的有限次数级数，它建立了一种单独定量频率信号特性并表示为持续信号测量的新方式。
拉普拉斯变换和索尔兹变换也属于线性时不变系统变换，它们可以在采样频率和持续频率之间进行转换，从而实现连续调制信号以模拟技术进行传输的功能。
时域转换到频率域转换是将时域（系统的输入/输出特性）中的特征映射到频
域的信号，通常采用傅里叶逆变换（Inverse Fourier Transform）和拉普拉斯逆
变换（Inverse Laplace Transform）来实现。
傅里叶逆变换可以将频率域的调制
信号的频率转换为传输信号的时域表示，而拉普拉斯逆变换可以将持续频率的输入信号转换为离散信号的频域形式，从而实现从持续信号到采样信号或从采样信号到持续信号的变换。
总之，频域和时域的转换技术在互联网领域有着广泛的应用，它可以帮助我们
更好的分析和理解信号的特性，从而提升信号传输的品质及改善相关科技的发展。
时域和频域的转换公式
时域和频域的转换公式时域和频域是信号处理中常用的两个概念。
时域描述了信号在时间轴上的变化情况，而频域描述了信号在频率轴上的变化情况。
两者之间存在着转换关系，通过转换公式可以将时域信号转换为频域信号，或者将频域信号转换为时域信号。
一、时域信号与频域信号的定义1.时域信号：时域信号是指信号在时间轴上的变化情况。
时域信号可以表示为x(t)，其中t表示时间，x(t)表示在时间t时刻信号的幅值。
2.频域信号：频域信号是指信号在频率轴上的变化情况。
频域信号可以表示为X(f)，其中f表示频率，X(f)表示在频率f上的信号功率。
二、傅里叶变换与傅里叶逆变换傅里叶变换是将时域信号转换为频域信号的数学工具，傅里叶逆变换则是将频域信号转换为时域信号的数学工具。
1.傅里叶变换：傅里叶变换可以将一个时域信号x(t)转换为频域信号X(f)，其公式为：X(f) = ∫[x(t) * e^(-j2πft)] dt其中，∫表示积分符号，e为自然对数的底数，f为频率，j为虚数单位。
2.傅里叶逆变换：傅里叶逆变换可以将一个频域信号X(f)转换为时域信号x(t)，其公式为：x(t) = ∫[X(f) * e^(j2πft)] df其中，∫表示积分符号，e为自然对数的底数，f为频率，j为虚数单位。
三、快速傅里叶变换快速傅里叶变换（FFT）是一种计算傅里叶变换和逆变换的高效算法，它可以大幅度减少计算量。
FFT算法将信号分解为多个频率块，通过对这些频率块进行傅里叶变换，最后将它们合并成一个完整的频域信号。
FFT算法的关键思想是将一个长度为N的离散时域信号转换为长度为N的离散频域信号。
FFT有两种形式：正向FFT和反向FFT。
正向FFT将时域信号转换为频域信号，而反向FFT则将频域信号转换为时域信号。
显示如下为正向FFT公式：X(k) = Σ[x(n) * e^(-j2πkn/N)]，其中k为频率索引，N为时域信号的长度，n为时间索引。
反向FFT公式：x(n) = (1/N) * Σ[X(k) * e^(j2πkn/N)]，其中k为频率索引，N为时域信号的长度，n为时间索引。
时域滤波器和频域滤波器的变换
时域滤波器和频域滤波器的变换卷积定理函数空间域的卷积的傅⾥叶变换是函数傅⾥叶变换的乘积。
对应地，频率域的卷积与空间域的乘积存在对应关系。
由卷积定理可知所有频域的滤波理论上都可以转化为空域的卷积操作。
给定频率域滤波器，可对其进⾏傅⾥叶逆变换得到对应的空域滤波器；滤波在频域更为直观，但空域适合使⽤更⼩的滤波模板以提⾼滤波速度。
因为相同尺⼨下，频域滤波器效率⾼于空域滤波器，故空域滤波需要⼀个更⼩尺⼨的模板近似得到需要的滤波结果。
空域卷积将模板在图像中逐像素移动，将卷积核的每个元素分别和图像矩阵对应位置元素相乘并将结果累加，累加和作为模板中⼼对应像素点的卷积结果。
通俗的讲，卷积就是对整幅图像进⾏加权平均的过程，每⼀个像素点的值，都由其本⾝和邻域内的其他像素值经过加权平均后得到。
在像素的处理上，是先将结果暂存在于⼀个副本，最后统⼀拷贝，故不会出现处理顺序不同⽽结果不同的情况。
⼆维连续卷积的数学定义：离散形式：频域滤波频率域是由傅⾥叶变换和频率变量 (u,v)定义的空间，频域滤波处理过程：先对图像进⾏傅⾥叶变换，转换⾄频率域，在频域使⽤滤波函数进⾏滤波，最后将结果反变换⾄空间域。
即：⾼斯函数公式：形状：空域⾼斯平滑滤波⾼斯模板的⽣成因为图像是离散存储的，故我们需要⼀个⾼斯函数的离散近似。
具体地，对⾼斯函数进⾏离散化，以离散点上的⾼斯函数值作为权值，组成⼀定尺⼨的模板，⽤此模板对图像进⾏卷积。
由于⾼斯分布在任意点处都⾮零，故理论上需要⼀个⽆穷⼤的模板，但根据" 准则"，即数据分布在的概率是0.9974，距离函数中⼼超过数据所占权重可以忽略,因此只需要计算的矩阵就可以保证对⾼斯函数的近似了。
假设⼆维模板⼤⼩,则模板上元素处的值为：前⾯的系数在实际应⽤中常被忽略，因为是离散取样，不能使取样和为1，最后还要做归⼀化操作。
程序：function filt=mygaussian(varargin)%参数初始化，使⽤varargin处理可变参数情况siz=varargin{1};%模板尺⼨if(numel(siz)==1)siz=[siz,siz];endstd=varargin{2};%⽅差centa = (siz(1)+1)/2;%此处不要取整centb = (siz(1)+1)/2;filt = zeros(siz(1),siz(2));summ=0;for i=1:siz(1)for j=1:siz(2)radius = ((i-centa)^2+(j-centb)^2);filt(i,j) = exp(-(radius/(2*std^2)));summ=summ+filt(i,j);endendfilt=filt/summ;%归⼀化测试：执⾏mygaussian(4,1)得：0.0181 0.0492 0.0492 0.01810.0492 0.1336 0.1336 0.04920.0492 0.1336 0.1336 0.04920.0181 0.0492 0.0492 0.0181执⾏fspecial('gaussian',4,1)得：0.0181 0.0492 0.0492 0.01810.0492 0.1336 0.1336 0.04920.0492 0.1336 0.1336 0.04920.0181 0.0492 0.0492 0.0181可以看出与Matlab结果相同。
频率域磁异常处理与转换
磁法勘探上机实验报告*名：***学号： ********** 指导教师：***日期：2020.4.20一、实验目的1、加深对磁性体磁异常在频率域处理转换原理与作用的认识；2、用Matlab语言编程实现频率域磁异常处理与转换，如向上延拓、导数计算、∆T磁异常化极等处理，培养程序开发与数据处理的动手能力。
二、实验内容1、利用两个大小与埋深不同的球体产生的平面磁异常数据（如ΔT、Za），选择频率域向上延拓算子、导数计算算子、化磁极算子，通过频率域实现磁异常的处理与转换；2、磁异常数据准备：自行准备，利用正演实验计算球体的磁异常数据（ΔT、Za）；3、频率域处理转换算子：向上延拓算子、垂向二阶导数算子、化磁极算子。
1）球体正演参数自行设定，也可参考以下参数设置：（1）正演2个球体的叠加异常（如ΔT、Za），用于向上延拓、垂向二阶导数计算、化磁极的数据源；（2）平面磁异常计算范围：x：-200m至200m，y：-200m至200m，地磁场总强度T=5*10−9T，磁倾角I=45°，磁偏角D=0°；（3）假设球体1的参数：半径r1=2m,球心坐标（50m,0m,5m）,磁化强度M1= 0.1 A/m；（4）假设球体2的参数：半径r1=15m,球心坐标（-50m,0m,30m）,磁化强度M2= 0.2 A/m；2）如上数据可以存储为grd文件，处理转换时直接调用即可。
三、实验要求球体某分量磁异常上延计算、导数计算、化磁极可仟选其一，具体要求如下：1、上延计算：对ΔT或Za向上延拓不同高度，画出对应结果图；2、导数计算：利用ΔT或Za异常计算垂向二阶导数，画出对应结果图；3、化磁极：利用ΔT数据进行化极处理，画出对应结果图；4、观察转换前后的异常特征，分析这些处理转换的作用。
四、实验原理在频率域实现磁异常处理转换计算简便，速度快。
通过频率域处理转换的基本过程为：1）先将对磁异常进行傅里叶变换求磁异常频谱；2）将磁异常频谱与频域处理转换的算子相乘，得到转换后的频谱；3）对转换后的频谱进行反傅里叶变换得到磁异常。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
! call search_unit(10,nunit_in)
OPEN(11,file=cmd_filename,status='old')
（2）观测面高度之差：height=3.3m
（3）低高度网格化数据：A20_mag.grd
（4）特征值：
（6）输出文件名：expound_a20_2D.grd
2.
（1）输入文件：input_filename
(2) 输出文件：out_filename
(3) 点数：mpoint
(4)线数：nline
(5)异常值：filed
SUBROUTINE Expound_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3,line)
integer m0,m1,m2,m3,n0,n1,n2,n3
real Field(m0:m3,n0:n3)
real::pi=3.141592654
do i=n1,n2
Field(m0,i)=(Field(m1,i)+Field(m2,i))/2
重磁实验报告
实验名称：频率域位场处理与转换
*******
专业名称：勘查技术与工程
学生学号：************
指导老师：巨鹏 王万银
一
（1）实验内容：
对一个网格化数据进行向上延拓，得到延拓后的结果，并计算延拓后的均方误差。
（2）要求：
进行软件设计、代码编写、软件测试、编写实验报告。
二．基本原理
分别介绍空间域和频率域位场延拓的计算公式。
七．源程序及代码
PROGRAM wavenumber_continuation
INTEGER mpoint,line,m0,m1,m2,m3,n0,n1,n2,n3,m,n
REAL height,eigval,Xmin,Xmax,Ymin,Ymax,dx,dy
CHARACTER*80 input_filename,output_filename,cmd_filename,Expound_a20_2D_filename
character*(*) input_filename
integer m0,m1,m2,m3,n0,n1,n2,n3
real eigval
real Field(m0:m3,n0:n3)
Field=eigval
! call search_unit(10,nunit_in)
open(13,file=input_filename,form='formatted')
CALL Input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field,eigval)
CALL Expound_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3,line)
CALL Putout_GRD_Expound_A20_2D(Field,Expound_a20_2D_filename,m0,m3,n0,n3,m,n,eigval,Xmin,Xmax,Ymin,Ymax)
end do
do j=n2+1,n3-1
Field(i,j)=Field(i,n2)+cos(pi/2.0*(n3-j)/(n3-n2))*(Field(i,n3)-Field(i,n2))
end do
end do
END SUBROUTINE
!*******************输出扩边Expound_a20_2D.grd***********************
end do
end do
do i=m0,m3
Field(i,n0)=(Field(i,n1)+Field(i,n2))/2
Field(i,n3)=Field(i,n0)
end do
do i=m0,m3
do j=n0+1,n1-1
Field(i,j)=Field(i,n0)+cos(pi/2.0*(n1-j)/(n1-n0))*(Field(i,n1)-Field(i,n0))
REAL,ALLOCATABLE::Field(:,:),Freal(:,:),Fimage(:,:)
REAL,ALLOCATABLE::U(:),v(:),w(:,:),Fcount_Real(:,:),Fcount_image(:,:)
cmd_filename='wave.par'
CALL Read_CMD(cmd_filename,height,input_filename,Expound_a20_2D_filename,output_filename,eigval)
! logical unit_open
! integer nunit_in,nunit_start
! nunit_in=nunit_start
! unit_open=.TRUE.
! do while (unit_open)
! nunit_in=nunit_in+1
! inquire(unit=nunit_in,opened=unit_open)
END SUBROUTINE
!……………………
!计算扩边点位
!……………………
SUBROUTINE Calculate_m1m2_dx(mpoint,m0,m1,m2,m3,m,Xmin,Xmax,dx)
integer mpoint,mtemp,m0,m1,m2,m3,m
real Xmin,Xmax,dx,factor_m,mu
ALLOCATE(U(m0:m3),V(n0:n3),W(m0:m3,n0:n3),Fcount_Real(m0:m3,n0:n3),Fcount_image(m0:m3,n0:n3))
ALLOCATE(Field(m0:m3,n0:n3),Freal(m0:m3,n0:n3),Fimage(m0:m3,n0:n3))
mu=int(log(float(mpoint))/0.693147+factor_m)
m=2**mu
end if
m1=1+(m-mpoint)/2
m2=m1+mpoint-1
m3=m
END SUBROUTINE
!………………………………
!输入源文件的重磁异常值
!………………………………
SUBROUTINE Input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field,eigval)
! call search_unit(10,nunit_in)
OPEN(12,file=input_filename)
READ(12,*)
READ(12,*) mpoint,line
READ(12,*) Xmin,Xmax
READ(12,*) Ymin,Ymax
! read(12,*)
close(12)
Deallocate(Field,Freal,Fimage,U,V,W,Fcount_Real,Fcount_image)
END PROGRAM
!……………………………………
!搜寻未打开的文件的通道号
!……………………………………
!SUBROUTINE SEARCH_UNIT(NUNIT_START,NUNIT_IN)
Freal=Field
Fimage=0.
CALL FFt2(Freal,Fimage,m3-m0+1,n3-n0+1,2) !2说明是正变换
CALL WAVE2D_sub(U,V,W,m,n,dx,dy)
CALL Factor_conti(height,m0,m3,n0,n3,U,V,W,Fcount_Real,Fcount_image) !计算延拓因子
CALL Multi_sub(Freal,Fimage,Fcount_Real,Fcount_image,m0,m3,n0,n3)
CALL FFt2(Freal,Fimage,m3-m0+1,n3-n0+1,1)
CALL Putout_GRD(Freal,output_filename,m0,m1,m2,m3,n0,n1,n2,n3,eigval,Xmin,Xmax,Ymin,Ymax)
m0=1
factor_m=1.5
dx=(Xmax-Xmin)/(mpoint-1)
mtemp=mpoint
do while((mod(mtemp,2).eq.0).and.(mtemp.ne.0))
mtemp=mtemp/2
end do
if(mtemp.eq.1) then
m=mpoint*2
else
CALL Read_mn(input_filename,mpoint,line,Xmin,Xmax,Ymin,Ymax)
CALL Calculate_m1m2_dx(mpoint,m0,m1,m2,m3,m,Xmin,Xmax,dx)
CALL Calculate_m1m2_dx(line,n0,n1,n2,n3,n,Ymin,Ymax,dy)
read(11,*) height
read(11,*) input_filename
read(11,*) Expound_a20_2D_filename
read(11,*) output_filename
read(11,*) eigval
close(11)
END SUBROUTINE