短时傅里叶变换

短时傅里叶变换（精选十篇）

短时傅里叶变换 篇1

关键词：Wigner-Ville分布,短时傅里叶变换,反褶积,时窗长度,地震主频

元坝气田是中国首个埋藏最深的生物礁气田,位于四川盆地川中低缓构造北斜坡与九龙山和通南巴构造带倾末端的过渡带,二、三叠系海相地层埋藏深,构造形变弱,断裂不发育。大面积发育的长兴组礁滩储层为元坝气田的主力储层,主要发育台地边缘礁滩复合体和台地边缘礁后浅滩两种类型储层,具有“早滩、晚礁,前礁、后滩”特征,受沉积微相控制,储层发育具有复杂的空间展布形态、较强的非均质性以及气藏的复杂气水分布。储层埋深大于6 500 m,单层厚度薄,尤其是台地边缘礁后浅滩储层,厚度相对礁滩复合体普遍薄,由滩心到滩缘厚度减薄、物性变差。受地震分辨率的影响,元坝气田超深层储层( 平均井7 048 m深) 特别是薄储层的精细描述和气水识别困难,深度挖掘地震信息,发展薄储层的油气检测技术,能提高元坝气田的勘探开发效益,对四川盆地及我国南方碳酸盐岩油气勘探具有较好的借鉴意义。

低频伴影分析是进行流体识别的主要方法之一。Ebrom[1]的研究指出了产生含气层地震低频伴影的10 个可能的影响因素。Castagna等[2]利用地震记录的时-频分析,将低频伴影识别碳氢化合物的应用效果提高到一个新水平。贺振华[6]根据低频伴影检测油气储层的效果,总结和概述了应用低频伴影特征作为流体识别标志的准则。陈学华等[10]将广义S变换用于单频剖面的提取,取得了一定的效果。但杨小江等[7]指出,薄储层的低频伴影现象不明显,因此针对礁后滩薄储层,想要有效地进行低频伴影分析,需要发展高分辨的时频分析技术。朱恒等[8]、Lu W K等[3]提出的反褶积短时傅里叶变换用于地震信号的时频分析,在提高低频伴影分析分辨率方面取得了一定的进展。但值得注意的一个问题是,时窗的选择对时频分辨率有较大的影响。通过建立模型以及实际地震数据实验,得到了最佳时间和频率分辨率的地震主频之间的关系。元坝地区的应用实例表明,通过该关系可获得较高时频分辨率,有利于薄储层的低频伴影分析。

1 反褶积短时傅里叶变换原理

1. 1 短时傅里叶变换

信号x( t) ∈ L2( R) 的短时傅里叶变换定义为:

式( 1) 中,ωt,f( τ) = ω( τ - t) e- i2πfτ是窗函数,常用的是高斯窗,* 是共轭转置。短时傅里叶变换谱定义为

信号短时傅里叶变换的时窗选取对分析结果有较大影响,时窗增大可以提高频率分辨率,但是会降低时间分辨率; 同理时窗变窄可以提高时间分辨率,但是会降低频率分辨率。因此,短时傅里叶变换中时间分辨率和频率分辨率存在相互受限。

1. 2 Wigner-Ville分布

对于一个信号x( t) ,它的Wigner-Ville分布定义为

Wigner-Ville分布具有较高的时频分辨率,但是其本身存在的交叉项对信号分析存在较大的干扰影响,它是由Wigner-Ville分布非线性引起的。

1. 3 反褶积短时傅里叶变换

信号x( t) 的短时傅里叶变换谱的表达式为:

信号x( t) 的短时傅里叶变换谱可以写成以下褶积的形式:

式( 5) 中WVDx和WV Dh分别为信号x( t) 和窗函数h( t) 的Wigner-Ville分布。短时傅里叶变换谱的交叉项在大部分情况为零值。式( 5) 也可写成

式( 6) 中,Sx是信号x( t) 的短时傅里叶变换谱,Wx是信号x( t) 的Wigner-Ville分布,Wh是窗函数h( t)的Wigner-Ville分布,**代表二维褶积。如果已知Wh和Sx就可以求出Wx。从短时傅里叶变换谱获取Wigner-Ville分布估计值是一个反褶积问题。

处理该反褶积问题一种效果较好的算法是Lucy-Richardson反褶积算法,该算法是一种非线性的迭代复原算法。从最大似然估计的角度出发,假设地震信号服从Poisson分布,通过k次迭代可以得到Wigner-Ville分布。

Lucy-Richardson算法表达式为

式( 7) 中,n + 1 是现在的迭代次数,Wx-( 0) = Sx。其中* 为相关算子,为褶积算子。反褶积短时傅里叶变换运算流程图如图1。

反褶积结果Wx-与Wigner-Ville分布有相近的时频分辨率,而且由于对短时傅里叶变换谱进行反褶积而减少了交叉项。因此,反褶积短时傅里叶变换综合了两种方法的优点,具有较高的时间频率分辨率。

实际运用中发现: 同一主频的地震剖面在反褶积短时傅里叶变换中采用不同长度时窗,其时频分析结果存在较明显差异。对含油层做时频分析中,采用不同的时窗长度,其“低频阴影”现象存在差异。因此,在使用该方法进行油气检测,很有必要对时窗选取进行探讨。

2 时窗选取对合成信号的影响

2. 1 合成交叉信号分析

设分析合成信号为: 采样点nsample: 256,采样时间间隔: dt = 1 /nsample。其信号表达式为

该信号加入10 d B的高斯白噪声如图2。对该信号采用不同长度的时间窗口时频分析结果如图3。

根据图3 分析可知: 采用不同长度时窗对信号做反褶积短时傅里叶变换的结果存在较为明显的差异,在上述图中表现为: 30 个采样点高频端的时间分辨率较好,但是交叉处低频端频率分辨率相对较差。40 个采样点高频端的时间分辨率较好,且交叉处低频端频率分辨率也较好。50 个采样点频率分辨率总体较好,但高频端的时间分辨率变差。60 个采样点频率分辨率与50 个采样点相近,但是高频端的时间分辨率更差。

根据以上分析可以确定时窗大小对反褶积短时傅里叶变换的时频分辨率存在影响。随着时窗长度的增加,频率分辨率提高,时间分辨率降低; 反之则频率分辨率降低,时间分辨率提高。需要选取合适时间窗口使得频率分辨率和时间时间分辨率达到相对最优化。

在油气低频阴影检测中,较高的时间和频率分辨率对检测结果至关重要,尤其是在对薄层分析中,频率分辨率和时间分辨率要求更高。为此,对两个反射层的地震记录加以分析以说明时窗选取的重要性。

2. 2 两层地震记录分析

设两层反射系数模型,中间地层含油。模型1:初始入射地震子波主频为80 Hz雷克子波,根据黏弹性衰减理论,设入射至下部地层的雷克子波主频衰减为60 Hz。其反射系数和地震记录1 如图4。

选取不同时窗长度的时频分析结果如图5。

从图5 中可看出,不同时窗长度对地震信号的时频定位效果存在明显差异。其规律是: 随着选择时间窗口长度增加,频率分辨率提高,时间分辨率明显降低。图5( a) 和图5( b) 两个时窗长度对应分析结果可以在时频域将两反射层较好区分开。对比图5( a) 和图5( b) 两结果可以发现两者的时间分辨率虽然大致相同( 即准确将上下反射层定位至30 和40 单位时间点) ,但是图5 ( a) 中10 个采样点分析结果的频带拉伸较严重,频率分辨率较差,这种频带拉伸影响对低频伴影的“高频上强下弱”识别是很不利的。图5( c) 和图5( d) 两个时窗长度对应分析结果的频率分辨率较高,但是时间分辨率很差。上下两层反射波信息无法在时间轴上区分开,且位于30 与40 单位时间上下反射层时间内也包含了两层的反射信息,造成虚假时频信息。这种层间干扰对薄层的识别影响较严重。综合对比分析表明: 15 个采样点时频分辨率较好,上下两层主频及时间定位准确且频带无明显拉伸。

为验证15 个采样点时窗是否适用于其他主频地震数据,改变上下两层入射子波主频,建立模型2: 顶层子波主频50 Hz,底层子波主频衰减为30Hz。其反射系数和地震记录2 如图6。

选取不同时窗长度的时频分析结果如图7。

从图7 中可看出,15 个采样点已不适用于该主频的地震记录。15 个采样点结果存在严重的高频带拉伸现象,且上层反射位于30单位时间点的主频错误地定位在60 Hz附近。对比表明:33个采样点时频分析结果较好,上下两层即(30与40单位时间点)与子波主频50 Hz和30 Hz很好地对应。综合以上分析,对于不同主频的地震数据需要选取对应合理的时窗长度,以达到低频伴影分析的要求:既需要较高的时间分辨率(即不希望目的层之外的信息影响目的层的“低频伴影”分析),又需要较高的频率分辨率(即不希望用于“低频伴影”分析的剖面频带过宽),因此应特别注意时窗的选择。

3 地震主频与时窗长度

通过模型分析表明: 时间窗长度选取中应当考虑到地震主频。随着地震主频升高,时间窗长度应该适当减少。笔者对按照第2 部分的方法,同时结合大量实际地震数据分析比较,获得了主频与时窗组合如表1。

采用高斯函数拟合该统计关系

曲线拟合如图8。

其函数表达式为

4 实际应用

实际分析的工区为元坝地区过Yb104 井的地震剖面,如图9。实钻证实该红色区域含气。利用反褶积短时傅里叶变换提取单频剖面进行“低频伴影”分析。

根据地震频谱分布可知,该地震剖面主频为20Hz附近。主频时窗关系式表明其适宜的时间窗长度为50 个采样点。为此,分别采用时间窗长度为25、50、75 对该剖面进行反褶积短时傅里叶变换对比分析,其分析结果单频剖面依次为图10、图11、图12。

以上分析结果表明: 针对同一主频的地震数据,选用不同时窗长度的单频剖面结果有存在较大差异。时窗为25 样点的单频剖面结果具有较高的时间分辨率,但是20 Hz与40 Hz两剖面对比表明没有明显衰减变化,其原因是相对于20 Hz主频的地震剖面采用25 个采样点长度过短,存在高频带拉伸从而严重影响其频率分辨率。时窗为50 个采样的反褶积短时傅里叶变换“低频阴影”效果明显,同时其时间分辨率也相对较好,能较好刻画层间分布信息。时窗为75个采样点的反褶积短时傅里叶变换没有“低频阴影”效果,且时间分辨率较差,时窗过长导致层间信息相互干扰是导致该结果的主要因素。

5 结论

MAtlab傅里叶变换实验报告 篇2

信工 142

学号

姓名

何岩

实验组别

实验日期

室温

报告日期

成绩

报告内容:(目得与要求, 原理, 步骤, 数据, 计算, 小结等)1、求信号得离散时间傅立叶变换并分析其周期性与对称性;给定正弦信号 x(t)=2＊ｃｏs(2＊pｉ*1０＊ｔ),ｆs=100ＨZ,求其ＤTFＴ。

(ａ)代码: ｆ=10;Ｔ=1／f;w=－1０:0、2:10;t1=0:0、000１:１;ｔ２＝0:0。0１:1;n1=-2;ｎ2=8;ｎ０=０;n＝n1:0。0１:n２;ｘ5＝［ｎ>=0.01];x１=2＊cos(2*ｆ＊pi＊ｔ1);x２＝2*coｓ(2＊f*pi＊t２);x3=(exp(—ｊ).^(t2’*w));x４=ｘ２*x３;subplot(2,2,1);ploｔ(t1,x1);axｉs(［０ １ 1、1*min(x2)1。1＊maｘ(ｘ2)]);xlaｂeｌ(’x(n)’);ylabeｌ(’x(n)“);titｌｅ(＇原信号 x1”);ｘlａbeｌ(“t”);ylabeｌ(“x1’);sｕbpｌot(2,2,3);sｔem(t2,ｘ2);axis(［０ 1 1、1＊min(x2)1。１＊max(x2)]);title(’原信号采样结果 x2＇);xlabel(＇ｔ’);ｙlaｂel(＇x2”);ｓubplot(2,2,２);steｍ(n,x5);axiｓ([0 1 1、1*miｎ(x5)1.1*mａx(x5)]);xｌabｅｌ(’n’);ｙｌabel(＇x2“);titlｅ(’采样函数ｘ２＇);ｓubｐlｏt(2,2,4);ｓtｅm(t2,ｘ4);axis(［0 1 —0、2＋1。1*min(x４)１、1*ｍａx(x4)］);xlabel(’t”);ylabel(＇x4“);title(”ＤＴFＴ结果 x4＇);(b)结果:

２、用以下两个有限长序列来验证 DTFT 得线性、卷积与共轭特性;x1(n)=［１ ２ 3 ４ 5 6 7 8 9 10 11 １２];x2(ｎ)=R 10(n)(1)线性:(a)代码: w=linｓｐaｃe(-8,8,10000);ｎｘ1＝[０:11］;nx2＝［0:9];ｘ1=［1 2 ３ 4 ５ ６ ７ 8 9 10 1１ １２］;x2=[1 １ 1 1 １ 1 1 1 １ １];x3=［x2,ｚeｒｏs(１,(leｎgth(x1)—lｅngth(x2)))］;x4=2＊x1+3＊x3;Ｘ1=x1*exｐ(-j＊ｎx1＇＊w);％频率特性 X3=x３*exp(-j*ｎｘ1＇＊ｗ);%频率特性 X4＝ｘ4*exp(—ｊ*ｎｘ1’*w);％频率特性

suｂplｏt(5,３,１),ｓｔem(nx1,x１),ａxiｓ([-1,１3,0,15］);tiｔle(＇x１’), ｙlabel(“ｘ(n)’);

suｂｐｌot(5,3,2),stem(nx２,x２),ａxis(［—１,１3,0,5］);ｔitｌe(”x2＇);

subplot(5,３,3),steｍ(nx1,x4),axis([-1,13,0,2６]);titlｅ(’x4=2＊ｘ1+3＊ｘ３“);

sｕbplｏt(５,3,4),ｐｌｏｔ(w,abs(X1));ylabｅl(＇幅度’)

subplot(5,３,7),ｐlot(w,aｎgle(Ｘ1));ｙｌabel(’相位＇)

subpｌｏt(5,３,10),ploｔ(w,rｅal(Ｘ１));yｌabel(’实部’)

sｕbｐlｏt(５,3,13),pｌot(w,imag(X１));yｌａbｅl(”虚部’)ｓuｂploｔ(5,3,5),pｌot(ｗ,aｂs(X３));

ｓubpｌoｔ(5,3,8),plot(w,angｌe(Ｘ3));

subｐlot(５,3,11),ploｔ(w,reａl(X3));subplｏt(5,3,１4),plot(w,imag(X3));

subplｏｔ(５,3,6),ｐlot(ｗ,abs(X4));

sｕbplot(5,3,9),plot(w,anｇｌe(Ｘ４));

subpｌot(5,3,１2),ploｔ(w,rｅal(X４));sｕbploｔ(5,３,15),plot(w,ｉmag(Ｘ4));

(b)结果:

(2)卷积:(a)代码: nx１=0:11;ｎｘ2＝０:9;ｎx3=０:20;

w=ｌｉnspace(-８,８,40);％w=[—8,8］分 10000 份

x1＝[1 2 3 ４ ５ 6 ７ ８ 9 10 11 12］;x２＝［1 1 1 1 1 1 1 1 1 1];x3=cｏnv(x１,ｘ２);％ x1 卷积 x2 ｘ4=x1＊exp(-ｊ＊nx1“*w);％ x１频率特性 x5=ｘ2*exp(－ｊ＊ｎｘ2’＊w);% x2 频率特性 x6=x3*exp(-j＊nｘ3”＊w);% x1 卷积 x２频率特性 x7＝x４、*x5;

sｕｂpｌｏt(2,2,１),stｅm(ｎx1,x1),axｉｓ(［—1,15,0,15］),ｔｉtｌe(’x1“);su ｂ plo ｔ(2,2,2), ｓ ｔ em(nx2, ｘ ２),ax ｉ s([—1, １

5,0,5］),title(’x2’);ｓuｂplot(2,1,２),stem(ｎx3,ｘ3),axiｓ([—1,25,0,80]);ｔitle(＇x１卷积ｘ2 结果 x3’);figuｒe,subｐlｏt(2,２,1),stｅｍ(x４,”filled’),titｌe(“ｘ１得DTFT 结果ｘ4’);

subploｔ(2,2,2),ｓtem(x５,”filled＇),titlｅ(’x２得 DTFT结果 x5’);

ｓｕｂplot(2,２,3),stｅｍ(ｘ6,＇filｌｅd’),ｔitｌｅ(’ｘ3得 DTFT 结果 x6’);

subｐlot(2,2,4),ｓteｍ(x7,“fｉlled＇),ｔitle(＇x4 得ＤTFT 结果ｘ7’);

figure,suｂｐlot(3,２,1),stem(ｗ,abｓ(ｘ6)), ylabel(”幅度’),ｔｉｔle(’x1 卷积 x2 得 DTFＴ＇);

ｓuｂplｏｔ(4,2,3),steｍ(w,aｎｇle(ｘ6)),yｌａbeｌ(“相位”)

sｕbpｌoｔ(4,2,5),stem(w,ｒeａl(ｘ6)),ylabｅl(“实部’)

subplｏｔ(4,2,7),sｔem(w,ｉmag(x6)),ylabeｌ(＇虚部’)

ｓubｐlot(４,2,2),stem(w,abs(x７)), tｉtlｅ(’x1 与 x2 得 DTFＴ得乘积’);

sｕbplot(４,2,4),ｓtem(w,ａngle(x７));

subplot(4,2,6),stｅｍ(w,real(ｘ7));

subpｌot(4,2,８),sｔｅｍ(w,ｉｍag(x7));

(ｂ)结果:

(3)共轭:(a)代码: x1n=［1 2 3 4 5 6 7 ８ 9 1０ １1 12］;w=—10:10;N１=ｌｅngth(ｘ1n);n１=0:N1—1;x1=real(x1n);x2=iｍaｇ(ｘ１n);x2n=ｘ1—j＊x2;

Ｘ1＝x2n*(exｐ(－j)、^(n1＇＊ｗ));X2=x１n*(exp(j)、＾(n1’＊w));x3=rｅal(X2);x４＝imａg(Ｘ2);Ｘ２=x３—j*x4;fｉguｒe,subpｌｏt(211);stem(w,Ｘ1,”.’);tiｔｌｅ(“x1ｎ共轭得 DＴFT’);

suｂplot(2１2);sｔem(ｗ,Ｘ2,”、’);tｉｔle(“ｘ1n 得 DＴFT 取共轭且反折”);(b)结果:

3。

求 LTI 系统得频率响应 给定系统 H(Ｚ)=B(Ｚ)/A(Z),A＝[0。９8777 －0。3１18３ 0、０256] B=［0.98997 0.989 ０。989９7］,求系统得幅频响应与相频响应、(要求使用ｆｉlｔer(B,Ａ,δ(n))求解。

(a)结果: A=[0、98777-０。3１1８3 0、０２５6］;Ｂ=[0。９899７ 0、９89 0、９899７];C=[1 0 0 ０ 0 0 0 0 0 0 0 0 0 0 ０ ０ 0 0 ０ 0 ０ ０ 0 0 0 0 0］ y=filtｅr(B,A,C);suｂpｌot(2,2,1);stem(ｙ,＇、’);titｌｅ(’原始序列“);

mag=abs(ｙ);ｐｈ=ａnｇle(ｙ);ph=ph*1８0／ｐi;sｕbplot(２,2,2);ｓｔｅm(mａｇ,”、＇);title(＇幅频特性＇);xlａbel(＇时间信号ｎ“);ylabel(＇信号幅度＇);suｂplｏt(2,2,３);stｅm(ph,”、’);ｔitlｅ(“相频特性”);xlaｂｅｌ(“时间信号 n＇);yｌabｅl(”信号相位“);(b)结果:

4.采样与频谱混叠 给定信号ｘ(t)=100*eｘp(-100*ｔ)＊coｓ(2*ｐi＊500＊t),求该信号得频谱;当采样频率分别为 fs１=200０ＨZ,fｓ2=1000HＺ;fs3=５00HZ;fｓ４=200HZ,时输出序列得 DＴFＴ。

(a)代码: x=100＊ｅxp(-100*t)、＊cos(2＊pｉ＊500＊t);ｔ=—2:０、1:2;w=－10:0。1:10;

y＝ｘ*(exｐ(-j)、^(ｔ’*ｗ));subpｌot(２,1,1),plot(t,ｘ);ｓubpｌot(2,1,2),ｐlｏt(w,ｙ);tｉtle(’原始信号得频谱＇);figｕre,ｆｓ1＝2000;Ts１＝1／ｆs1;n1=-2:Tｓ1:2;

fs2＝10００;Ｔs2=1/fs２;n2=－2:Ts2:2;

fs3=500;Ts3=1/fｓ３;n3＝-2:Ts3:2;

fs4=200;Ts4=1/ｆｓ４;n4＝—2:Tｓ4:2;x1＝１00。＊eｘp(—10０＊n1)。*cos(２*pi*500*ｎ１);ｙ１=ｘ1＊(eｘp(-j)。^(n1”＊w));sｕbploｔ(221);plｏt(w,ｙ１);title(“经 2000Hz 采样后信号得 DTFT”);x2=100。*exｐ(-１0０*ｎ2)、*ｃos(２*ｐｉ*500＊n２);ｙ2＝x2＊(eｘｐ(-j)、^(n2＇*w));subplｏｔ(222);ｐｌot(w,y２);tｉｔle(’经 1000Hｚ采样后信号得 DTＦT’);x3=１00、＊ｅxp(—１00*n3)、＊coｓ(２＊pi＊5０0＊n3);

ｙ3＝x3*(exｐ(—j)、＾(ｎ３“＊w));suｂpｌｏt(22３);ｐlｏt(ｗ,y3);titlｅ(’经５00Hz 采样后信号得 DTFＴ”);x4=100.＊exp(—100*n4)。＊cｏs(２*pi*5００*ｎ4);y4=x４*(exp(—ｊ)、^(n4’＊ｗ));suｂplｏt(２24);ploｔ(w,y4);tiｔle(’经 2００Hｚ采样后信号得 DTＦT＇);(b)结果:

收获及感想: DFＴ针对得就是有限长数字信号得傅立叶变换或傅立叶时频分析问题。但 以前得傅立叶变换就是定义在整个时间轴上得,而且一般针对得就是连续信号 ,获得得就是一个连续得频谱。

离散傅里叶变换(DFT),就是傅里叶变换在时域与频域上都呈现离散得形式,将时域信号得采样变换为在离散时间傅里叶变换(DTＦT)频域得采样。在形式上,变换两端(时域与频域上)得序列就是有限长得,而实际上这两组序

列都应当被认为就是离散周期信号得主值序列。即使对有限长得离散信号作DＦT,也应当将其瞧作经过周期延拓成为周期信号再作变换。在实际应用中通常采用快速傅里叶变换以高效计算 DＦT。

关于从傅里叶变换到拉普拉斯的变换 篇3

【关键词】傅里叶变换 拉普拉斯变换 复频域

一、引言

利用傅里叶变换只能求系统函数的零状态响应，而不能求系统的零输入响应。在需要求零输入响应时，还得利用其它方法，例如经典的方法。由于傅里叶变换在工程上受到一些限制，所以现今在研究线性系统问题时引入了拉普拉斯变换。

二、傅里叶变换和拉普拉斯变换的意义

傅里叶变换简单通俗的理解就是把看似杂乱无章的信号考虑成由一定振幅、相位、频率的基本正弦（余弦）信号组合而成，傅里叶变换的目的就是找出这些基本正弦（余弦）信号中振幅较大（能量较高）信号对应的频率，从而找出杂乱无章的信号中主要振动频率的特点。而拉普拉斯变换是为简化计算而建立的实变量函数和复变量函数间的一种函数变换。对一个实变量函数作拉普拉斯变换，并在复数域中作各种运算，再将运算结果作拉普拉斯反变换来求得实数域中的相应结果，这种计算往往比直接在实数域中求出同样的结果简单的多。

三、从傅里叶变换到拉普拉斯变换

当函数 不满足绝对可积条件时，可采取给 乘以因子 （ 为任意实常数）的办法，这样即得到一个新的时间函数 。若能根据函数 的具体性质，恰当的选取 的值，从而使当 时，函数 ，既满足条件

则函数 即满足绝对可积的条件了，因而它的傅里叶变换一定存在。可见因子 起着使函数 收敛的作用，故称 为收敛因子。

设函数 满足狄里赫利条件且绝对可积（这可通过恰当的选取σ的值来达到），根据 可以得到 ，在上式中， 是以 的形式出現的。令 ，s为一复数变量，称为复频率。 的单位为 ， 的单位为 。这样上式变为 由于上式中的积分变量为t，故积分结果必为复变量s的函数，故应将 改写成 ，即

复变函数 称为时间函数 的单边拉普拉斯变换。 称为 的像函数， 称为 的原函数。一般记为

符号 为一算子，表示对括号内的时间函数 进行拉普拉斯变换。利用 （t>0）或 可推导出 反变换的公式，即 对上式等式两边同时乘以 ，并考虑到 不是 的函数而可置于积分号内。于是得 由于上式中被积函数是 ，而积分变量却是实变量 ，所以欲进行积分，必须进行变量代换。因 故 （ 因为任意实常数）故 且当 时， ；当 时， 。将以上这些式子带入 中即得 （t>0）

或写成 称为拉普拉斯变换，可以已知的像函数 求与之对应的原函数 。一般记为 符号 为一算子，表示对括号内的像函数 进行拉普拉斯变换。式子 与 （t>0）或 构成拉普拉斯变换对，一般记为 或

若 不是因果信号，则拉普拉斯变换式 的积分下限应改写为（ ），即 称为双边拉普拉斯变换。因为一般常用信号均为因果信号（即有始信号），所以我们一般主要讨论和应用单边拉普拉斯变换。由上述可知，傅里叶变换是建立了信号的时域与频域之间的关系，即 而拉普拉斯变换则是建立了信号的时域与复频域之间的关系，即

复频率平面是以复频率 的实部 和虚部 为相互垂直的坐标轴而构成的平面，称为复频率平面，简称s平面，如下图5-1所示。

复频率平面（即s平面）上有三个区域： 轴以左的区域为左半开平面； 轴以右的区域为右半开平面； 轴本身也是一个区域，它是左半开平面与右半开平面的分界轴。故将s平面划分为这样的三个区域，为以后研究提供了很大的方便。

四、结论

在线性微分方程的已知输入求输出时，若系统稳定，且输入信号的傅里叶变换存在，在零初始条件下（t<0，x（t），y（t）及各阶导数为零），两种方法求解结果相同，傅里叶变换同样可用于对瞬时过程的求解。

【参考文献】

[1]严晋强，机械工程测试技术，北京，机械工业出版社，1990

[2]关正毅，信号处理及信号变换，北京，清华大学出版社，1989

快速傅里叶变换的c++实现 篇4

傅里叶变换是一种谱分析的方法, 在数学与工程技术分析中有着广泛的应用。本文从傅里叶变换的原理介绍开始, 然后介绍适合计算机上运算的离散傅里叶变换即DFT (discrete Fourier transform) , 而由于普通离散傅里叶变换在计算机上进行多点运算时, 运算量过大, 人们开始从算法中进行研究, 发明了效率更高的计算傅里叶变换的方法, 即FFT (Fast Fourier transform) , 为了使读者更好的理解FFt, 本文给出了一个基2的N点FFt程序。

1 离散傅里叶变换的定义

我们的计算机只能处理离散的数据, 在机械工程上数据采集卡采集来的数据也都是离散的, 要对这些数据进行分析就要用到离散的傅里叶变换, 离散的傅里叶变换的定义如下:

对长度为N的复数序列A0, A1, AN-1称

为序列{Ak}的离散傅里叶变换DFT (discrete Fourier transform) 。离散傅里叶变换有时也称为有限傅里叶变换。这里, WN=exp (2πi/N) 。

2 快速傅里叶变换

显然按由{Ak}按插值的方法求{xj}需要N2次复数乘法运算。由于一般情况下人们可以主动选择N使之满足一定的条件, 以此为基础建立的快速傅里叶变换 (Fast Fourier transform) FFt算法可以大大减少复数乘法的计算量。比如当取N=2r时, N2=22r=4r, 建立的FFt算法的复量运算量为O (Nlog2N) =r2r。当N很大时, 运算量的节省是显著的。

FFT算法有效地利用WNk=exp (2πik/N) 的周期性。它具有运算量少, 稳定性好和精确度高等优点。由于WNj N=1, j为整数WNk+l=WNkWNl设N可表示为N=r×s, r, s为整数 (2.1)

简记j= (j1, j0) , 其中j=j1×r+j0j1=0, 1, …s-1;j0=0, 1, …r-1 (2.2)

简记k= (k0, k1) , 其中k=k1×s+k0

这时

记xj=x (j1, j0) , Ak= (k1, k0) , 则

我们知道直接计算{xj}需要N2个复数运算, 若分两步计算, 在 (2.6) 中k0是固定的, 可将WN- (j0k0+j0k1) 看成一个复数完成 (2.6) 共需要r2s=Nr次复数运算。从序列A1 (j0, k0) 计算序列x (j1, j0) , 即完成 (2.7) , 共需要r2s=Ns次复数运算。故由{Ak}求{xj}共需N (r+s) 次复数运算。如果将N分解成N=r1r2…rm逐次重复上述过程可以看出共需要复数运算为N (r1+r2+…+rm) 若考虑N=rm, ri=r, i=1, 2, …m, 则复数运算总量为

特别当r=2时, 则复数运算总量为2Nlog2N当N充分大时, N2和2Nlog2N相比差别是很大的。比如当N=216=65536

即该算法的运算量只有N2的2048分之一, 因此计算量的节约是巨大的。

3 基2的FFT算法c++程序实现

基2的FFT的算法的讲解在计算方法的书上有详细的讲解, 这里不再累述。以下给出c++的完整程序, 此程序在vc6.0中可以直接应用。

结论

本文主要介绍了实现FFT的c++算法的, 本程序可以直接在vc6.0上运行现在代入八个点进行测试, 把1, 1+i, 2+i, 3+2i, 1+2i, 0, 2, -1+i代入, 得到9+7i, -2.1213+0.5355i, 2, 0.7071+5.9497i, 1+i, 2.1213-0.5355i, -4, -0.7071-3.9497i。结果正确。

摘要：傅里叶变换是一种谱分析的方法, 在数学与工程技术分析中有着广泛的应用。为了使读者更好的理解FFt, 文章给出了一个基2的N点FFt完整程序, 该程序经测试可以直接在vc6.0中运行。

关键词：傅里叶变换,FF,c++

参考文献

[1]蒋长锦, 蒋勇.快速傅里叶变换及c程序[M].合肥:中国科技大学出版社, 2004.

读《活着》有感：生命的傅里叶变换 篇5

——读《活着》有感

最近读了余华的《活着》。这是一部读完让人心情沉重的作品。小说中的主人公福贵是民国时期的一个地主家的少爷，年轻时由于嗜赌放荡，输尽家财。父亲被气死后，福贵一家成为佃农，并很快被国军抓壮丁卷入国共内战。随着内战、三反五反，大跃进，“”等社会变革，他的人生和家庭也不断经受着苦难，所有亲人都先后离他而去，仅剩下年老的他和一头老牛相依为命。富贵的一生犹如坐滑梯，只不过是从云端滑到地面。

而在现在相对稳定、公平的社会环境下，大多数人是在爬滑梯，虽然最终只有少数人爬到云端，可只要努力，今天终归是比昨天站的高些。

看了一篇写傅里叶变换的文章：任何周期函数都可以看成不同振幅、不同相位正弦波的叠加，而正弦波就是一个圆周运动在直线上的投影。所以每个生命都可以看成不同领域(事业、家庭等)、不同关系(同事、亲朋等)的正弦波的叠加，而正弦波就是每个人在某一方面的生命投影(事业的波折、情感的经历等)。简单来说：想谱写好人生的周期函数就要在有限的生命长度里尽量拓展生命的宽度。

我在特别年幼时就爱看译制片，每到周末，一堆男女老少挤在学校一间十几平的小房子里，看一台小小的电视，每次我都能坚持到最后。说“特别年幼”是因为那时真的很小，小到不能理解那些外国电影的情节，但我记住了那些车水马龙、那些穿着连衣裙、踩着高跟鞋、拿着公事包的摩登女郎，心里好生羡慕，希寄有一天能成为其中一员。当我真的成为Office Lady又特别羡慕那些有事业、有家庭、有财富的Office Queen。看过一篇文章说：当平凡的你在一个场合看到气场强大的女王出现而心生羡慕时，还真的没有什么好办法一蹴而就，生活是努力、是磨炼，更是沉淀与积累。

就像见识了腾格里沙漠后就再也不会回味乌丹的玉龙沙漠一样，就像登过3500米贺兰山就不会再恐高1500米的泰山一样，就像跑完了10公里就再也不会抵触1000米一样，就像读了老舍就不再迷恋鲁迅一样，就像试了奶酪金砖就不在理睬提拉米苏一样，就像吃过切糕王子就不会再买三只松鼠一样，见识到更好的才不会做自大的夜郎。做一个懂生活、会生活的人，或者读书、或者健身、或者美食，尽可能尝试谱写不同的“正弦波”，并使其在你生命中留下不同的印记。不会的东西现在就学起来，毕竟每一天都是余生里的最年轻的一天。

短时傅里叶变换 篇6

关键词：傅里叶变换光谱仪； 光谱反演； 拟合算法； 可见光谱

中图分类号： TP 911.73 文献标志码： A doi： 10.3969/j.issn.10055630.2015.06.007

Abstract：In the Fourier transform spectrometer （FTS）， obtaining the optical path difference of each sampling points accurately is the key of data processing. The velocity fluctuation of moving mirror will cause equaltime sampling error， and traditional equal optical path difference sampling method can not be used for visible wavelengths. The influence of speed fluctuation is analyzed and a method to calculate the optical path difference of each sampling points based on fitting algorithm is developed. The method is verified by processing the interference signal of visible light with large velocity fluctuation. The experimental results verify its correctness and effectiveness. This method can also be used in any kind of Fourier transform spectrometer.

Keywords： FTS； spectrum recovery； fitting algorithm； vis spectrum

引 言

傅里叶变换光谱仪（FTS）因其具有高光谱分辨率、高光通量[12]、杂散光影响小、波数精确度高等优点在化学分析、大气探测等领域应用越来越广泛。其中可以获得非常高的光谱分辨率是傅里叶变换光谱仪最重要的优点，而光谱分辨率虽然受限于仪器自身，但是干涉信号采样误差和反演算法带来的误差同样会影响最终的光谱分辨率[34]。在傅里叶变换光谱仪中，随着动镜的移动，两干涉光束产生光程差，从而产生一个随光程差变化的干涉信号。干涉信号和光谱数据是一对傅里叶变换对[4]，对获取到的干涉信号以光程差为变量做傅里叶变换即可得到入射光的光谱数据，因此要进行光谱反演就需要在采样干涉信号的同时能精确地获取每个干涉信号采样点处的光程差。

如果动镜运动速度为绝对匀速，那么对干涉信号进行等时间采样等同于等光程差采样。但在实际情况中，动镜速度不可能是绝对匀速的，如果继续采用等时间采样则不可避免会导致非常大的采样误差，从而在得到的光谱图中产生误差[5]。常用的方法是采用激光作为参考信号，对其进行滤波整形后作为干涉信号采样的触发信号，这种方法只适用于目标光源波长大于参考激光两倍的傅里叶变换光谱仪中，例如中长波红外傅里叶变换光谱仪，当目标光源在近红外和可见光波段时，继续采用上述方法则会因不满足奈奎斯特采样定理而无法从离散的干涉信号中复原出光谱数据。本文基于曲线拟合算法提出了一种精确获取每个干涉信号等时间采样点处光程差的方法，仍然采用激光作为参考信号，对于波长在可见光波段甚至紫外波段的目标光源，即使速度波动较大，只要等时间采样频率满足采样定理，本方法都可以获取每个数据点处的光程差，从而对干涉信号进行精确反演获取光谱数据。

1 傅里叶变换光谱仪干涉信号分析

傅里叶变换光谱仪属于调频的干涉光谱仪，主要用来观测光谱信息，常用的傅里叶变换光谱仪采用的是经典迈克尔逊结构，光程差是动镜行程的两倍，其光学结构如图1所示。动镜采用角镜，这样可以消除动镜倾斜造成的影响[6]，而且在相同的动镜行程下可以获取四倍的光程差，从而得到更高的光谱分辨率。入射光被分束器分成强度相等的两束光，经分束器反射的光束经过左侧的定镜1和定镜2到达角镜左侧，然后反射回来，经分束器反射和投射分别到达探测器2和探测器1；而经分束器透射的光束经右侧的定镜3和定镜4到达角镜右侧，然后反射回来，经分束器反射和投射到达探测器1和探测器2。当动镜沿着光轴来回做直线运动时，在探测器上汇合的两光束的光程差发生周期性变化，形成干涉信号。在傅里叶变换光谱仪中，动镜作匀速扫描的过程，相当于在整个采样过程中以速度u对干涉光信号进行调制，其结果是把频率很高的光波调制成频率很低的电信号[7]。干涉信号被红外探测器转换为电信号，经过放大器和模拟滤波器，进入模数转换器进行等光程差采样，得到数字干涉图。

式（6）中后两项就是所谓的鬼线。如果采用等时间采样采集上述干涉信号，再直接对其做傅里叶变换，在所得到的光谱数据中会有额外的尖峰出现。如果入射光源为连续光谱，则所得到的光谱图会有很多鬼线出现。常用的方法是采用稳频激光作为参考信号，这是因为激光单色性好而且波长λ0已知，动镜每移动λ0/4便产生一个干涉信号周期，即每当激光干涉信号出现过零点时就知道光程差变化了λ0/2，因此只要以参考激光干涉信号过零点作为采样触发信号即可实现等光程差采样[9]。然而从采样点中复原出原谱线必须保证采样频率满足奈奎斯特采样定理，如果目标光源的波长λ全都大于激光波长的两倍，那么当光程差每变化一个λ，激光干涉信号已经至少出现了4个过零点，此时满足奈奎斯特采样定理；但是当目标光源的波长小于激光波长的两倍时，光程差每变化一个λ，激光干涉信号只出现一个过零点，甚至不出现。即对于工作波段为0.4～1.0 μm的傅里叶变换光谱仪，继续采用上述方法进行等光程采样则会造成干涉信息丢失，无法得到完整干涉信号，更无法得到光谱信息，因此如何对可见光干涉信号进行等光程差采样是能否反演出可见光光谱的关键。

2 利用曲线拟合算法获取光程差

在现在智能的控制方案下，虽然动镜往复运动的实际速度是不断变化的，但其速度不会存在突变的现象，因此在傅里叶变换光谱仪中，不可能存在动镜移动0.10～0.25 μm（λ/4）的过程中速度发生很大变化的情况，一般在动镜移动这么短的距离过程中，速度基本可以看作是匀速的。那么在光程差变化一个激光波长的过程中，产生的激光干涉信号应该是一个标准的带有初始相位的单频率正弦信号。

对于采集到的离散信号，很少能直接求得其函数表达式，一般是采用插值和拟合的方法，利用采集到的离散数据点来得出一条近似正确的连续的曲线。如果采样得到的离散信号点没有误差则采用插值方法，如果采样得到的离散信号点与真实值有差距则采用曲线拟合方法。常用的获取曲线的方法有样条插值、多项式插值、多项式拟合、基于遗传算法的拟合算法等，其中多项式拟合是最小二乘拟合的一种常用形式。最小二乘法是应用最广泛的曲线拟合算法，其核心思想是寻找合适的函数参数使得函数与所采样得到的数值之间的误差平方和达到最小[10]，即使下式达到最小值：

式中：N为采样总数。当用标准正弦信号去拟合时，即f（x）=Asin（2πfx+φ）+B，其中幅值A、频率f、相位φ、直流分量B四个参数都是未知的[11]，而且误差平方和是这四个参数的非线性函数，因此无法求出误差平方和的闭合解，只能用迭代法求出每个参数的局部最优解。

本文基于最小二乘法用f（x）=a+bsin（cx）+dcos（cx）对每个激光干涉信号周期进行正弦拟合，求出每个周期对应正弦信号的四个参数和拟合误差，如果能用该方法拟合出结果而且拟合误差非常小，则可近似认为：在光程差变化一个激光波长的过程中，产生的激光干涉信号是一个标准的单频率正弦信号。在傅里叶变换光谱仪中，由于参考激光干涉信号和目标光源干涉信号所经过电路的延迟不同，即使是采用等光程差采样仍然会有采样误差产生，根据文献[1213]的计算，只要速度波动的相对误差小于2%，这种误差对光谱的影响基本可以忽略。因此只要每个干涉信号周期的拟合误差小于2%，则基本可以忽略。对一个完整干涉信号的拟合流程图如图2所示。

3 实验验证

为验证本文所提出的方法，对所用傅里叶变换光谱仪进行速度开环控制，在速度有较大波动的情况下进行干涉信号采样。动镜采用无刷直流电机驱动，电机本身速度较快，经过减速箱减速，动镜运动的平均速度控制在3 cm/s，根据干涉信号频率和速度、波数的关系，可得干涉信号的频率范围为100～400 kHz，因此选择采样频率为2 MS/s，采样位数为16位。

参考激光采用波长为0.685 2 μm的稳频氦氖激光器，因对可见光波段气体的吸收率较低，为方便进行结果对比，目标光源分别采用波长为0.659 8 μm的红色激光和波长为0.532 μm的绿色激光。对动镜移动一个单程所获得的数据进行处理可以得到实际速度波动情况，如图3（a）所示，计算其波动的相对误差值为10.04%。拟合误差如图3（b）所示，最大拟合误差为1.30%，拟合误差的平均值为0.37%，满足误差小于2%的要求。图4和图5分别为用本文提出的算法对波长为0.659 8 μm和波长为0.532 μm的光源进行光谱反演的结果，并给出峰值光谱处的细节，从实验结果中可以看出，当速度变化较大时，采用本文所用的方法所得到的反演结果波数准确度高，而且基本可以达到仪器自身的光谱分辨率0.05 cm-1。

4 结 论

傅里叶变换光谱仪的动镜速度不可避免地会有波动，对等时间采样的干涉信号直接进行傅里叶变换将会产生很大的光谱误差。当目标光源为中长波红外光时，可以用参考激光干涉信号触发采样以实现等光程差采样，但是当目标光源的波长较短时，无法直接进行等光程差采样。本文提出了一种基于曲线拟合的光程差获取算法，能在速度波动较大时仍然能准确地反演出光谱图，而且适用于任何波段的目标光源，同时对速度均匀性无严格要求，降低了对控制系统精度的要求。

参考文献：

[1]夏翔.干涉光谱仪多速度模式扫描技术研究[D].上海：中国科学院研究生院（上海技术物理研究所），2014：1115.

[2]吴航行，华建文，王模昌.新型红外空间遥感用傅里叶变换光谱仪[J].红外与激光工程，2004，33（4）：397400.

[3]盛灏.傅里叶变换光谱仪干涉信号数据获取研究[D].上海：中国科学院研究生院（上海技术物理研究所），2014：1421.

[4]范世福.光谱技术和仪器的新发展[J].光学仪器，2000，22（4）：3540.

[5]王明艳.基于DSP的弹光调制非均匀干涉信号获取及反演算法研究[D].太原：中北大学，2013：3135.

[6]王文桂.干涉光谱仪[M].北京：宇航出版社，1988：180190.

[7]杨隆梓，段星辉，魏焕东，等.星载傅立叶光谱仪数据采集技术研究[J].科学技术与工程，2008，8（8）：20342037.

[8]胡盛雯，殷德奎.高光谱大气探测仪干涉信号的获取电路设计[J].红外与激光工程，2006，35（S1）：528533.

[9]胡玲，潘征宇，洪扁，等.双拼激光干涉仪中的数字相位计设计[J].光学仪器，2014，36（3）：258262.

[10]齐国清，吕健.正弦曲线拟合若干问题探讨[J].计算机工程与设计，2008，29（14）：36773680.

[11]梁志国，朱济杰，孟晓风.四参数正弦曲线拟合的一种收敛算法[J].仪器与仪表学报，2006，27（11）：15131519.

[12]WILLIAMS C S，Mirror misalignment in fourier spectroscopy using a Michelson interferometer with circular aperture[J].Applied Optics，1966，5（6）：10841085.

[13]COHEN D L，Performance degradation of a Michelson interferometer when its misalignment angle is a rapidly varying，random time series[J].Applied Optics，1997，36（18）：40344042.

[14]姚涛，吕群波，相里斌，等.非均匀采样干涉数据光谱反演技术研究[J].光谱学与光谱分析，2010，30（5）：14301433.

[15]吕群波，姚涛，相里斌，等.干涉数据光谱反演方法研究[J].光谱学与光谱分析，2010，30（1）：114117.

常见傅里叶变换的滤波性能分析 篇7

在电网发生故障时,交流信号中含有谐波分量和衰减直流分量,因此不能通过以正弦函数模型的算法来计算,如半周积分法、采样和导数法、两点(三点)乘积算法等[1]。而在一般输电线路(指中、低压网络或不长的高压输电线)上发生短路时,可以不用考虑非整数倍的高频分量;否则需要考虑采用基于随机模型的最小二乘法或卡尔曼算法。由于卡尔曼滤波算法在数据窗暂态条件下能给出基波分量的最优估计,但计算过于复杂,限制了实际应用[2]。这种情况下就可以通过采样周期函数模型的算法,即假设输入信号可分解为基波的整倍数频率分量,如傅里叶(以下简称傅氏)变换算法或沃尔什(Walsh)函数算法,来计算需要分解的基波分量。沃尔什函数尤其在相位计算时,计算量较大、耗时较多。相比而言傅氏变换算法具有计算量小、精度高、相位求取方便等优点,因此得到了广泛的应用,同时也针对其缺陷做了相应的改进。

1 傅氏变换简介

傅氏变换是数字信号处理的一种重要工具,来源于傅氏级数[3]。采用正交函数作为样品函数,将此正交样品函数与待处理的周期函数进行相应的积分变化,以求出与样品函数频率相同分量的实部和虚部,从而可以求出需分析的周期函数中该频率的分量模值和相位[4]。

1.1 基本傅氏变换

傅氏算法的基础是假定输入信号是周期函数,可以分解为整倍数频率分量之和,其中包括恒定的直流分量[5],即

式中m为谐波的次数,m=0,1,…,n;ω为基波角频率;Xcos、Xsin为第m次谐波的余弦和正弦分量的幅值。

根据傅氏变换的原理,可求出基波的正弦分量Xsin和余弦分量Xcos:

全周傅氏算法的滤波系数为可事先得到的常数,故算法的实时计算量不大[6]。

继电保护采样模件对连续的输入信号i(k)等间隔采样,把采样值依次转换成数字序列,采用离散傅氏变换DFT(Discrete Fourier Transform)。为了提高运算速度,一般都采用快速傅氏变换FFT(Fast Fourier Transform)。通过矩形积分可得全周傅氏算法:

式中N为一个周期的采样次数;k表示从故障开始时的采样点序号。

但这种算法需要一个周期内的N个采样数据,其数据窗为一个整周期,即20 ms[7]。实际上,无论采用何种算法或数字滤波器,要提高滤波性能,都不可避免地需要延长它们的数据窗[6]。

全周傅氏算法虽不能完全滤除输入信号中非整数次的高次谐波,但具有很大的抑制能力,滤波效果也是较好的。但其抑制低频分量(非周期分量引起)的能力较差[8]。

为了后一次计算可以利用前一次的计算结果,使用如下方式[9]:

该算法由于将与故障分量无直接联系的故障前数据作为故障后数据处理,算法的收敛速度较慢,精度较差[10]。但由于它主要用于微机保护等实时性要求较高的场合,因此在满足准确性要求的基础之上,计算量因素更显重要[11]。

有时为了提高响应速度,需要缩短数据窗采用半周傅氏算法:

但也会带来很多不利的影响。

1.2 改进傅氏变换

当故障信号中存在衰减直流分量时,半周傅氏算法的误差非常大,根本不能应用。全周傅氏算法的误差也很大,也不能直接使用[12]。为了减少这种误差,就提出了许多改进的傅氏变换。下面就介绍2种比较常用的方法。

1.2.1 差分傅氏变换

采用差分滤波的主要目的就是尽量减少衰减直流分量带来的误差。为了更好地抑制衰减直流分量的不利影响,需要选取尽量短的数据窗。

式中M≥1,为确定的常数,称为差分步长(当M=1时,抑制衰减直流分量的效果最好)。

计算量因每点均要计算差值而增加许多,且增加了算法对高频分量的敏感度[13]。

1.2.2 任意衰减时间常数的补偿

假设输入的信号为

式中I0表示衰减的直流分量;τ为时间常数;Im表示m次谐波的幅值;ω表示基波角频率;φm表示m次谐波的相角。

如果衰减直流分量的时间常数τ可以确定,就可以通过预先计算出补偿量来抑制误差,并能减少部分运算量。但是,对于实际的故障信号,时间常数是不能预知的,难以确定。所以需要在一个周期的数据窗的基础上增加一个采样点,就能够在事先未知衰减时间常数的情况下对衰减直流分量进行精确补偿,理论上属于一种精确算法[14]。需要变换得到:

经傅氏变换即得正弦分量和余弦分量补偿系数Ksin、Kcos:

从而得到经过补偿后基波的正弦分量和余弦分量:

2 仿真分析

假设输入含有衰减直流分量和5次谐波的交流信号,如图1所示,图中i(k)为输入交流信号波形,w(k)为实际基波信号波形;取N=24,τ=30 s,图2~4同。

其中基波信号

由于有衰减直流分量和高次谐波的叠加,输入信号相比基波波形已经发生明显的畸变。

2.1 全周和半周傅氏算法的比较

通过全周和半周傅氏算法滤波仿真计算,可以得到变换后的波形,如图2所示。图中Q(k)、B(k)分别为通过全周和半周傅氏算法滤波得到的基波波形,w(k)同图1。

显然,全周傅氏算法可以完全滤除高次谐波分量,仅受衰减直流分量的影响。而半周傅氏算法不能滤除偶次谐波分量,故与实际的基波相比还有较大的误差。

2.2 差分算法的比较

为了抑制衰减直流分量的影响,引入差分算法得到变换后的波形,如图3所示。图中,Q′(k)、B′(k)分别为通过差分全周和半周傅氏算法滤波得到的基波波形,w(k)同图1。

在差分步长取1时,可以很有效地抑制衰减直流分量的影响。

2.3 衰减直流分量补偿方法

对衰减直流分量进行精确补偿:一个周期的数据窗上增加一个采样点,在事先未知衰减时间常数的情况下对衰减直流分量补偿,补偿后的结果如图4所示。图中,Q″(k)、B″(k)分别为通过全周和半周补偿傅氏算法滤波得到的基波波形,w(k)同图1。

对衰减直流分量的补偿效果也很显著,缺点就是相对增加一些运算量。虽然相位发生了一点偏差,但是只要统一采用相同的算法,这种偏差是可以完全抵消的。

根据上述仿真后的图形可以定性地发现改进后的算法要远比改进前的性能提高了很多,再定量比较上述6种算法滤波性能及周期内不同采样次数的影响,如表1所示。

%

因此采样频率越高,算法的滤波特性周期越大,算法的整体滤波特性越好[15]。随着N的增加运算量也随之增加,但滤波的性能并不是呈线性地提高。

3 结语

综上所述,衰减直流分量补偿方法的滤波精度是最好的,但是其计算量也是最大的。相比而言差分算法计算量增加不大,而且精度提高了许多,值得推荐。同时,发现在采样次数N=24时可以达到速度和精度的最优效果。

对比采用全周和半周数据窗时,半周不能消除偶次谐波分量,同时也不能滤除直流分量,产生很大的误差;全周不仅可消除各次谐波分量和恒定的直流分量,而且还可有效地抑制衰减直流分量的影响。

为了兼顾保护的动作速度和测量精度,可以在故障一开始用半周傅里叶算法为灵敏度较低而安全裕度较大的保护提供数据,一周后改用全周算法相应地提高保护的灵敏度[3]。

摘要：为选取一种满足继电保护可靠性和灵敏度要求的滤波算法,对傅里叶变换算法的2种常见改进算法——差分算法和补偿算法进行分析。研究当输入模拟信号叠加谐波和衰减直流分量时,在不同的采样频率下,这2种改进算法的性能和效率。衰减直流分量实时补偿的算法精度较高,但计算量过大。差分算法能在计算量增加不大的情况下,大幅度地提高滤波精度。通过MathCAD仿真计算,分析上述变换的误差,并对比采用全周和半周不同数据窗造成的差异。在一周期24点采样时,即可达到精度和速度的最优效果。考虑兼顾保护的动作速度和测量精度,可以在故障初期采用半周数据窗的傅里叶算法为保护提供数据,后期改用全周数据窗算法。

基于分数傅里叶变换的数字水印技术 篇8

目前的数字水印算法大多为空域和频域水印两种类型, 其抵抗攻击的能力尚不是很令人满意, 而且在诸多性能指标之间难以平衡。介于空域和频域之间有效的数字水印算法大多集中于小波分解的方法, 利用小波的空频局域化特征表达能力, 将水印信息嵌入到小波域中。小波的双域分析能力使得该算法在数字水印的实际应用中起着重要的作用。

本文介绍另一种基于双域分析的数字水印算法, 即分数傅里叶变换谱域的水印算法。分数傅里叶变换谱也具有空域和频域双域信息表达能力, 但与小波变换不同的是分数傅里叶变换是全域的。该算法根据分数傅里叶变换谱的双域信息表达能力, 对水印载体图像进行某个级次的分数傅立叶变换。将加密水印嵌入到载体图像的分数傅立叶谱中, 然后再进行逆变换得到含水印的图像。

(一) 分数傅里叶变换

傅立叶变换 (Fourier Transform, 简称FT) 是一种重要的信号分析工具, 它将空 (时) 域信息变换为频域信息, 便于后续的图像处理与分析。分数Fourier变换是对经典Fourier变换的改进, 同小波变换把研究对象变换成维数更高的新对象来处理 (其实质是“升维”) 一样, 分数傅立叶变换FRFT (Fractional Fourier Transformation) 是经典Fourier变换在分数级次上的推广, 是一种介于函数与其傅里叶变换之间的信号双域描述。

其积分核

其中

a是FRFT的幂次, 可取任何实数。

二维分数傅立叶变换可以通过先作x轴 (或沿y轴) 的一维分数傅立叶变换, 再作沿y轴 (或沿x轴) 的一维分数傅立叶变换来实现。

当FRFT的幂次a从0连续增长达到1时, FRFT的结果相应地从原始信号的纯时间 (空间) 形式开始逐渐变换为它的纯频域 (谱) 形式, 幂次a在0到1之间的任何时刻对应的FRFT是包含时 (空) 域信息和频 (谱) 域信息的混合信号。显然, FRFT是一种时 (空) 频描述和分析的工具。

(二) 数字水印嵌入算法

水印信号的嵌入过程如图1所示, 输入包括水印和载体图像。水印可以是序列号、文字或图像等。对水印加密主要是用于提高水印的安全性。输出是含水印的图像。

对载体图像进行分数傅立叶变换, 然后将加密后的水印嵌入到变换域载体图像幅度谱中, 再将所得数据结合傅立叶变换相位谱作分数傅立叶逆变换, 得到含水印的图像。

1. 对载体图像进行二维分数傅立叶变换

图像的二维离散分数傅立叶变换 (DFRFT) 可通过对图像进行两次一维分数傅立叶变换来实现, 即先对每行进行级次为px的一维分数傅立叶变换, 再对变换后图像的每列进行级次为py的一维分数傅立叶变换。

对于256×256大小的256灰度级载体图像baboon图2进行px=0.4, py=0.1的256×256二维离散分数傅立叶变换 (DFRFT) , 得到其分数傅立叶变换幅度谱和相位谱。对DFRFT幅度谱矩阵图2 (b) 分别进行列和行降序排列, 因此, 幅值较大的点集中在256×256矩阵的左上角, 如图2 (c) 所示。变换域的载体图像幅度谱矩阵与排序后的矩阵之间存在一定的一一映射关系。

2. 水印加密

Arnold变换, 又称猫脸变换, 可看作剪切和拼接的过程, 通过这一过程将离散化的数字图像矩阵中的点重新排列。对一幅大小为N*N的图像有如下描述:

水印的加密算法采取Arnold变换的图像置乱算法进行多次加密, 本实验加密次数T=6次。加密结果如图3所示。

3. 嵌入水印

将经Arnold变换加密后的水印图3 (c) 采用叠加的方法嵌入到分数傅立叶变换并排序的矩阵中, 水印嵌入强度alfa=0.1, 得到DFRFT域含水印的幅度谱矩阵图4。

(三) 仿真实验

由于含水印图像可能会受到很多种方式的攻击, 侵权者将试图从含水印图像中检测并除去水印。对算法进行攻击测试是对水印鲁棒性检测的一种必要手段, 一个好的水印算法必须经过各种攻击测试才能对之做出客观的评价。为了分析这种水印算法抵抗攻击的性能, 本文对几种典型的攻击形式进行了数字仿真实验。

1. 加入噪声攻击后的实验结果

2. 剪切攻击后的实验结果

3. JPEG压缩处理后的实验结果

(四) 结束语

本文提出一种基于分数傅里叶变换的图像水印算法, 分数傅里叶变换具有时域和空域双域表达能力, 结果表明该算法能够很好地隐藏水印信息, 并在常见的攻击如噪声、剪切和JPEG压缩下具有较好的鲁棒性, 因而可实现数字图像的版权保护。

参考文献

[1]刘正君, 赵海发, 朱邦和.分数傅里叶域数字水印算法[J].光子学报, 2003, 32[3]:332-335.

[2]田震, 张晓娟.数字水印及其鲁棒性研究[J].测绘学院学报, 2004, 21[1]:54-56.

[3]袁占亭, 张秋余, 陈宁.数字水印的鲁棒性分析与研究[J].计算机工程与设计, 2005, 26[3]:614-616.

[4]周四清, 余英林.频域图像水印算法模型及其分析[M].北京:科学出版社, 2002.

多参数离散分数傅里叶变换的应用 篇9

关键词：离散分数傅里叶变换,加密,译码,模型

自1965年J.W.Tuky和T.W.Coody提出快速傅里叶变换算法以来,傅里叶变换迅速的渗透到现代科学技术的很多领域,并且发挥着越来越重要的作用,然而随着深入的研究和广泛的使用[1,2,3,4,5],傅里叶变换在处理一些问题时的局限性逐渐显现出来,于是人们对傅里叶变换进行各种推广,提出了小波变换、时频分析、Gabor变换、线性调频小波变换等新的信号分析理论。作为傅里叶变换的推广,V.Namias于1980年首先提出了分数傅里叶变换的概念,作为一种全新的信号时频分析工具,近年来,分数傅里叶变换(FRFT)变换已经在微分方程求解,量子力学,衍射理论和光学传输、光学系统和信号处理等方面有了广泛的应用[6,7,8]。其中将傅里叶变换应用于图像处理是分数傅里叶变换的重要应用之一,通常包括图像加密和数字水印技术。文中提出一种多参数分数傅里叶变换在图像加密解密中的应用方案,具有算法简单,容易实现,保密性高的特点。

傅里叶变换是一种线性算子,它提供了一种从时域到频域的转换方法,可以看作是从时间轴旋转π/2到频率轴,则如果找到一个可以作任意角度α=pπ/2旋转的算子就可以将研究对象当作维数更高的对象来处理。将离散分数傅里叶变换引入到图像加密解密中可以提高数据的安全性。

1 离散分数傅里叶变换的定义

离散分数傅里叶变换(DFRFT)是基于DFT矩阵的特征分解定义的,它具有两个主要特征:首先,离散分数傅里叶是具有一个附加阶参数的DFT的推广,并且具有FRT的所有性质;其次,DFR-FT的输出类似于连续FRT的采样结果。

FRT定义为

其中变换内核

式中α=απ/2,Ψn(t)是第n阶连续HermiteGaussian函数。

离散分数傅里叶变换可以由连续分数傅里叶变换得出,根据Candan等人提出的离散分数傅里叶变换的定义为

X(n)是带有N个矢量元素的输入信号,Fα是变换核矩阵,α是分数阶,当α=1时,离散分数傅里叶变换就退化成DFT,即离散傅里叶变换是离散分数傅里叶变换的特殊形式。求解离散分数傅里叶变换的关键是构造变换核矩阵F的特征分解。SooChang Pei等人将离散分数傅里叶变换核矩阵定义为

式中T表示矩阵的转置,当N为奇数时,矩阵V=[V0|V1|…|VN-2|VN-1],当N为偶数时,矩阵V=[V0|V1|…|VN-2|VN],Λ是一个对角矩阵,其对角线上的元素是V中的每列特征向量Vk的特征根,

下面来求解V和Λ。由于分数傅里叶变换和傅里叶变换具有相同的特征函数,因此可以用离散傅里叶变换的实变换矩阵S求出离散分数傅里叶变换的特征向量。实际上矩阵S并不是离散傅里叶变换的变换核矩阵F,然而,矩阵S和F满足交换律,SF=FS,因此S的特征向量也就是F的特征向量。由于矩阵S是对称的,特征矢量都是实数并且相互正交,这些特征矢量构成了离散分数傅里叶变换的Hermite-Gaussian多项式,因此可以由矩阵S得到的特征矢量矩阵V,V由N列特征向量{V1,V2,…,VN}组成。

已知N×N DFT矩阵F的特征值当N为偶数时为{(-i)n:n=0,…,N-1},当N为奇数时为{(-i)n:n=0,…,N-2,N}。因此,N×N DFT矩阵F定义为

DFT矩阵F有4个不相同的特征值:1,-1,j,-j。现在定义一个三角矩阵S,其非零元素为

矩阵F和S具有相同的特征向量,但特征根不同,根据对矩阵F的特征分解,可以将第α阶DFRFT矩阵定义为

已知连续分数傅里叶变换的特征值记为λk=exp(-jkπ/2),对于离散分数傅里叶变换的特征值与连续分数傅里叶变换的特征值相同,仅仅是k取有限值k=0,1,…,N。这些特征值构成一个的对角矩阵A如下

基于离散傅里叶变换的离散分数傅里叶变换的特征向量和特征值方法产生的定义不是唯一的[12],对特征值和特征向量的不同选择,导致了离散傅里叶变换的不同定义形式。如果用不同的分数次幂替代DFT矩阵的特征值λk=exp(-jkπ/2),则将FRFT推广到了MPDFRFT。

N点N×N MPDFRFT矩阵定义为

式中N分别为奇数和偶数,diag(r1,r2,…,rN)表示N×N对角矩阵的对角元素是r1,r2,…,rN,是1×N型参数矢量构成的具有N个独立阶参数的MPDFRFT

为了简化表示,令

Λ是N×N对角矩阵DFT的特征值,则式(6)可以记作

则带有参数矢量的数据矢量N×1 MPDRFT可以由式(8)计算

2 MPDRFT域中进行图像加密解密模型

文中在MPDFRFT域中采用双自由度编码进行数字图像加密解密,加密解密过程如图1和图2所示。由于将2D-MPDFRFT的阶参数扩展成了附加密钥,因此,这种加密方法能够极大地提高数据安全性。

首先将1D-MPDRFT扩展到2D-MPDRFT,设有一幅N×M点阵的数字图像P,则图像P的带参数矢量(,)2D-MPDFRFT可以表示为

其中,分别是N点和M点MPDPFRFT矩阵。分别是1×N和1×M参数矢量。

用[exp(jα(n,nm))]和[exp (jβ(n,m))]表示两个N×M自由相位矩阵,式中α(n,m)和β(n,m)当0≤n≤N,0≤m≤M时均匀分布在[0,2π]内,并且α(n,m)和α(n,m)相互独立。则加密图像Q和输入图像P的关系可以表示为

式中C=AB表示矩阵A和B矩阵的对应元素相乘,则加密图像Q的复共轭可以表示为

其中假设输入图像P非负实数。因此,解密图像R可以表示为

在以上加密过程中,参数矢量和自由相位码构成了MPDFRFT域中双自由相位编码加密的密钥。

3 计算机仿真

对上述加密解密过程在Matlab软件中进行仿真,设以下加密解＿密过程中采用同一个随机相位矩阵,用和表示在加密解密过程中的参数矢量。图3(a)表示用来加密的256×256的点阵原始图像,图3(b)表示按照上述加密原理进行加密在二维多参数离散傅里叶变换域内的双自由相位加密后的输出图形,其中1×256的加密矢量互独立且随机的选自[0,2]内。然后用正确的矢量参数按照上述原理进行解密,解密输出如图3(c)所示,从图中可以看出输出图像和原图像基本一样。如果解密矢量错误将得不出正确的输出图形。图3(d)表示在错误的解密矢量下解密结果。

这里用均方误差来表示这种算法的可靠性,均方误差可以定义为[10]

式中N×M表示图像的大小,I1(i,j)和I2(i,j)分别表示原图像和解密图像。对于任何一个分数阶独立误差为|△>0.005|时,MSE>2 000,得到的解密图像人眼将很难辨认与原图像的差异。

采用多参数分数傅里叶变换加密的另外一个重要特征是抗攻击性强。设加密参数矢量和误差为均匀分布,为了获得较高的抗攻击性,参数误差一般要求较小,这里假定为[-0.01,0.01],在已知其他两个加密矢量和以及自由相位密钥的情况下,根据多参数分数傅里叶变换的参数的周期性可知对于和的第阶元素,其概率为min(0.02/(4/k),1),因此成功破译的概率低于P=[(1/(4/0.02)2,(2/(4/0.02)2,…,(166/(4/0.02)2],如果采用PC机来破译,大约需要6.0×10134年。

4 结束语

一种基于傅里叶变换的形状描述方法 篇10

关键词：傅里叶描述符,闭合边界,多尺度,形状分类

0 引言

二维目标的形状分析与识别在图像分析、计算机视觉和目标识别等领域得到了广泛的发展, 广泛应用于工业监控、染色体分类, 目标识别和场景分析等场合中。在日常生活中, 人的视觉系统可以轻松地根据物体的边界信息来识别它的形状, 但机器视觉很难达到这种程度。从生理学的角度, 人的视觉系统对图像信息的分析和识别就是多尺度的。目前, 用于形状描述的方法主要有:基于区域的描述方法、基于轮廓边界的描述方法和自回归模型。前者利用一定区域内像素的分布信息来描述形状[1,2,3,4,5,6], 信息量大, 计算过程较复杂。而后者主要是利用物体形状的边界信息, 比较直观。本文将采用轮廓边界的描述方法。最常用的轮廓边界描述方法有链码, 傅里叶描述子和曲率尺度空间 (CSS) 描述方法。链码及其相关技术对噪声比较敏感[2], 自回归模型很难选择合适的自适应参数[3], CSS方法计算较复杂, 而且只考虑了边界的局部信息[4]、傅里叶描述的方法概念简单, 便于理解, 而且考虑了边界的全局和局部信息, 对噪声不敏感[5]。本文在分析了有关形状描述方法的基础上, 结合多尺度方法和傅里叶描述符, 提出了一种多尺度的傅里叶描述方法, 针对该方法提取了形状特征, 并进行了形状分类的研究和分析。

1 目标形状特征提取

1.1 物体边界函数的提取

物体的边界函数是使用某种特征函数来描述形状的边界线。位置函数是以边界线上的某个起始点, 沿逆时针或顺时针方向跟踪边界线, 以边界线上的点到起始点的弧长的函数。对二维目标进行轮廓跟踪, 得到了边界点对边界中心的坐标函数s (k) =[x (k) -Ox, y (k) -Oy], 其中k=0, 1, 2, …, N, N是边界的长度, (Ox, Oy) 是目标的中心。这样的边界函数是具有平移不变性的。然后对闭曲线s (k) 使用尺度为σ的高斯核进行线性卷积, 得到不同尺度下的进化曲线s (σ, k) 。高斯分布具有良好的特性, 对称于均值且随远离均值而衰减, 便于微分和积。σ值较大时, 边界曲线较平滑, 使得边界曲线呈现在较粗的尺度空间。图一显示了不同尺度的边界函数。

再有, 每对坐标对可以看成一个复数, 即:

对于复数序列, x轴表示实轴, y轴作为虚轴。这种表示法的优点是它将一个二维问题简化成一个一维问题, 这样就得到了一维的复数曲线簇s (σ, k) 。

1.2 不变性特征的提取

显然, s (k) 是具有平移不变性的, 经过高斯核线性卷积的边界曲线簇s (σ, k) 也具有平移不变性, 它不具有旋转不变性和尺度不变性, 而且对边界的起始点敏感。考虑到离散傅里叶变换具有的循环移位特性[8], 即序列的周期移位不改变序列频谱的幅值特性, 可以将边界曲线簇s (σ, k) 进行离散傅里叶变换 (DFT) , 得到不变性特征, 也就是将边界函数的多尺度描述转换到频域进行分析, 如公式 (3) 所示。

其中n=0, 1, 2, …, N-1。令不同尺度因子σ下的不变性特征FD (σ) 表示为[7]:

对于不同尺度下的特征, 令F={FD (σ) │σ奂A}, 其中A={σ1, σ2, …σR}可以很好地描述形状。对于频域的z (n) , 它的高频部分描述信号的细节特征, 低频部分描述信号的整体特征, 所以可以取FD (a) 的一个包含L个低频信号的子集来近似地描述边界信号。

由于边界函数与高斯函数卷积, 本身就有滤噪的特性。另外噪声的频谱特性主要反映在高频部分, 对低频部分影响较小, 而对形状的描述主要是通过边界序列的低频部分, 低频部分对噪声不太敏感。所以DFT得到的不变性特征, 其低频部分受噪声影响很小, 它具有良好的抗噪性能。而且DFT也有快速算法FFT, 计算过程耗时短。

1.3 形状的相似性

对于提取出傅立叶描述特征F的两个形状A和B, 通过计算它们之间的相似程度来识别是否属于同一类。由于傅立叶变换的各频率分量互相正交, 可以采用相关性函数计算它们傅立叶描述子F间的形状相关性, AB间的相关性计算表示为:

其中E是数学期望, cov是协方差, D是方差, σ是标准差, X、Y分别是AB的傅里叶描述特征。如果R值越接近1, 表明AB两物体相关程度越高, 形状越接近。反之, 它们的差距越大。

2 实验结果及分析

图二所示的6片叶子属于同一形状, 但分别具有不同的尺寸、大小和旋转方向。A1、A2、A3具有相同的尺寸, 但是具有不同的旋转方向;A4、A5、A6与A1的尺寸、大小和旋转方向都不同。令L=24, A={1, 15, 30}来计算傅立叶描述特征, 通过本文的方法计算出来的相似性数据显示本方法有较好的平移、尺度和旋转不变性, 如表一所示。但是, 由于高斯核函数本身具有平滑特性, 对于同一形状在不同尺度下的边界函数, 它们的低频成分是相似的。所以对于图一中的 (a) 、 (c) 、 (d) , 尽管它们的形状相差很大, 但是它们的低频成分任然很相近, 无法区分开来。

图三所示的7个形状属于三类, B1、B2、B3是五瓣, B4、B5是六瓣, B5相对B4有些变形, B6和B5是十瓣, 它们的边界都有噪声。令L=24, A={1, 15, 30}来计算傅立叶描述特征, 它们的相似性数据如表二所示。B1B2B3, B4B5, B6B7分别属于一类形状, 类内之间的相似程度都比较高, 三种形状之间的形似程度都比较低。由于B1B2B3是五瓣, B6B7是十瓣, 这两种形状之间的相似性较B3B4的高。相似数据的大小近似反映了形状差别的变化, 而且本方法具有较强的抗噪性能。

本文从MPEG-7标准体系的形状图像数据库中选取了20个类, 共400张图片组成了图像库, 图四是图像库中包含的20种形状。下面用这种傅里叶特征来进行分类能力的测试。令尺度因子A={1, 15, 30}来计算傅立叶描述特征, 使用KNN算法进行分类, 然后用留一交叉验证法来评估分类法的准确率, 即每一形状轮流作为测试样本, 剩余的组成训练样本, 使用KNN算法进行训练, 通过分类法准确率来体现分类能力。本文用复坐标的傅里叶描述子与本方法进行对比。由于MPEG-7标准库中的形状相差比较大, 使得分类的平均准确率都比较好。从表三可以看出, 对于不同的L值, 本文方法的平均准确率较高, 也就是说本文方法提取的特征是能够有效地作为形状分类的特征。

3 结束语

本文使用数字曲线的离散傅立叶变换 (DFT) 来计算形状的傅立叶描述特征, 这种方法主要适用于处理具有封闭轮廓的二维物体。我们提取了多尺度的傅立叶描述特征并分析了它的平移、旋转和尺度不变的特性以及对不同形状进行分类的有效性, 实验结果表明这种方法比较有效。但是傅立叶描述子对于边界出现遮挡的情况效果不理想。这种方法可以与SVM、贝叶斯网络或者神经网络结合起来识别物体, 这样可以增强识别方法的鲁棒性, 这是以后进一步研究的内容。

参考文献

[1]S.O.Relkasim, M.Shridhar, and M.Ahmadi.Pattern recognition with moment invariants:A comparative study and new results.Pattern Recognition, 1991, vol.24:1117-1138.

[2]J.Iivarinen, M.Peura, J.Sarela, A.Visa.Compari-son of Combined Shape Descriptors For Irregular Ob-jects[J].The8th British Machine Vision Corcfer-ence, 1997, 2:430-439.

[3]Sekita I, Kurita T, Otsu N.Complex autoregres-sive model for shape recognition[J].IEEE Trans on PAMI, 1992, 14 (4) :489-496.

[4]D.Zhang, G.Lu.A Comparative Study of Curva-ture Scale Space and Fourier Descriptors for Shape-Based image Retrieval[J].Journal of Visual Communication and Image Representation, 2003, 14 (1) :41-60.

[5]D.S.Zhang and G.Lu.Shape Based Image Re-trieval Using Generic Fourier Descriptors.Interna-tional Journal of Signal Processing:Image Communi-cation, 2001, 8.

[6]B.M.Mehtre, M.S.Kankanhalli, W.F.Lee.Shape Me-asures for Content Based Image Retrieval:A Compari-son[J].Information Processing&Management, 1997, 33 (3) :319-337.

[7]Gonzalez, R.C..Digital Image Processing[M].Second Edition, published by Pearson Education.Inc.