实验五傅立叶变换与频率域滤波
一、实验目的
1.理解傅立叶变换；
2.熟悉MATLAB中各种傅立叶变换相关的函数；
3.掌握频域滤波的步骤以及MATLAB的实现方法；
4.理解频域滤波器与空域滤波器的关系。
二、实验内容及步骤
1、傅立叶变换及傅立叶反变换
（1）傅立叶变换相关函数
MATLAB提供了几个和傅立叶变换相关的函数。其说明如下：
F=fft2(f); 二维傅立叶变换
real(F); 傅立叶变换的实部
imag(F); 傅立叶变换的虚部
abs(F); 获得傅立叶频谱
fftshift(F); 将变换的原点移至频率矩形的中心
ifft2(F); 二维傅立叶反变换
iffshift(F); 反中心平移
（2）傅立叶频谱
傅立叶频谱反映了图像的频率成分。幅值谱的能量往往集中于中低频部分，并且中低频部分的能量反映了图像的实体。图像的噪声往往集中于高频部分。
下面的例子对课本中123页的图Fig4.03(a).jpg进行傅立叶变换，得到傅立叶频谱。
例：x=imread('Fig4.03(a).jpg');
F=fft2(x); %二维傅立叶变换
FP=sqrt(real(F).^2+imag(F).^2); %计算傅立叶频谱，或者使用abs()函数
imshow(uint8(FP)) %显示傅立叶频谱，直流成分分布在四个边角
figure(2); imshow(uint8(fftshift(FP))) % 中心平移的频谱图
思考题1：对课本125页的图Fig4.04(a).jpg进行傅立叶变换，得到傅立
叶频谱，为清楚地显示该谱，将其进行对数变换处理，增强其灰度细节。结果类似于图5_1。
图5_1 Fig4.04(a)的傅立叶谱
x=imread('Fig4.04(a).jpg');
F=fft2(x);
FP=sqrt(real(F).^2+imag(F).^2);
Image=log(1+double(fftshift(FP)));
imshow(x);
figure(2);imshow(Image,[]);
（3）傅立叶变换对
将一幅图像进行傅立叶变换，再进行傅立叶反变换，可以得到原始图像。下面的例子对课本中134页的图Fig4.11(a)进行傅立叶变换，然后再进行傅立叶反变换，观察并了解实现过程。
理论上，经过傅立叶变换和傅立叶反变换，应该得到原始图像。比较fa 和f ，基本差不多。查看变量C 的内容，可以看到他们之间还是有细微的差别的。 例： close all
f=imread('Fig4.11(a).jpg');
F=fft2(f); %傅立叶变换
ff=ifft2(F); %反傅立叶变换
fa=real(ff); % 取实部
subplot(2,2,1),imshow(f);
subplot(2,2,2),imshow(F); %只显示变换的实部
subplot(2,2,3),imshow(ff,[]);
subplot(2,2,4),imshow(fa,[]);
figure(2)
C=fa-double(f); %反变换回来的图像和原图相减
imshow(C,[]); %两者是有差别的
2、频率域滤波
按照频域滤波的步骤，在MATLAB 中是很容易编程实现频域滤波。由于滤波器就是频率域中的函数，关键是如何构造滤波器函数。频域滤波的步骤为： 对图像进行傅立叶变换DFT ， 即F （u,v ）；
傅立叶变换原点中心平移；
用滤波器函数H(u,v)乘以F （u,v ）；
将原点反中心平移；
傅立叶反变换；
取上一步结果中的实部，即频域滤波后的结果。
（1） 低通滤波
低通滤波是使低频部分通过，而使高频部分受到抑制，从而使图像变得平滑。常用的低通滤波有理想低通滤波、巴特沃斯低通滤波和高斯低通滤波。
下面的例子实现了对课本135页图Fig4.11(a)进行理想低通滤波处理，截止频率D 0＝30。
理想低通滤波函数：
00
1((,)0H u v ≤?=??，D u,v)D ，D(u,v)>D ，221/2(,)[(/2)(/2)]D u v u M v N =-+- 例：
close all
clear all
f=imread('Fig4.11(a).jpg');
% ------构造理想低通滤波器----------
[Hh Hw]=size(f); %以图像的行列值作为滤波器的行列
H(1: Hh,1: Hw)=0; %滤波器的初值为0
x0=Hh/2; y0=Hw/2; % 滤波器的中心点
for x=1:Hh
for y=1:Hw
if(sqrt((x- x0)*(x- x0)+(y-y0)*(y- y0))<30) %理想低通滤波器D0=30 H(x,y)=1;
end
end
end
% -----------傅立叶变换及频域滤波---------
F=fft2(f);
F=fftshift(F); %中心平移
FF=F.*H; % 频域滤波。注意必须是点乘
FFS=ifftshift(FF);
g=ifft2(FFS); % 傅立叶反变换
g=uint8(real(g)); %取变换后的实部
% 显示并比较结果
figure(1), imshow(f); title('原图')
figure(2),imshow(H);title('理想低通滤波器')
figure(3), imshow(log(abs(F)+1),[ ]);title('频率谱') %对数变换，便于显示结果figure(4), imshow(log(abs(FF)+1),[ ]);title('频域滤波')
figure(5),imshow(g,[]);title('滤波结果')
思考题2：参考上例，使用理想低通滤波器对图Fig4.11(a).jpg进行频域滤波。将理想低通滤波器的截止频率D0分别设为15、30、80，对结果进行分析与
比较。
clear all
f=imread('Fig4.11(a).jpg');
% ------构造理想低通滤波器----------
[Hh Hw]=size(f); %以图像的行列值作为滤波器的行列
H(1: Hh,1: Hw)=0; %滤波器的初值为0
x0=Hh/2; y0=Hw/2; % 滤波器的中心点
for x=1:Hh
for y=1:Hw
H(x,y)=1/(1+(sqrt((x- x0)*(x- x0)+(y-y0)*(y- y0))/15)^4); %巴特沃斯低通滤波器D0=15
end
end
end
% -----------傅立叶变换及频域滤波---------
F=fft2(f);
F=fftshift(F); %中心平移
FF=F.*H; % 频域滤波。注意必须是点乘
FFS=ifftshift(FF);
g=ifft2(FFS); % 傅立叶反变换
g=uint8(real(g)); %取变换后的实部
% 显示并比较结果
imshow(f); title('原图')
figure(2);imshow(H);title('理想低通滤波器')
figure(3);imshow(g,[]);title('滤波结果')
思考题3：构造巴特沃思低通滤波器，并使用它对图Fig4.11(a).jpg 进行低通滤波处理。其中阶数n=2，截止频率D0分别为15、30、80，和思考题2使用理想低通滤波器滤波结果进行比较和分析。
注：巴特沃思低通滤波器（BLPF ）的传递函数为：
201(,)1[(,)/]
n H u v D u v D =+，221/2(,)[(/2)(/2)]D u v u M v N =-+-
f=imread('Fig4.11(a).jpg');
% ------构造理想低通滤波器----------
[Hh Hw]=size(f); %以图像的行列值作为滤波器的行列
H(1: Hh,1: Hw)=0; %滤波器的初值为0
x0=Hh/2; y0=Hw/2; % 滤波器的中心点
for x=1:Hh
for y=1:Hw
H(x,y)=1/(1+(sqrt((x- x0)*(x- x0)+(y-y0)*(y- y0))/30)^4); %巴特沃斯低通滤波器 D0=30
end
end
end
% -----------傅立叶变换及频域滤波---------
F=fft2(f);
F=fftshift(F); %中心平移
FF=F.*H; % 频域滤波。注意必须是点乘
FFS=ifftshift(FF);
g=ifft2(FFS); % 傅立叶反变换
g=uint8(real(g)); %取变换后的实部
% 显示并比较结果
imshow(f); title('原图')
figure(2);imshow(H);title('巴特沃斯低通滤波器')
figure(3);imshow(g,[]);title('滤波结果')
（2）高通滤波
高通滤波是使高频部分通过，而使低频部分受到抑制，从而使图像边缘锐化。常用的高通滤波有理想高通滤波、巴特沃斯高通滤波和高斯高通滤波。
思考题4：参考理想低通滤波器的构造方法，构造理想高通滤波器，并使用它对图Fig4.11(a).jpg进行理想高通滤波处理。截止频率D0分别为15、30、80，对结果进行分析和比较。
三、频域滤波器的空间形式
由卷积定理我们知道，空间域滤波和频域滤波之间存在对应关系。频域滤波器在空间域存在对应的空间域滤波器。。对频域滤波器进行反傅里叶变换、中心平移，然后取其实部，就得到对应的空间滤波器。
例：close all
clear
%-------构造理想低通滤波器－－－－
H(1:256,1:256)=0;
x0= 256/2; y0= 256/2;
for x=1:256
for y=1:256
if(sqrt((x- x0)*(x- x0)+(y-y0)*(y- y0))<30) % 注1
H(x,y)=1;
end
end
end
figure
imshow(H,[])
fh=ifft2(H); %频域滤波器进行傅立叶反变换
fh=fftshift(fh);
figure
imshow(log(1+abs(real(fh))),[])
title('理想低通滤波器的空间形式')
思考题5：观察截止频率分别为15、30、80时，理想低通滤波器的空间形式。了解：频域滤波器越窄，在空域的滤波器越宽，即空域滤波模板越大，滤除的低频成分越多，图像越模糊；以及理想低通滤波器的截止频率越小振铃现象越明显。
四、实验报告及要求
1、按照上述实验步骤把实验结果截图，并附上文字说明。
2、回答思考题，写出实现的命令及实验结果截图。
3、使用word文档写报告，交电子文档；
波阻抗相位（FDEM） MT/AMT/CSAMT频率域电磁法勘探反演所用的波阻抗反演方法，测量点必须位于波区（又叫做平面波区或远区）同时测量相互正交的电场分量和磁场分量，电场与磁场的比值具有阻抗的量纲，称为波阻抗，用符号Z来标示，x方向的电场与y方向的磁场比值记为Z xy。 注意： Zxy：是复数 K：波数，是复数 ω：角频率 μ：磁化率 σ：电导率 ρ：电阻率 均匀介质中电场相位角落后于磁场，这个角度就是MT/AMT/CSAMT勘探数据处理过程中所给出的振幅和相位曲线中的相位曲线。 视电阻率计算公式如下：
当平面电磁场垂直入射均匀大地时，即使不知道场源强度，只要测量出大地表面相互正交的一对电场和磁场，便可以确定大地的电阻率，而选用不同的频率可达到不同的勘探深度，这就是天然场源MT/AMT 或人工场源CSAMT的波阻抗反演的理论基础。 大地电磁测深一般要测量相互正交的两个水平电场Ex，Ey和相互正交的两个水平磁场Hx，Hy（MT测量过程中还要测量垂直磁场Hz）。测量两个水平电场是用两对不极化电极，电极距一般为100～200米。因为AMT和MT的天然电磁场信号较弱，应该采取措施避免测量电线晃动切割地球磁场产生的噪声。测量磁场则是用两个相互正交的匝数很多的高导磁芯线圈。 MT/AMT/CSAMT波阻抗反演数据处理流程电磁场的测量是在时间域进行的，再用傅里叶变换将测量信号转换为频率域信号。测量电磁场信号的采样时间间隔应使截止频率高于所需的最高频率，采样时窗宽度应大于所需的最低频率对应的周期。为了避免数据量太大，当需要测量的频带范围较宽时，一般分为几个频段采样，并分段作傅里叶变换。测量电磁场的频率范围应使最高频率对应的穿透深度为所需探测的第一层厚度的几分之一，最低频率对应的穿透深度为最大勘探深度的数倍。为了去除局部电磁场的影响，现在实际测量中采用所谓的“远参考系统”，除测点外，还在距离测点数十公里以外的地方设立一个参考点，同时进行测量。测量数据中属于平面电磁场的信号应该是互相关的，而局部干扰电磁场的信号是互不
1 基本原理简介 1.1 MATLAB 简介 MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品 语言和交互式环境，主要包括MATLAB和Simulink两大部分。 MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程 大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并 己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。
1.2傅立叶变换基本原理 傅里叶变换能将满足一定条件的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。最初傅里叶分析是作为热过程的解析分析的工具被提出的[1]。傅里叶变换属于谐波分析。傅里叶变换的逆变换容易求出，而且形式与正变换非常类似。正弦基函数是微分运算的本征函数，从而使得线性微分方程的求解可以转化为常系数的代数方程的求解。在线性时不变的物理系统内，频率是个不变的性质，从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取。卷积定理指出：傅里叶变换可以化复杂的卷积运算为简单的乘积运算，从而提供了计算卷积的一种简单手段。离散形式的傅里叶变换可以利用数字计算机快速的实现（其算法称为快速傅里叶变换算法（FFT））。 MATLAB中提供的变换函数 （1） fft2:用于计算二维快速傅立叶变换，语句格式： B=fft2(I,m,n) 按指定的点数计算m,返回矩阵B的大小为m×n,不写默认为原图像大小。（2）fftn:用于计算n维快速傅立叶变换 （3）fftshift:用于将变换后的图像频谱中心从矩阵的原点移到矩阵地中心，语法格式：B=fftshift(I) (4)ifft2:用于计算图像的二维傅立叶反变换，语法格式 B=ifft2(i) (5)ifftn:用于计算n维傅立叶变换 快速卷积实验：傅立叶变换一个重要特性是可以实现快速卷积。 设A为M×N矩阵，B为P×Q的矩阵，快速卷积方法如下： *对A和B补0，使其大小都为(M+P-1)×(N+Q-1) *利用fft2对矩阵A和B进行二维变换 *将两个FFT结果相乘，利用ifft2对得到的乘积进行傅立叶反变换
第五章 频率域方法
第5章频域分析法 基本要求 5－1 频率特性 5－2 典型环节的频率特性 5－3 系统的开环频率特性 5－4 频率稳定判据 5－5 系统闭环频率特性与阶跃响应的关系5－6 开环频率特性与系统阶跃响应的关系
基本要求 1. 正确理解频率特性的概念。 2. 熟练掌握典型环节的频率特性，熟记其幅相特性 曲线及对数频率特性曲线。 3. 熟练掌握由系统开环传递函数绘制系统的开环对 数幅频渐近特性曲线及开环对数相频曲线的方法。 4. 熟练掌握由具有最小相位性质的系统开环对数幅 频特性曲线求开环传递函数的方法。
5. 熟练掌握乃奎斯特稳定判据和对数频率稳定判据及其 它们的应用。 6. 熟练掌握稳定裕度的概念及计算稳定裕度的方法。 7. 理解闭环频率特性的特征量与控制系统阶跃响应的定 性关系。 8. 理解开环对数频率特性与系统性能的关系及三频段的 概念，会用三频段的分析方法对两个系统进行分析与比较。
一、控制系统在正弦信号作用下的稳态输出 5－1 频率特性 ()sin r r t A t ω=输入信号： 2 2)(ω ω+=s A s R 其拉氏变换式
输出 1()n i i i C B D C s s s s j s j ωω==++ -+-∑1 ()() ()() i n s t j t j t i i t s c t C e De Be c t c t ωω-== ++=+∑拉氏反变换得[()]2()2 j j r j A e π φωφω∠-=?22 ()()()2r s j r A D s s j s A j j ω ω φωω φω==?-+=?其中
目录 1.题目........................................................ .................................. .2 2.要求........................................................ (2) 3.设计原理........................................................ . (2) 3.1 数字滤波器基本概念......................................................... (2) 3.2 数字滤波器工作原理......................................................... (2) 3.3 巴特沃斯滤波器设计原理 (2) 3.4脉冲响应不法......................................................... .. (4) 3.5实验所用MATLAB函数说
明 (5) 4.设计思路........................................................ .. (6) 5、实验内容........................................................ . (6) 5.1实验程序......................................................... . (6) 5.2实验结果分析......................................................... . (10) 6.心得体会........................................................ . (10) 7.参考文献........................................................ . (10) 一、题目：巴特沃斯数字低通滤波器 二、要求：利用脉冲响应不变法设计巴特沃斯数字低通滤波器，通带截止频率100HZ,采样频率1000HZ，通带最大衰减为0.5HZ，阻带最小衰减为10HZ，画出幅频、相频相应相应曲线。并假设一个信号x(t)=sin(2*pi*f1*t)+sin(2*pi*f2*t),其中f1=50HZ,f2=200HZ。用此信号验证
作业5 1、用理想低通滤波器在频率域实现低通滤波 程序代码如下： clear; A=imread('picture4.jpg'); I=rgb2gray(A); figure(1); imshow(I); title('原图像'); g = imnoise(I, 'gaussian' ,0 ,0.01); J = I+g; figure(2); imshow(J); title('加高斯噪声后图像'); s=fftshift(fft2(I)); figure(3); imshow(abs(s),[]); title('图像傅里叶变换所得频谱'); [a,b]=size(s); a0=round(a/2); b0=round(b/2); d=150; for i=1:a for j=1:b distance=sqrt((i-a0)^2+(j-b0)^2); if distance<=d h=1; else h=0; end; s(i,j)=h*s(i,j); end; end; s=uint8(real(ifft2(ifftshift(s)))); figure(4); imshow(s); title('低通滤波后所得图像');
得到的图像如下：
2、用理想高通滤波器在频率域实现高频增强 程序源代码如下： clrar; A=imread('picture5.jpg');
I=rgb2gray(A); figure(1); imshow(I); title('原图像'); s=fftshift(fft2(I)); figure(2); imshow(abs(s),[]); title('图像傅里叶变换所得频谱'); figure(3); imshow(log(abs(s)),[]); title('图像傅里叶变换取对数所得频谱'); [a,b]=size(s); a0=round(a/2); b0=round(b/2); d=150; p=0.2;q=0.5; for i=1:a for j=1:b distance=sqrt((i-a0)^2+(j-b0)^2); if distance<=d h=0; else h=1; end; s(i,j)=(p+q*h)*s(i,j); end; end; s=uint8(real(ifft2(ifftshift(s)))); figure(4); imshow(s); title('高通滤波所得图像'); figure(5); imshow(s+I); title('高通滤波所得高频增强图像'); 得到的图像如下：
《重磁资料处理与解释》实验二频率域位场处理和转换实验 学院：地测学院 专业名称：勘查技术与工程 学生姓名： 学生学号： 指导老师： 提交日期：2018年1月9日 二0一八年一月
目录 1 基本原理 (2) 1.1位场的方程 (2) 1.2二维傅里叶变换及卷积性质 (2) （1）傅里叶变换 (2) （2）卷积性质 (2) 1.3频率域位场延拓原理 (3) 2 输入/输出数据格式设计 (3) 2.1 输入数据格式设计 (3) 2.2 输出数据格式设计 (3) 2.3 参数文件数据格式设计 (3) 3 总体设计 (4) 3.1频率域位场处理与转换的一般步骤 (4) 3.2软件总体设计结果流程图 (4) 4 测试结果 (5) 4.1 测试参数 (5) （1）向上延拓 (5) （2）向下延拓 (5) 4.2 测试结果 (6) 5 结论及建议 (7) 附录：源程序代码 (8)
1 基本原理 1.1位场的方程 由场论知识可知，位场方程分为 两大类：有源的Possion 方程()02 ≠?U ，以及无源的Laplace 方程()02 =?U 。 Laplace 方程的第一边值问题()1|f U S =通常为Dirichlet 问题，第二边值问题 ?? ? ??=??2|f n U s 通常称为Nueman 问题。若P 点在S 平面内称为内部问题，反之称为外部问题。由唯一性定理可知，Dirichlet 的内部和外部问题的解是唯一的，而Nueman 内部问题的解不是唯一的，有一常数差，但其外部问题解是唯一的。 外部问题的解的唯一性的原因：。 0; 0=??=∞ →∞ →r r n U U 无源区域位场可以表示为： ds n G W n W G p W ??? ? ?????-??= π41)( （1-1） ()() ()()()[] ()() z y x h W d d z y x W z z y W -=-+-+--=??+∞∞-+∞ ∞ -ξξηεη εξηεξηεπξ,,*,,,,2,,x 2 3 22 2 （1-2） 1.2二维傅里叶变换及卷积性质 （1）傅里叶变换 []??+∞∞-+∞ ∞ -+-= =dxdy y x g y x g F v u G e vy ux i ) (2),(),(),(π （1-3） []? ?+∞∞-+∞ ∞ -+-= =dudv v u G v u G y x g e F vy ux i ) (21 ),(),(),(π （1-4） （2）卷积性质 ()()[]()()v u P v u G y x p y x F ,*,,*,g = （1-5） ()()[]()()y x p y x v u P v u G F ,*,g ,*,1=- （1-6）
燕山大学 课程设计说明书题目：等波纹低通滤波器的设计 学院（系）：里仁学院 年级专业：仪表10-2 学号： 学生姓名： 指导教师： 教师职称：
燕山大学课程设计（论文）任务书 院（系）：电气工程学院基层教学单位：自动化仪表系 2013年7月5日
摘要 等波纹最佳逼近法是一种优化设计法，它克服了窗函数设计法和频率采样法的缺点，使最大误差（即波纹的峰值）最小化，并在整个逼近频段上均匀分布。用等波纹最佳逼近法设计的FIR数字滤波器的幅频响应在通带和阻带都是等波纹的，而且可以分别控制通带和阻带波纹幅度。这就是等波纹的含义。最佳逼近是指在滤波器长度给定的条件下，使加权误差波纹幅度最小化。与窗函数设计法和频率采样法比较，由于这种设计法使滤波器的最大逼近误差均匀分布，所以设计的滤波器性能价格比最高。阶数相同时，这种设计法使滤波器的最大逼近误差最小，即通带最大衰减最小，阻带最小衰减最大；指标相同时，这种设计法使滤波器阶数最低。实现FIR数字滤波器的等波纹最佳逼近法的MATLAB信号处理工具函数为remez和remezord。Remez函数采用数值分析中的remez多重交换迭代算法求解等波纹最佳逼近问题，求的满足等波纹最佳逼近准则的FIR数字滤波器的单位脉冲响应h（n）。由于切比雪夫和雷米兹对解决该问题做出了贡献，所以又称之为切比雪夫逼近法和雷米兹逼近法。 关键词：FIR数字滤波器 MATLAB remez函数 remezord函数等波纹
目录 摘要---------------------------- ----------------------------------------------------------------2 关键字------------------------------------------------------------------------------------------2 第一章第一章数字滤波器的基本概-------------------------------------------------4 1.1滤波的涵义----------------------------------------------------------------------4 1.2数字滤波器的概述-------------------------------------------------------------4 1.3数字滤波器的实现方法-------------------------------------------------------4 1.4 .数字滤波器的可实现性------------------------------------------------------5 1.5数字滤波器的分类-------------------------------------------------------------5 1.6 FIR滤波器简介及其优点----------------------------------------------------5- 第二章等波纹最佳逼近法的原理-------------------------------------------------------5 2.1等波纹最佳逼近法概述-------------------------------------------------------9 2.2.等波纹最佳逼近法基本思想-------------------------------------------------9 2.3等波纹滤波器的技术指标及其描述参数介绍---------------------------10 2.3.1滤波器的描述参数-----------------------------------------------------10 2.3.2设计要求-----------------------------------------------------------------10 第三章matlab程序------------------------------------------------------------------------11 第四章该型滤波器较其他低通滤波器的优势及特点--------------------12 第五章课程设计总结---------------------------------------------------------------------15 参考文献资料-------------------------------------------------------------------------------15
设计题目频率域滤波的MATLAB设计与实现
目录 摘要...................................................................................................................... - 3 - 1. 数字图像处理. (1) 1.1发展概况： (1) 1.2关键技术： (1) 2．频率域滤波的产生背景及意义 (3) 2.1傅立叶级数和变换简史： (3) 2.2频率域滤波的意义： (3) 3. 频率域滤波的常用方法 (4) 3.1低通滤波 (4) 3.1.1理想低通滤波器的截面图 (5) 3.2高通滤波 (7) 3.3带阻滤波 (9) 3.4带通滤波 (10) 4．原理及实现 (10) 4.1频率域增强基本理论 (10) 4.2傅立叶变换 (11) 4.3频率域理想低通（ILPF）滤波器 (12) 4.3.1理想低通滤波器的截面图 (12) 4.3频率域巴特沃兹(Butterworth)低通滤波器 (13) 4.4频率域高斯(Gaussian)低通滤波器 (14) 5．程序设计 (15) 5.1算法设计(程序设计流程图) (15) 5.2 对灰度图像进行Fourier变换的程序 (15) 5.3频率域理想低通滤波器 (15) 5.4 二阶巴特沃斯(Butterworth)低通滤波程序 (16) 5.5 高斯(Gaussian)低通滤波程序 (17) 6．结果与分析 (19) 6.1 对灰度图像进行Fourier变换后的频谱图 (20) 6.2二阶巴特沃斯（Butterworth）低通滤波结果与分析 (20)
数字图像处理
本章包含的主要内容
傅立叶变换 卷积和卷积定理 频率域低通滤波 频率域高通滤波
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 )dxdy
6
1、低通滤波器（LPF） 低通滤波器是用来通过低频信号，衰减或抑制高频信号。 如图13－2（a）所示，为典型的二阶有源低通滤波器。它由两级RC滤波环节与同相比例运算电路组成，其中第一级电容C接至输出端，引入适量的正反馈，以改善幅频特性。 图13－2（a）二阶低通滤波器电路图 图13－2（b）二阶低通滤波器电路仿真图 电路性能参数： 二阶低通滤波器的通带增益
截止频率，它是二阶低通滤波器通带与阻带的界限频率。 品质因数，它的大小影响低通滤波器在截止频率处幅频特性的形状。 2、高通滤波器（HPF） 与低通滤波器相反，高通滤波器用来通过高频信号，衰减或抑制低频信号。 只要将图13－2低通滤波电路中起滤波作用的电阻、电容互换，即可变成二阶有源高通滤波器，如图13－3所示。高通滤波器性能与低通滤波器相反，其频率响应和低通滤波器是“镜象”关系，仿照LPH分析方法，不难求得HPF的幅频特性。 图13－3 二阶高通滤波器电路图 电路性能参数A uf、f0、Q各量的函义同二阶低通滤波器 3、带通滤波器（BPF）
图13-4 二阶带通滤波器 这种滤波器的作用是只允许在某一个通频带范围内的信号通过，而比通频带下限频率低和比上限频率高的信号均加以衰减或抑制。这种滤波器的作用是只允许在某一个通频带范围内的信号通过，而比通频带下限频率低和比上限频率高的信号均加以衰减或抑制。 典型的带通滤波器可以从二阶低通滤波器中将其中一级改成高通而成。如图13－4所示。 电路性能参数： 通带增益中心频率 通带宽度选择性 的比例就可改变频宽而不影响中心频率。 此电路的优点是改变R f和R 4 4、带阻滤波器（BEF） 如图13－5所示，这种电路的性能和带通滤波器相反，即在规定的频带内，信号不能通过（或受到很大衰减或抑制），而在其余频率范围，信号则能顺利通过。
截止频率 | | (2013-10-07 23:50:04) 转载▼ 分类：Vision 在物理学和电机工程学中，一个系统的输出信号的能量通常随输入信号的频率发生变 化（频率响应）。截止频率（英语：Cutoff frequency —）是指一个系统的输出信号能量开始大幅下降（在带阻滤波器中为大幅上升）的边界频率。 概述 电子滤波器等信号传输通道中的诸如低通、高通、带通、带阻等频带特性都应用了截止频率的概念。截止频率有时被定义为电子滤波器的导通频带和截止频带的交点，例如电路标称输出信号减3分贝的位置的频率。在带阻滤波器中，截止频率则被定义在输出信号能量大幅上升（或大幅下降）、失去“阻止”（或失去“通过”）信号效果的位置。在波导管或者天线的例子中，截止频率通常包括上限频率和下限频率。 截止频率的概念除了在电子工程有广泛应用，截止频率的概念还在等离子区振荡中 有所应用。 电子学 参见：波德图及分贝 在电子学中，截止频率是电路（例如导线、放大器、电子滤波器）输出信号功率超出或低于 传导频率时输出信号功率的频率。通常截止频率时输出功率为传导频率的一半，在波德图上
相当于为降低3分贝的位置所表示的功率，因为此时功率比例它:①将能m传到频带上的 Slk^fwdl 输出功率［2］。 低通滤波器的截止频率 右图所示为一个一阶的低通滤波器。它的截止频率由下式决定： 严—翻字 当信号频率低于这个截止频率时，信号得以通过；当信号频率高于这个截止频率时, 信号输出将被大幅衰减。这个截止频率即被定义为通带和阻带的界限。 高通滤波器的截止频率 右图所示为一个一阶的高通滤波器。它的截止频率由下式决定： JQ_2TT RC固 当信号频率高于这个截止频率时，信号得以通过；当信号频率低于这个截止频率时, 信号输出将被大幅衰减。这个截止频率即被定义为通带和阻带的界限。
2 化极原理 空间域位场转换复杂的褶积关系, 在频率域表现为简单的乘积形式. 即由实测异常的傅里叶变换频谱乘上相应的转换因子, 再反变换, 就是需要的转换结果, 其中转换因子可以是单个, 也可以是多种转换因子的组合, 这是频率域处理转换的突出优点例如实际资料的化极计算, 转换因子就应该包括去除高频干扰的滤波因子与化极因子的组合, 这类组合在频率域实现起来非常方便。化极计算涉及到磁化方向转换与测量方向转换, 该方向转换因子一般形式为H(u,v)=─—, ⑴ 其中q k=i(ul k+ vmk)+ nk√u2+ v2,(k=0,1,2,3),i=√-1,u,v为x,y方向的圆频率；lk=cosIk·cosDk，mk=cosIk·sinDk，nk=sinIk为方向余弦，Ik，Dk分别为磁化方向（和测量方向）的倾角和偏角；q0和q1分别为测量方向和磁化强度方向的频率域因子；q2和q3分别为转化后的测量方向和磁化强度方向因子。 当为化极时：I2=I3=90o,q2=q3= u2+ v2,且现在经常测量的是总场磁异常T , 其对应的测量方向是地磁场方向. 假设磁化强度方向与地磁场方向一致( 特别在稍大一点测区, 总是这样考虑) , 因此有q0= q1, 具体化极因子可简化为 用u= rcosθ, v= rsinθ代入( 2) 式, 即得极坐标系下的转换因子H( r,θ ) 为 其中r= u2+ v2, = arctan( u /v ) . 可以清楚看出频率域化极因子H( r, ) 是角度的单一函数, 与频率的高低无关, 因而可写成H(θ) . 上述频率域化极因子为扇形放大因子, 其数值直接依赖于磁倾角. 在I0= 0的极端情况下, 即磁赤道附近, 化极因子为 当θ= D0±90o时, H (θ) →-∞, 其特点见图1所示. 当磁倾角I0较小时, 化极因子的放射状线的极大值近似与磁倾角的平方成反比, 即 在接近该线较窄的扇形区域, 化极因子幅值升幅很快. 由( 5) 式可知, 在θ接近D0±90o时, H (θ) 数值很大, 造成计算结果很不稳定, 表现为化极结果沿磁偏角方向D0条带明显, 这是化极因子在θ= D0±90o方向由低频到高频的放大造成的. 为此, 需要对化极因子进行改造, 压制沿D0±90o方向的放大作用, 使计算稳定, 减少甚至消除条带现象. 然而, 化极因子沿D0±90方向的放大作用是其重要特征,改造得太严重, 就会失去其特征, 同样得不到理想的化极结果. 从理论上讲, 化极因子的所有特征都保留, 对应的必然是理论的化极结果. 但实际计算中,一方面数值必须是有限的, 超过则计算会溢出或误差很大. 另一方面, 数据是有限的、离散的, 其频谱必然与理论谱有误差, 该误差必然会被化极因子传递, 由于化极因子是放大因子, 沿某一方向一定范围内由低频到高频都放大, 计算中的误差就会放大传递, 对计算结果必然带来影响, 有时影响是巨大的.因此, 实际计算过程中应该对化极因子中不致于造成溢出的部分保留( 可逆部分) , 而对会引起数值极度放大的部分( 不可逆部分) 进行压制. 本文为此提出针对性措施压制因子法. 2. 1压制因子法 根据低纬度化极因子的平面、剖面特征( 图1) , θ0= D0±90o为死亡线 ,θ0 ,α0的扇形区域为死亡地带 ,α0为一较小的角度. 为了压制靠近D0±90o附近的过度放大效应, 设计一个压制因子, 该因子在D0±90o附近趋于零, 即压制作用最强; 一定范围以外等于1, 即不压制. 另外要求压制因子足够光滑. 余弦函数具有很好的特点, 对其加以改造, 可以满足上述要求. 为此设计如下压制因子: F( u, v) = F0 ( ) , 那么该滤波因子的特征( 图2) 应有
航磁数据位场转换处理及效果 航磁T 测量数据是不同深度、不同形态、规模的磁性地质体磁场信息在观测面上的综合反映。由于场的叠加效应，使得某些具有一定地质意义的异常变得复杂，在原始图件上很难识别，给地质解释工作带来了难度。为了提高对航磁异常的分辨能力，突出更多有用信息，根据测区航磁异常特征和地质解释需要，对原始测量数据进行了原平面化极、上延、垂向一阶导数以及剩余异常提取等几种位场转换处理。 第一节位场转换处理及效果 航磁平面网格数据位场转换处理采用表达式简单、运算速度快捷的频率域算法，进行化极、导数换算、解析延拓等处理。频率域转换的过程是：首先对异常资料进行傅立叶正变换，以得到异常资料的频谱；而后把异常的频谱和与转换相应的频率相应函数点积，得到处理后异常的频谱；最后对处理后异常的频谱进行傅立叶反变换，从而得到处理后的异常。 位场转换处理使用的软件是中国国土资源航空物探遥感中心自主开发的 WINDOWS系统下地球物理数据处理解释软件(GeoProbe Mage)及航空物探彩色矢量成图系统( AgsMGis)。 一、原平面化极处理 化极，即化磁极，就是把斜磁化异常转变为垂直磁化异常，相当于在磁北极观测异常。测区处于中纬度地区，由于倾斜磁化的影响，造成磁异常中心不是正好对应在地质体的正上方，而是相对于地质体的中心向南部产生一定的偏移。这对于确定磁性地质体的空间位置、形态、分布范围以及对磁异常的定性定量解释均带来一定的困难。化极可用于消除由于非垂直磁化引起的异常不对称性，在剩磁很小或感磁远大于剩磁且两者方向一致的情况下，将实测的斜磁化异常转化为垂直磁化异常，这样可以较为准确的确定异常的场源位置，提高异常解释的定位精度。从而使异常形态简化，并与磁性体位置保持一致，有利于圈定磁性体边界和走向。 作化极处理时要注意剩磁的影响，化极处理一般都假定磁化方向与地磁场方向一致，对于那些剩磁远远大于感磁且剩磁方向与地磁场方向不一致的磁性体就不符合这一假设条件，特别是测区中的火山岩分布区，由于剩磁较大会出现磁场畸变现象，使用时应注意甄别。从项目组野外物性测量结果看，区内多数岩石以感磁为主，剩磁方向与感磁方向接近，符合化极的前提条件。 全区采用"频率域偶层位变倾角磁方向转换方法"实现磁场全变倾角化极。在观 测面上建立笛卡尔直角坐标系，使x轴志向磁北，z轴垂直向下。假设观测场T是
航磁数据位场转换处理及效果 ?测量数据是不同深度、不同形态、规模的磁性地质体磁场信息在观测航磁T 面上的综合反映。由于场的叠加效应，使得某些具有一定地质意义的异常变得复杂，在原始图件上很难识别，给地质解释工作带来了难度。为了提高对航磁异常的分辨能力，突出更多有用信息，根据测区航磁异常特征和地质解释需要，对原始测量数据进行了原平面化极、上延、垂向一阶导数以及剩余异常提取等几种位场转换处理。 第一节位场转换处理及效果 航磁平面网格数据位场转换处理采用表达式简单、运算速度快捷的频率域算法，进行化极、导数换算、解析延拓等处理。频率域转换的过程是：首先对异常资料进行傅立叶正变换，以得到异常资料的频谱；而后把异常的频谱和与转换相应的频率相应函数点积，得到处理后异常的频谱；最后对处理后异常的频谱进行傅立叶反变换，从而得到处理后的异常。 位场转换处理使用的软件是中国国土资源航空物探遥感中心自主开发的WINDOWS系统下地球物理数据处理解释软件（GeoProbe Mager）及航空物探彩色矢量成图系统（AgsMGis）。 一、原平面化极处理 化极，即化磁极，就是把斜磁化异常转变为垂直磁化异常，相当于在磁北极观测异常。测区处于中纬度地区，由于倾斜磁化的影响，造成磁异常中心不是正好对应在地质体的正上方，而是相对于地质体的中心向南部产生一定的偏移。这对于确定磁性地质体的空间位置、形态、分布范围以及对磁异常的定性定量解释均带来一定的困难。化极可用于消除由于非垂直磁化引起的异常不对称性，在剩磁很小或感磁远大于剩磁且两者方向一致的情况下，将实测的斜磁化异常转化为垂直磁化异常，这样可以较为准确的确定异常的场源位置，提高异常解释的定位精度。从而使异常形态简化，并与磁性体位置保持一致，有利于圈定磁性体边界和走向。 作化极处理时要注意剩磁的影响，化极处理一般都假定磁化方向与地磁场方向一致，对于那些剩磁远远大于感磁且剩磁方向与地磁场方向不一致的磁性体就不符合这一假设条件，特别是测区中的火山岩分布区，由于剩磁较大会出现磁场畸变现象，使用时应注意甄别。从项目组野外物性测量结果看，区内多数岩石以感磁为主，剩磁方向与感磁方向接近，符合化极的前提条件。 全区采用"频率域偶层位变倾角磁方向转换方法"实现磁场全变倾角化极。在观 ?是测面上建立笛卡尔直角坐标系，使x轴志向磁北，z轴垂直向下。假设观测场T
实验四傅立叶变换与频率域滤波 实验目的 通过本次实验，实现以下几个目标： 1.理解傅立叶变换； 2.熟悉MATLAB中各种傅立叶变换相关的函数； 3.掌握频域滤波的步骤以及MATLAB的实现方法； 4.理解频域滤波器与空域滤波器的关系。 实验内容 一、傅立叶变换及傅立叶反变换 1.傅立叶变换相关函数 MATLAB提供了几个和傅里叶变换相关的函数。其说明如下： F=fft2(f); 二维傅立叶变换 abs(F); 获得傅立叶频谱 fftshift(F); 将变换的原点移至频率矩形的中心 ifft2(F); 二维傅立叶反变换 real(ifft2(F)); 提取变换后的实部 imag(ifft2(F)); 提取变换后的虚部 2.傅里叶频谱 傅里叶频谱反映了图像的频率成分。下面的例子对课本中123页和 125页的图Fig4.03(a) 和图Fig4.04(a)进行傅立叶变换，得到傅立叶 频谱。显示傅立叶频谱时，使用了对数变换以获得更好效果。 f=imread('Fig4.03(a).jpg'); F=fft2(double(f)); F=fftshift(F); figure(1), imshow(f); figure(2), imshow(log(abs(F)+1),[ ]); f=imread('Fig4.04(a).jpg'); F=fft2(double(f)); F=fftshift(F); figure(1), imshow(f); figure(2), imshow(log(abs(F)+1),[ ]); 为了更好地理解频谱，显示下面三个图像（x6.jpg，x60.jpg，y6.jpg）的傅里叶频谱，观察并比较、分析结果。 显示频谱时使用下面的语句来做灰度变换可找出其主要的频率成分。 figure(2), imshow(log(abs(F)+1).^4,[ ]); %先对数、再幂次变换
陈建友，王晓凯，杨长春.一种均匀网格反泄露傅里叶变换的频率域高效实现方法[J ].石油物探，2018，57（1）:45 -49 CHEN Jian y ou ，WANG Xiaokai ，YANG Chan g chun.Efficient im p lementation of re g ular g rid antileaka g e Fourier transform in the fre q uenc y domain [J ].Geo p h y sical Pros p ectin g for Petroleum ，2018，57（1）:45 -49收稿日期:2017-01-15;改回日期:2017-11-09. 作者简介:陈建友（1989 ），男，硕士，研究方向为地震信号处理和地震储层解释.Email :x j tuc j y 2009＠1６https://www.doczj.com/doc/356526056.html, 基金项目:国家自然科学基金（41504092，41774135）二中国博士后科学基金（201６T90925，2015M5725６7） 及中央高校基本科研业务费专项资金（x jj 201６0６5）联合资助.This research is financiall y su pp orted b y the National Natural Science Foundation of China （Grant Nos.41504092，41774135），China Postdoctoral Science Foundation （Grant Nos.201６T90925，2015M5725６7），and the Fundamental Research Funds for the Central Universities （Grant No.x jj 201６0６5）. 一种均匀网格反泄露傅里叶变换的频率域高效实现方法 陈建友1，2，王晓凯1，3，杨长春3 （1.西安交通大学，陕西西安710049;2.酒泉卫星发射中心，甘肃酒泉732750;3.中国科学院地质与地球物理研究所，北京100029） 摘要:反泄露傅里叶变换已在地震数据规则化二多维数据去噪及数据压缩等领域得到广泛应用，但其多在时问域实现，计算效率不高.对均匀网格反泄露傅里叶变换的常规实现方法进行了改进，在频率域给出一种均匀网格反泄露傅里叶变换的高效实现方法.该方法在频率域寻找傅里叶系数极大值及其位置，通过在频率域减去SINC 函数，避免了常规实现方法在每次迭代过程中所使用的离散傅里叶变换及逆离散傅里叶变换，大幅减少了均匀网格反泄露傅里叶变换的运算量，提高了均匀网格反泄露傅里叶变换的计算效率.利用合成数据对方法进行了验证，结果表明，该方法不仅能够明显减少频谱泄露现象，提高稀疏表示性能，而且较常规实现方法提高了计算效率. 关键词:反泄露傅里叶变换;均匀网格;SINC 函数; 稀疏表示;频谱泄露;数据压缩中图分类号:P６31文献标识码:A 文章编号:1000-1441（2018）01-0045-05DOI :10.39６9/j . issn.1000-1441.2018.01.00６Efficient im p lementation of re g ular g rid antileaka g e Fourier transform in the fre q uenc y domain CHEN Jian y ou 1，2，WANG Xiaokai 1，3，YANG Chan g chun 3 （1.Xi an J iaoton g Universit y ，Xi an 710049，China ;2.J iu q uan Satellite Launch Center ，J iu q uan 732750，China ;3.Institute o f Geolo g y and Geo p h y sics ，Chinese Academ y o f Sciences ，Bei j in g 100029，China ）Abstract :The antileaka g e Fourier transform （ALFT ）has been widel y used in seismic data re g ularization ，multi -dimensional de -noisin g ，and data com p ression.However ，the com p utational efficienc y of re g ular g rid ALFT is not hi g h enou g h in the time domain to p rocess massive seismic data ，because it im p lements one discrete Fourier transform to convert the si g nal from the time domain to the fre q uenc y domain in each iteration.In this p a p er ，we focus on an im p lementation of the ALFT in the fre q uenc y domain to im -p rove its efficienc y .The method searches for the maximum Fourier coefficient and subtracts the SINC function in the fre q uenc y do -main ，to avoid the discrete Fourier transform and the inverse discrete Fourier transform in each iteration.Tests on s y nthetic data show that the method can im p rove com p utational efficienc y ，while reducin g discernible fre q uenc y leaka g e and im p rovin g s p arse re p resentation.Ke y words :antileaka g e Fourier transform ，re g ular g rid ，SINC function ，s p arse re p resentation ，s p ectrum leaka g e ，data com p ression 近年来，压缩感知（com p ressive sensin g ，CS ） 技术[1]已在地震勘探领域引起极大关注并取得了诸多 成功应用[2-６].多维数据规则化二多维数据去噪及压缩是地震数据处理中的重要步骤[7-8]，其中最为重要的是寻找合适的稀疏变换，这同时也是压缩感知的重 要步骤.XU 等[9]于2005年提出了反泄露傅里叶变 万方数据