实验四位场边缘识别程序设计实验
边缘化提取实验报告
一、实验目的1. 理解图像边缘检测的基本原理和过程。
2. 掌握常用的边缘检测算法，如Roberts算子、Sobel算子、Prewitt算子、Laplacian算子和Canny算子。
3. 通过实验验证不同边缘检测算法的效果，并分析其优缺点。
4. 了解特征提取的基本原理和方法，对图像边缘进行特征提取。
二、实验原理图像边缘是图像中灰度值或颜色值发生突变的地方，是图像分割和特征提取的基础。
边缘检测的目的是找到图像中灰度值变化明显的区域，即边缘。
边缘检测算法可以分为两类：基于微分算子的边缘检测算法和基于二值化的边缘检测算法。
1. 基于微分算子的边缘检测算法：- 利用一阶导数或二阶导数检测图像边缘。
- 常见的算子有Roberts算子、Sobel算子、Prewitt算子、Laplacian算子等。
2. 基于二值化的边缘检测算法：- 利用图像的二值化处理，将图像分为前景和背景两部分。
- 常见的算法有Otsu算法、Sauvola算法等。
三、实验内容1. 实验材料：- OpenCV库- Python编程环境2. 实验步骤：（1）读取图像：使用OpenCV库读取待检测的图像。
（2）灰度化：将图像转换为灰度图像，以便进行边缘检测。
（3）边缘检测：- 使用Roberts算子检测边缘。
- 使用Sobel算子检测边缘。
- 使用Prewitt算子检测边缘。
- 使用Laplacian算子检测边缘。
- 使用Canny算子检测边缘。
（4）特征提取：对检测到的边缘进行特征提取，如计算边缘长度、宽度、方向等。
（5）结果展示：将检测到的边缘和提取的特征进行可视化展示。
四、实验结果与分析1. Roberts算子：- 效果：Roberts算子对图像噪声敏感，边缘检测效果较差。
- 分析：Roberts算子对图像局部区域进行检测，容易受到噪声的影响。
2. Sobel算子：- 效果：Sobel算子对图像噪声有一定的抑制能力，边缘检测效果较好。
- 分析：Sobel算子使用高斯滤波器对图像进行平滑处理，然后计算图像的一阶导数。
机器视觉应用技术实验05图像边缘检测
实验5 图像边缘检测一、实验目的1.掌握图像边缘检测的方法。
2.掌握AiCam框架的部署和使用。
二、实验环境硬件环境：PC机Pentium处理器双核2GHz以上，内存4GB以上操作系统：Windows7 64位及以上操作系统开发软件：MobaXterm实验器材：人工智能边缘应用平台实验配件：无三、实验内容1.算法原理1.1 基本描述边缘检测是图像处理和计算机视觉中的基本问题，边缘检测的目的是标识数字图像中亮度变化明显的点。
图像属性中的显著变化通常反映了属性的重要事件和变化。
图像边缘检测大幅度地减少了数据量，并且剔除了可以认为不相关的信息，保留了图像重要的结构属性。
本实验中使用的是canny边缘检测算子，除此之外还有Sobel、Laplacian算子等。
1.2 专业术语Canny边缘检测Canny边缘检测器由John F. Canny于1986年开发，是当今最流行的边缘检测方法之一，因为它非常强大和灵活。
该算法主要处理过程如下：1）降噪：原始图像像素通常会导致噪声边缘，因此在计算边缘之前减少噪声很重要在Canny边缘检测中，高斯模糊过滤器用于从本质上去除或最小化可能导致不良边缘的不必要细节。
2）计算图像的强度梯度：图像平滑（模糊）后，将使用Sobel核心进行水平和垂直过滤。
然后使用这些过滤操作的结果来计算每个像素的强度梯度3）抑制假边：在降低噪声和计算强度梯度后，这一步的算法使用一种称为边缘非最大抑制的技术来过滤掉不需要的像素（实际上可能并不构成边缘）。
为此，在正负梯度方向上将每个像素与其相邻的像素进行比较。
如果当前像素的梯度幅度大于其相邻像素，则保持不变。
否则，当前像素的大小设置为零。
4）滞后阈值：在Canny边缘检测的最后一步中，梯度幅度与两个阈值进行比较。
1.3 常用方法OpenCV实现Canny图像边缘检测方法如下：# image：原始图像# threshold1：阈值1 (minVal)# threshold2：阈值2 (maxVal)# 推荐的高低阈值比在2:1到3:1之间edges = cv2.Canny(image, threshold1, threashold2)2.功能设计2.1 功能描述AiCam人工智能轻量化应用框架是一款面向于人工智能边缘应用的开发框架，采用统一模型调用、统一硬件接口、统一算法封装和统一应用模板的设计模式，实现了嵌入式边缘计算环境下进行快速的应用开发和项目实施。
AVR学习笔记十九、4X4矩阵键盘实验
A VR学习笔记十九、4X4矩阵键盘实验19.1 实例功能在前面的实例中我们已经学习了在单片机系统中检测独立式按键的接口电路和程序设计，独立式按键的每个按键占用1位I/O口线，其状态是独立的，相互之间没有影响，只要单独测试链接案件的I/O口线电平的高低就能判断键的状态。
独立式按键电路简单、配置灵活，软件结构也相对简单。
此种接口方式适用于系统需要按键数目较少的场合。
在按键数量较多的情况下，如系统需要8个以上按键的键盘时，采用独立式接口方式就会占用太多的I/O口，这对于I/O口资源不太丰富的单片机系统来说显得相当浪费，那么当按键数目相对较多的时候，为了减少I/O口资源的占用，应该采取什么样的方式才能够既满足多按键识别，又减少I/O口的占用呢？当然我们可以采用端口扩展器件比如串并转换芯片实现单片机I/O口的扩展，但是这种方式既增加了电路的复杂性，又增加了系统的成本开销。
有没有比较经济实惠的方法呢？事实上，在实际引用中我们经常采用矩阵式键盘的方式来节约I/O口资源和系统成本。
在这个实验中，我们采用4X4矩阵键盘来实现使用8个I/O口识别16个按键的实验，本实例分为三个功能模块，分别描述如下：●单片机系统：利用A Tmega16单片机与矩阵键盘电路实现多按键识别。
●外围电路：4X4矩阵键盘电路、LED数码管显示电路。
●软件程序：编写软件，实现4X4矩阵键盘识别16个按键的程序。
通过本实例的学习，掌握以下内容：●4X4矩阵键盘的电路设计和程序实现。
19.2 器件和原理19.2.1 矩阵键盘的工作原理和扫描确认方式当键盘中按键数量较多时，为了减少对I/O口的占用，通常将按键排列成矩阵形式，也称为行列键盘，这是一种常见的连接方式。
矩阵式键盘接口见图1所示，它由行线和列线组成，按键位于行、列的交叉点上。
当键被按下时，其交点的行线和列线接通，相应的行线或列线上的电平发生变化，MCU通过检测行或列线上的电平变化可以确定哪个按键被按下。
使用计算机视觉技术进行图像边缘检测的步骤和注意事项
使用计算机视觉技术进行图像边缘检测的步骤和注意事项计算机视觉技术是一门研究如何使机器“看见”并理解图像或视频的技术。
其中一项重要的任务是图像边缘检测。
图像边缘是图像中像素灰度值变化明显的区域，边缘检测是在图像中找到这些边缘的过程。
本文将介绍使用计算机视觉技术进行图像边缘检测的步骤和注意事项。
图像边缘检测的步骤通常包括以下几个关键步骤：1. 预处理：首先，对输入的图像进行预处理。
预处理的目的是消除噪声、增强图像的对比度，以便更好地检测边缘。
常用的预处理方法包括高斯滤波、中值滤波和直方图均衡化等。
2. 灰度转换：将彩色图像转换为灰度图像。
这是因为大多数边缘检测算法在灰度图像上运行。
可以使用加权平均法或者取红、绿、蓝三个通道的平均值的方法将彩色图像转换为灰度图像。
3. 计算梯度：通过计算图像中每个像素点的梯度来确定边缘的位置。
梯度指的是图像灰度值的变化程度。
常用的方法有Sobel、Prewitt和Laplacian等算子。
这些算子可以检测水平、垂直和对角线方向上的边缘。
4. 非极大值抑制：在计算梯度之后，可能会出现多个边缘候选点。
非极大值抑制的目的是在提取出的边缘候选点中选取局部最大值，以得到更准确的边缘。
5. 双阈值处理和边缘连接：通过设置合适的阈值将边缘分为强边缘和弱边缘。
强边缘即明显的边缘，而弱边缘则可能是噪声或非边缘。
通常选择两个阈值进行分割，边缘像素灰度值大于高阈值的被标记为强边缘，灰度值介于低阈值和高阈值之间的被标记为弱边缘。
然后可以使用边缘连接的方法将弱边缘连接到强边缘，得到完整的边缘。
6. 后处理：根据应用需求进行后处理，如边缘修复、边缘精化等。
在进行图像边缘检测时，还需要注意以下几个事项：1. 选择合适的边缘检测算法：根据不同应用的需求选择适合的边缘检测算法。
常用的边缘检测算法包括Canny算法、Sobel算子、Laplacian算子等。
2. 调整算法参数：不同的边缘检测算法有不同的参数需调整。
图像边缘检测各种算子MATLAB实现以及实际应用
《图像处理中的数学方法》实验报告学生姓名：***教师姓名：曾理学院：数学与统计学院专业：信息与计算科学学号：********联系方式：139****1645梯度和拉普拉斯算子在图像边缘检测中的应用一、数学方法边缘检测最通用的方法是检测灰度值的不连续性，这种不连续性用一阶和二阶导数来检测。
1.（1）一阶导数：一阶导数即为梯度，对于平面上的图像来说，我们只需用到二维函数的梯度，即：∇f=[g xg y]=[ðf ðxðfðy],该向量的幅值：∇f=mag(∇f)=[g x2+g y2]1/2= [(ðf/ðx)2+(ðf/ðy)2]1/2,为简化计算，省略上式平方根，得到近似值∇f≈g x2+g y2；或通过取绝对值来近似，得到：∇f≈|g x|+|g y|。
（2）二阶导数：二阶导数通常用拉普拉斯算子来计算，由二阶微分构成：∇2f(x,y)=ð2f(x,y)ðx2+ð2f(x,y)ðy22.边缘检测的基本思想：（1）寻找灰度的一阶导数的幅度大于某个指定阈值的位置；（2）寻找灰度的二阶导数有零交叉的位置。
3.几种方法简介（1）Sobel边缘检测器：以差分来代替一阶导数。
Sobel边缘检测器使用一个3×3邻域的行和列之间的离散差来计算梯度，其中，每行或每列的中心像素用2来加权，以提供平滑效果。
∇f=[g x2+g y2]1/2={[(z7+2z8+z9)−(z1+2z2+z3)]2+[(z3+2z6+z9)−(z1+2z4+z7)]2}1/2（2）Prewitt边缘检测器:使用下图所示模板来数字化地近似一阶导数。
与Sobel检测器相比，计算上简单一些，但产生的结果中噪声可能会稍微大一些。
g x=(z7+z8+z9)−(z1+z2+z3)g y=(z3+z6+z9)−(z1−z4−z7)（3）Roberts边缘检测器：使用下图所示模板来数字化地将一阶导数近似为相邻像素之间的差，它与前述检测器相比功能有限（非对称，且不能检测多种45°倍数的边缘）。
实验四位场边缘识别程序设计实验
《重磁资料处理与解释》实验四位场边缘识别程序设计实验专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-1-31、基本原理地质目标体边缘时指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线，在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。
现在有重、磁位场边缘识别方法分为数理统计、数值计算以及其他三大类。
数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征。
在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值。
故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置，确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置和真实的位置有一定的偏差。
该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化。
因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。
(1)垂向导数：垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可以直接使用，对磁力异常必须转化为磁源重力异常或化极磁力异常才可以使VD R x y z：g(x, y, z)用VD Rx y去‘az ( 1.1)(2)解析信号振幅：解析信号振幅也是利用极大值位置来确定地质体的边缘位置适用于重、磁力异常ASM二THDR2 VDR2(3)总水平导数(THDR)2、输入/输出数据格式设计依据上述原理，现在对上述各种边缘识别方法进行程序设计。
(1-2)THDR(x, y,z)(1-3)2.1 输入数据格式设计2.3 参数文件数据格式设计将以上部分量保存在一个文件中，该文件名变量为 cmdfile ，字符串变量， 长度不超过 80，全路径名。
在该文件中保存的参数如下：输入数据文件名 input_file ，字符串变量，长度不超过 80；输出 vdr 数据文件名 output_file_vdr ，字符串变量，长度不超过 80； 输出 thdr 数据文件名 output_file_thdr ，字符串变量，长度不超过 80; 输出asm 数据文件名output_file_asm ，字符串变量，长度不超过 80factor_m ：扩边比例因子，实型变量（＞1）;本次实验给了正演的重力异常数据，为例如： .GRD 格式，均为实型变量DSAA201 -1000.000000 -1000.000000 5.549671E-01 5.549671E-01 5.897312E-012011000.000000 1000.000000 23.539846 5.634658E-01 5.987253E-015.721339E-01 5.808522E-016.078691E-01 6.171604E-012.2 输出数据格式设计 计算结果输出数据格式与输入格式对应，例如： 格式为.GRD 格式，均为实型变量DSAA201 -1000.000 -1000.000 -0.1465084 -3.6523044E-02 -3.2688729E-02 -3.2654848E-02-3.3138681E-02201 1000.000 1000.000 0.3190881 -3.3485338E-02 -3.2723978E-02-3.3061244E-02 -3.2748722E-02 -3.2787599E-02 -3.2927759E-023.总体设计此次程序采用IPO结构设计，首先通过读取cmd文件，得到相关输入参数：输入数据文件名gravity.grd、输出vdr文件名field_vrd.grd、输出thdr文件名field_thdr.grd、输出asm文件名field_asm.grd、扩边比例因子factor_m ;然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。
实验四 四位加法器的设计
实验四四位加法器的设计一、实验目的1、掌握VHDL程序设计流程。
二、实验内容设计四位加法器，并在数码管上显示计算结果。
三、实验仪器1、ZY11EDA13BE型实验箱通用编程模块，配置模块，开关按键模块，数码显示模块。
2、并口延长线，JTAG延长线。
3、安装MAX+PLUSII 10.2软件的PC机。
四、实验原理用VHDL编辑四位加法器程序，用拨码开关输入，加数、被加数、进位，用数码管显示计算结果和，用发光二极管显示进位。
五、实验步骤：步骤1：输入VHDL程序，编译，仿真，锁定引脚并下载到目标芯片。
步骤2：验证设计结果。
六、实验报告1、列出VHDL源程序。
七、思考题1、怎样实现8位加法器？2、对于数码管显示，如何实现循环显示不同数字？library ieee;use ieee.std_logic_1164.all;use ieee.std_logic_unsigned.all;use ieee.std_logic_arith.all;entity add8 isport(a,b:in std_logic_vector(3 downto 0);a0,b0,c0,c1:in std_logic;output:out std_logic_vector(6 downto 0);c2:out std_logic);end add8;architecture art of add8 issignal e:std_logic_vector(4 downto 0);beginprocess(a,b,c)begine<=a+b+c1;c2<=e(4);case e(3 downto 0) iswhen "0000" => output <="1111110";when "0001" => output <="0110000";when "0010" => output <="1101101";when "0011" => output <="1111001";when "0100" => output <="0110011";when "0101" => output <="1011011";when "0110" => output <="1011111";when "0111" => output <="1110000";when "1000" => output <="1111111";when "1001" => output <="1111011";when "1010" => output <="1110111";when "1011" => output <="0011111";when "1100" => output <="0001101";when "1101" => output <="0111101";when "1110" => output <="1001111";when "1111" => output <="1000111";when others => output <="0000000"; end case;a0 <='1';b0 <='1';c0<='1';end process;end art;entity scan_led isport(clk : in std_logic;sele : out std_logic_vector(1 downto 0);time1,time2 : out std_logic_vector(2 downto 0)); end;architecture a of scan_led issignal a : std_logic_vector(1 downto 0):="00";beginpro1:process(clk)beginif clk'event and clk='1' then a<=a+1; end if;end process pro1;pro2:process(a)begincase a iswhen "00" => time1<="111";time2<="001";when "01" => time1<="110";time2<="010";when "10" => time1<="101";time2<="100";when "11" => time1<="011";time2<="110";when others => null;end case;sele<=a;end process pro2;end a;。
Canny边缘检测算法的流程
Canny边缘检测算法的流程介绍边缘检测的⼀般标准包括：1) 以低的错误率检测边缘，也即意味着需要尽可能准确的捕获图像中尽可能多的边缘。
2) 检测到的边缘应精确定位在真实边缘的中⼼。
3) 图像中给定的边缘应只被标记⼀次，并且在可能的情况下，图像的噪声不应产⽣假的边缘。
在⽬前常⽤的边缘检测⽅法中，Canny边缘检测算法是具有严格定义的，可以提供良好可靠检测的⽅法之⼀。
由于它具有满⾜边缘检测的三个标准和实现过程简单的优势，成为边缘检测最流⾏的算法之⼀。
Canny边缘检测算法的处理流程Canny边缘检测算法可以分为以下5个步骤：1) 使⽤⾼斯滤波器，以平滑图像，滤除噪声。
2) 计算图像中每个像素点的梯度强度和⽅向。
3) 应⽤⾮极⼤值（Non-Maximum Suppression）抑制，以消除边缘检测带来的杂散响应。
4) 应⽤双阈值（Double-Threshold）检测来确定真实的和潜在的边缘。
5) 通过抑制孤⽴的弱边缘最终完成边缘检测。
下⾯详细介绍每⼀步的实现思路。
1 ⾼斯平滑滤波为了尽可能减少噪声对边缘检测结果的影响，所以必须滤除噪声以防⽌由噪声引起的错误检测。
为了平滑图像，使⽤⾼斯滤波器与图像进⾏卷积，该步骤将平滑图像，以减少边缘检测器上明显的噪声影响。
⼤⼩为(2k+1)x(2k+1)的⾼斯滤波器核的⽣成⽅程式由下式给出：下⾯是⼀个sigma = 1.4，尺⼨为3x3的⾼斯卷积核的例⼦（需要注意归⼀化）：若图像中⼀个3x3的窗⼝为A，要滤波的像素点为e，则经过⾼斯滤波之后，像素点e的亮度值为：其中*为卷积符号，sum表⽰矩阵中所有元素相加求和。
重要的是需要理解，⾼斯卷积核⼤⼩的选择将影响Canny检测器的性能。
尺⼨越⼤，检测器对噪声的敏感度越低，但是边缘检测的定位误差也将略有增加。
⼀般5x5是⼀个⽐较不错的trade off。
2 计算梯度强度和⽅向图像中的边缘可以指向各个⽅向，因此Canny算法使⽤四个算⼦来检测图像中的⽔平、垂直和对⾓边缘。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《重磁资料处理与解释》实验四位场边缘识别程序设计实验专业名称：地球物理学学生姓名：学生学号：指导老师：王万银、纪新林、纪晓琳、邱世灿提交日期：2016-1-31、基本原理地质目标体边缘时指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线，在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。
现在有重、磁位场边缘识别方法分为数理统计、数值计算以及其他三大类。
数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征。
在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值。
故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置，确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置和真实的位置有一定的偏差。
该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化。
因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。
(1)垂向导数：垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可以直接使用，对磁力异常必须转化为磁源重力异常或化极磁力异常才可以使用(,,)(,,)g x y zV D R x y zz∂=∂( 1.1)（2）解析信号振幅：解析信号振幅也是利用极大值位置来确定地质体的边缘位置适用于重、磁力异常(1.2)（3）总水平导数（THDR）(1.3)2、输入/输出数据格式设计依据上述原理，现在对上述各种边缘识别方法进行程序设计。
2.1 输入数据格式设计本次实验给了正演的重力异常数据，为.GRD格式，均为实型变量。
例如：DSAA201 201-1000.000000 1000.000000-1000.000000 1000.0000005.549671E-01 23.5398465.549671E-01 5.634658E-01 5.721339E-01 5.808522E-015.897312E-01 5.987253E-016.078691E-01 6.171604E-01…2.2 输出数据格式设计计算结果输出数据格式与输入格式对应，格式为.GRD格式，均为实型变量。
例如：DSAA201 201-1000.000 1000.000-1000.000 1000.000-0.1465084 0.3190881-3.6523044E-02 -3.3485338E-02 -3.3061244E-02 -3.2748722E-02 -3.2688729E-02-3.2654848E-02 -3.2723978E-02 -3.2787599E-02 -3.2927759E-02 -3.3138681E-02…2.3 参数文件数据格式设计将以上部分量保存在一个文件中，该文件名变量为cmdfile，字符串变量，长度不超过80，全路径名。
在该文件中保存的参数如下：输入数据文件名input_file，字符串变量，长度不超过80；输出vdr数据文件名output_file_vdr，字符串变量，长度不超过80；输出thdr数据文件名output_file_thdr，字符串变量，长度不超过80;输出asm数据文件名output_file_asm，字符串变量，长度不超过80factor_m：扩边比例因子，实型变量（>1）；3.总体设计此次程序采用IPO结构设计，首先通过读取cmd文件，得到相关输入参数：输入数据文件名gravity.grd、输出vdr文件名field_vrd.grd、输出thdr文件名field_thdr.grd、输出asm文件名field_asm.grd、扩边比例因子factor_m；然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。
下一步，进行二维余弦扩边，将扩完边的数据进行快速二维傅里叶变换，转换到频率域；接下来在频率域求出在x,y,z方向的导数并反变换；最后求出VDR、THDR、ASM数据。
最后去除扩边部分后输出。
总体设计见表1。
4.测试结果分析：由图4.3可看出，VDR 方法的零值线可较准确识别模型体的边界；由图4.4可看出，ASM 的极大值点边界可大致识别模型体边界，但精度不是很高。
对比图4.2到4.5可以看出，THDR 方法对模型边界的识别效果是最好的。
5 结论及建议经测试，VDR 与THDR 对模型体的边界位置识别效果较好，而ASM 对模型体边界识别效果较差。
三种方法中，THDR 效果最好。
图4.3 垂向导数（VDR ）图4.2 重力异常原始图附录：边缘识别程序源代码!******************************************************************************************* *********! 程序功能：实现频率域位场导数运算进行边缘识别! cmd文件参数:! cmdfile:存放有关参数的文件名变量! input_file:观测面位场数据文件! output_file_vdr：场值的水平导数数据文件! output_file_thdr:场值垂向导数数据文件! output_file_asm:场值的解析信号振幅数据文件! factor_m：扩边因子! .grd文件参数：! N_point,N_line：点数(x方向)、线数(y方向)! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! Ur:初始观测面场值! 扩边参数：! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! 求导参数：! field_re:初始观测面信号的实部! field_im：初始观测面信号的虚部! Px_re:x方向导数信号的实部! Px_im:x方向导数信号的虚部! Py_re:y方向导数信号的实部! Py_im:y方向导数信号的虚部! Pz_re:z方向导数信号的实部! Pz_im:z方向导数信号的虚部! W(m,n)：径向圆频率! field_vdr:对场值作水平导数的结果! field_thdr:对场值作垂向导数的结果! field_asm:场值的解析信号振幅的结果!******************************************************************************************* *********program deviationparameter(eigval=3.701411*1e5)character*(80)cmdfilecharacter*80 input_file,output_file_vdr,output_file_thdr,output_file_asmreal,allocatable:: field_re(:,:),field_im(:,:)real,allocatable:: Px_re(:,:),Px_im(:,:),Py_re(:,:),Py_im(:,:),Pz_re(:,:),Pz_im(:,:)real,allocatable:: field_vdr(:,:),field_thdr(:,:),field_asm(:,:)real,allocatable:: U(:),V(:),W(:,:)integer N_point,N_lineinteger m,n,m1,m2,n1,n2real factor_mreal xmin,xmax,ymin,ymax,dx,dycmdfile='deviation.cmd'call read_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)call read_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)call calculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(field_re(1:m,1:n),field_im(1:m,1:n))allocate(Px_re(1:m,1:n),Px_im(1:m,1:n),Py_re(1:m,1:n),Py_im(1:m,1:n),Pz_re(1:m,1:n),Pz_im(1:m,1:n)) allocate(field_vdr(1:m,1:n),field_thdr(1:m,1:n),field_asm(1:m,1:n))allocate(U(1:m),V(1:n),W(1:m,1:n))call input_grd(field_re,input_file,m1,m2,n1,n2,m,n)call expand_cos_2D(m1,m2,m,n1,n2,n,field_re,field_im)call FFT2(field_re,field_im,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,U,V,W)call factor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)call FFT2(px_re,px_im,m,n,1)call FFT2(py_re,py_im,m,n,1)call FFT2(pz_re,pz_im,m,n,1)call deviration(m,n,px_re,py_re,pz_re,field_vdr,field_thdr,field_asm)call OUTPUT_GRD(field_vdr,output_file_vdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)call OUTPUT_GRD(field_thdr,output_file_thdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)call OUTPUT_GRD(field_asm,output_file_asm,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)deallocate(field_re,field_im,px_re,px_im,py_re,py_im,pz_re,pz_im,field_vdr,field_thdr,field_asm,u,v,w) end program!****************************************************************************! 子程序：read_cmd! 功能：读取参数文件! 输入参数说明：! cmdfile：参数文件名! 输出参数说明：! input_file:观测面位场数据文件! output_file_vdr：对场值作水平导数处理后的数据文件! output_file_thdr:对场值作垂向导数处理后的数据文件! output_file_asm:对场值作总导数处理后的数据文件! factor_m：扩边因子!**************************************************************************** Subroutine read_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm) implicit nonecharacter*80 strcharacter*(*)cmdfilecharacter*(*)input_file,output_file_vdr,output_file_thdr,output_file_asmreal factor_mopen(10,file=cmdfile,status='old')read(10,*)str,input_fileread(10,*)str,output_file_vdrread(10,*)str,output_file_thdrread(10,*)str,output_file_asmread(10,*)str,factor_mclose(10)end Subroutine read_cmd!***************************************************************************! 子程序：read_grd! 功能：从原始观测.grd文件中读取相关参数! 输入参数说明：! filename_obser：输入文件名! 输出参数说明：! N_point,N_line：点数、线数! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值!*************************************************************************** subroutine read_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)implicit nonecharacter*(*)input_fileinteger N_point,N_linereal Xmin,Xmax,Ymin,Ymaxopen(10,file=input_file,status='old')Read(10,*)Read(10,*)N_line,N_pointRead(10,*)Xmin,XmaxRead(10,*)Ymin,YmaxClose(10)end subroutine read_grd!************************************************************************** ! 子程序：calculate_mn! 功能：确定扩边数据点号位置! 输入参数说明：! factor_m: 扩边比例因子（>1.0）! a,b：点数、线数! 输出参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)!**************************************************************************subroutine calculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicit noneinteger a,b,m,n,m1,m2,n1,n2integer mtemp,mu,nureal factor_mmtemp=aDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENm=a*2ELSEmu=int(log(float(a))/0.693147+factor_m)m=2**muEND IFm1=1+(m-a)/2m2=m1+a-1write(*,*)m,apausemtemp=bDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENn=b*2ELSEnu=int(log(float(b))/0.693147+factor_m)n=2**nuEND IFn1=1+(n-b)/2n2=n1+b-1write(*,*)m1,m2,n1,n2,m,npauseend subroutine calculate_mn!************************************************************************* ! 子程序：INPUT_GRD! 功能：读取grd文件中的数据! 输入参数说明：! filename_obser：输入文件名! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! 输出参数说明：! A：存放输出数据的二维数组名!************************************************************************* SUBROUTINE INPUT_GRD(A,input_file,m1,m2,n1,n2,m,n)character*(*)input_fileinteger m1,m2,n1,n2,m,nreal xmin,xmax,ymin,ymaxreal A(1:m,1:n)real i,j,kdo j=1,n,1do i=1,m,1A(i,j)=3.701411*1e10enddoenddoOpen(20,file=input_file,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)END SUBROUTINE INPUT_GRD!************************************************************************* ! 子程序：expand_cos_2D! 功能：二维扩边子程序并为信号虚部赋值! 输入参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部! 输出参数说明：! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部!*************************************************************************subroutine expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicit noneinteger m,n,m1,m2,n1,n2real Ur(1:m,1:n),Ui(1:m,1:n)real,allocatable::u(:),r(:)integer j,i,kallocate(u(1:m))do j=n1,n2,1do i=1,m,1u(i)=Ur(i,j)enddocall expand_cos_1d(1,m1,m2,m,u(1))do i=1,m,1Ur(i,j)=u(i)enddoenddodeallocate(u)allocate(r(1:n))do i=1,m,1do j=1,n,1r(j)=Ur(i,j)enddocall expand_cos_1d(1,n1,n2,n,r(1))do j=1,n,1Ur(i,j)=r(j)enddoenddodeallocate(r)do i=1,mdo j=1,nUi(i,j)=0enddoenddoend subroutine expand_cos_2D!************************************************************************** ! 子程序：expand_cos_1d! 功能：一维扩边子程序! 输入参数说明：! n0,n3：扩边后数据起点位置和终点位置! n1,n2：实际数据起点位置和终点位置! feild(i),(i=n1,n1+1,...,n2)：实际数据! 输出参数说明：! field(i),(i=n0,...,n1-1)：扩边后左边的数据! field(i),(i=n2+1,...,n3)：扩边后右边的数据!************************************************************************** Subroutine expand_cos_1d(n0,n1,n2,n3,Field)Real Field(n0:n3)pi=3.141592654Field (n0)=(Field (n1)+Field (n2))/2.0Field (n3)=Field (n0)i1=n0i2=n1DO i=i1,i2-1,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doi1=n2i2=n3DO i=i1+1,i2,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doEnd subroutine expand_cos_1d!************************************************************************ ! 功能：FFT2! 功能：复数组2-D快速Fourier变换! 输入参数说明：! m0,m3：x方向的起点和终点! n0,n3：y方向的起点和终点! field：输入信号(需要赋值给Freal，实部)! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Freal：信号的实部! Fimage：信号的虚部(对于实信号而言，赋值为0)! 对应频率分布说明：! 数据Freal(m,n)和Fimage(m,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!************************************************************************ SUBROUTINE FFT2(Freal,Fimage,m,n,nf)implicit noneINTEGER m,n,nfREAL Freal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE::Treal(:),Timage(:)integer nmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0) STOPDO i=1,m,1IF (n.ne.1) THENdo j=1,n,1Treal(j)=Freal(i,j)Timage(j)=Fimage(i,j)end docall FFT(Treal,Timage,n,nf)do j=1,n,1Freal(i,j)=Treal(j)Fimage(i,j)=Timage(j)end doEND IFEND DODO j=1,n,1IF(m.ne.1) THENdo i=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)end docall FFT(Treal,Timage,m,nf)do i=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)end doEND IFEND DOdeallocate(Treal,Timage,STAT=ierr)END SUBROUTINE FFT2!****************************************************************** ! 子程序：FFT! 功能：复数组1-D快速Fourier变换! 输入参数说明：! Xreal(n): 输入数据实部! Ximage(n): 输入数据虚部! N: 点数（N 必须为2 的幂次方）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Xreal(n): 变换后频谱实部! Ximage(n): 变换后频谱虚部! 对应频率分布说明：! 数据Xreal(n)和Ximage(n)对应的频率分布位置为:! 0,1,.......,n/2-1,n/2,-(n/2-1),......,-1!***************************************************************** SUBROUTINE FFT(Xreal,Ximage,n,nf)implicit noneINTEGER n,nfREAL Xreal(1:n),Ximage(1:n)integer nu,n2,nu1,k,k1,k1n2,l,i,ibitrreal f,p,arg,c,s,treal,timagenu=int(log(float(n))/0.693147+0.001)n2=n/2nu1=nu-1f=float((-1)**nf)k=0DO l=1,nu,1DO while (k.lt.n)do i=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+1end dok=k+n2END DOk=0nu1=nu1-1n2=n2/2END DODO k=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)=timageend ifEND DOIF(nf.ne.1) THENdo i=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)end doEND IFEND SUBROUTINE FFTFUNCTION IBITR(J,NU)implicit noneinteger ibitr,j,nuinteger j1,itt,i,j2j1=jitt=0do i=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2end doEND FUNCTION IBITR!*************************************************************************** ! 子程序：cal_dxdy! 功能：计算x,y方向的点距! 输入参数说明：! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! N_point,N_line：点数(x方向)、线数(y方向)! 输出参数说明：! dx,dy：x,y方向的点距!*************************************************************************** subroutine cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicit nonereal xmin,xmax,ymin,ymaxinteger N_POINT,N_LINEreal dx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEend subroutine cal_dxdy!******************************************************************! 子程序：WAVE2D! 功能：计算2D径向圆频率W! 输入参数说明：! dx: x 方向点距! dy: y 方向线距! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! 输出参数说明：! W(m,n): 径向圆频率!******************************************************************SUBROUTINE WAVE2D(m,n,dx,dy,U,V,W)implicit noneINTEGER m,nREAL dx,dyREAL W(1:m,1:n),U(1:m),V(1:n)real pi,delx,delyinteger midm,midn,i,j,xx,yymidm=m/2+1midn=n/2+1delx=float(m)/dxdely=float(n)/dydo j=1,n,1yy=jif(yy.gt.midn) yy=yy-nv(j)=dely*(yy-1)end dodo i=1,m,1xx=iif(xx.gt.midm) xx=xx-mu(i)=delx*(xx-1)end dodo j=1,n,1do i=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j))end doend doEND SUBROUTINE WAVE2D!*********************************************************************************** ! 子程序：factor! 功能：计算x,y,z方向导数的实部和虚部!! 输入参数说明：! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! field_re:初始观测面信号的实部! field_im：初始观测面信号的虚部! W(m,n)：径向圆频率! 输出参数说明：! Px_re:x方向导数信号的实部! Px_im:x方向导数信号的虚部! Py_re:y方向导数信号的实部! Py_im:y方向导数信号的虚部! Pz_re:z方向导数信号的实部! Pz_im:z方向导数信号的虚部!*********************************************************************************** subroutine factor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)implicit noneinteger m,nreal field_re(1:m,1:n),field_im(1:m,1:n)real px_re(1:m,1:n),px_im(1:m,1:n),py_re(1:m,1:n),py_im(1:m,1:n),pz_re(1:m,1:n),pz_im(1:m,1:n) real U(1:m),V(1:n),W(1:m,1:n)integer i,jreal pipi=3.1415926do i=1,m,1do j=1,n,1px_re(i,j)=field_im(i,j)*u(i)*(-1)*pi*2.0px_im(i,j)=field_re(i,j)*u(i)*pi*2.0py_re(i,j)=field_im(i,j)*v(j)*(-1)*pi*2.0py_im(i,j)=field_re(i,j)*v(j)*pi*2.0pz_re(i,j)=field_re(i,j)*w(i,j)*pi*2.0pz_im(i,j)=field_im(i,j)*w(i,j)*pi*2.0end doend doend subroutine factor!***************************************************************************! 子程序：deviration! 功能：计算异常的水平导数、垂向导数及总导数!! 输入参数说明：! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! Px_re:x方向导数信号的实部! Px_im:x方向导数信号的虚部! Py_re:y方向导数信号的实部! Py_im:y方向导数信号的虚部! Pz_re:z方向导数信号的实部! Pz_im:z方向导数信号的虚部! 输出参数说明：! field_vdr:对场值作水平导数的结果! field_thdr:对场值作垂向导数的结果! field_asm:对场值作总导数的结果!************************************************************************** subroutine deviration(m,n,px_re,py_re,pz_re,field_vdr,field_thdr,field_asm)implicit noneinteger m,nreal px_re(1:m,1:n),py_re(1:m,1:n),pz_re(1:m,1:n)real field_vdr(1:m,1:n),field_thdr(1:m,1:n),field_asm(1:m,1:n)integer i,jdo i=1,m,1do j=1,n,1field_vdr(i,j)=pz_re(i,j)field_thdr(i,j)=sqrt(px_re(i,j)**2+py_re(i,j)**2)field_asm(i,j)=sqrt(field_vdr(i,j)**2+field_thdr(i,j)**2)end doend doend subroutine deviration!*************************************************************************** ! 子程序：OUTPUT_GRD! 功能：输出结果!! 输入参数说明：! A：存放输出数据的二维数组! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! m1,m2：实际数据x方向起点位置和终点位置! n1,n2：实际数据y方向起点位置和终点位置! 输出参数说明：! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! filename:存放输出结果的文件名! field_vdr:对场值作水平导数的结果! field_thdr:对场值作垂向导数的结果! field_asm:场值的解析信号振幅的结果!************************************************************************** SUBROUTINE OUTPUT_GRD(A,filename,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax) real eigvalinteger m1,m2,n1,n2,m,nreal A(1:m,1:n)character*(*) filenamereal amin,amaxamin=HUGE(amin)amax=-HUGE(amax)write(*,*)m1,m2,n1,n2,m,npauseeigval=3.701411*1e5DO j=n1,n2,1do i=m1,m2,1if(ABS(A(i,j)).lt.eigval) thenamin=MIN(amin,A(I,j))amax=MAX(amax,A(I,j))end ifend doEND DOOPEN(20,file=filename,status='unknown',form='formatted') write (20,'(A)')'DSAA'write (20,*)m2-m1+1,n2-n1+1write (20,*)xmin,xmaxwrite (20,*)ymin,ymaxwrite(20,*)amin,amaxDo j=n1,n2,1write(20,*) (A(i,j),i=m1,m2)End doClose(20)END subroutine OUTPUT_GRD。