第5期2000年5月
电 子 学 报
ACT A E LECTRONICA SINICA V ol. 28 N o. 5
M ay 2000
用卷积运算实现反卷积
刘明亮, 蔡永泉, 饶 敏, 刘 弟
(北京工业大学计算机学院, 北京100022)
摘 要: 本文首先简介了用卷积运算实现反卷积的原理、方法和实例; 其次, 提出了该方法的病态问题, 并给出
了解决的办法; 最后, 给出了这种方法的估计误差及实例.
关键词: 卷积; 反卷积; 病态问题; 误差中图分类号: T N911172 文献标识码: A 文章编号: 037222112(2000) 0520111202
Deconvolution Via Convolution Operation
LI U Ming 2liang ,C AI Y ong 2U Di
(Beijing Polytechnic Univer sity , , )
Abstract : olution operation are introduced at first briefly. The ill 2posed problem of method is given next. Finally ,the estimate error and exam ple of this technique is given.
K ey words :conv olution ;deconv olution ;ill 2posed problem ;error
1 引言
根据卷积定理, 有关信号恢复和系统辨识的反卷积问题
借助FFT 可以很方便地求解, 即
x (n ) =IFFT [X (k ) ]=IFFT [Y (k ) /H (k ) ]
(1)
h (n ) =IFFT [H (k ) ]=IFFT [Y (k ) /X (k ) ]
但是, 由于X (k ) 或H (k ) 有零点, 形成所谓的病态反卷积. 为此, 人们进行了不少有益的研究工作[1, 2]. 用卷积运算实现反卷积的方法直观上可以绕过X (k ) 或H (k ) 的零点[3]. 但是, 我们发现这种方法同样存在所谓的“病态问题”; 对此提出相应的解决方法. 另外, 在存在测量误差或噪声污染时, 用卷积运算实现反卷积这一方法与其它算法一样将产生很大误差, 有时甚至造成难以接受的反卷积结果.
成不必要的位移. 通常数组下标是0开始的, 即为0, 1, …, N
δ(0) , 更为一般化. -1, 这样x p (n ) =[A , 0, 0, …, 0]=A
例1:当系统的输入和输出分别为x (n ) =[1, 3, 2, 1, 0, 0, 0, 0], y (n ) =[2, 13, 31, 38, 31, 14, 4, 0], 用卷积运算实现反卷积的方法求h (n ) .
解:依原理可得
x 3(n ) =[-875, 0, 0, 0, 0, 0, 0, 0]=-875δ(0)
y 3(n ) =[-1750, -6125, -5250, -3500, 0, 0, 0, 0]
所以h (n ) =y 3(n ) /-875=[2, 7, 6, 4, 0, 0, 0, 0].显然, 很容易验证y (n ) =x (n ) h (n ) .
3 用卷积运算实现反卷积的病态问题
用卷积运算实现反卷积时, 我们发现有时会出现x p (n ) =[0, 0, …, 0].这样也无法求得h (n ) =y p (n ) /x p (n ) , 产生所谓“病态问题”. 分析表明, 出现所谓病态问题是由于x p -1
33(a 0) =A , x p -1(a N/2) =±A , x p -1(a 0) =A , x p -1(a N/2) = A
2 用卷积运算实现反卷积
这一方法和原理, 文献[3]已详细论述. 在这里仅给出必
要的修正. 通常的圆周卷积的公式为
N -1
y (n ) =h (n ) x (n ) =
m =0
∑h (m ) x (n -
m ) N R N (n ) ,
n =0, …, N -1 (2)
如果对式(2) 进行操作使其变为
δ(0) y p (n ) =h (n ) x p (n ) =h (n ) A 则
h (n ) =y p (n ) /A
(3)
δ(0) . 但是, 由于文献[3]中数组下其中, p =log 2N , x p (n ) =A
δ(N ) , 造标是从1开始, 从而推成x p (n ) =[0, 0, 0, …, A ]=A
造成的. 下面我们举例来说明解决这类问题的方法.
例2:已知系统的输入和输出分别为x (n ) =[1, 1, 1, 1, 0, 0, 0, 0], y (n ) =[1, 7, 4, 3, 3, -3, 0, 1], 求系统的单位抽样响应h (n ) .
解:依原理可求得
3
x (n ) =[1, -1, 1, -1, 0, 0, 0, 0]x 1(n ) =[1, 0, 1, 0, -1, 0, -1, 0]
y 1(n ) =[3, 7, -3, 5, -3, -7, 3, -5]
收稿日期:1999202222; 修订日期:1999205226
112
3
x 1(n ) =[1, 0, -1, 0, -1, 0, 1, 0]
电 子 学 报2000年
-011, 0104, 01015, -0103].与例1h (n ) =[2, 7, 6, 4, 0, 0, 0, 0]
x 2(n ) =[4, 0, 0, 0, -4, 0, 0, 0]
y 2(n ) =[0, 24, -12, -4, 0, -24, 12, 4]
3
x 2(n ) =[4, 0, 0, 0, 4, 0, 0, 0]
相比, 可知当|e x |=|e y |=1%, (e ′h ) m =10%.按式(9) 计算可得e h =16%, 且(e ′h ) m
5 结束语
本文对用卷积运算实现反卷积这一方法进行了一点修正, 并给出了误差分析. 特别是指出了本方法仍存在所谓的病态问题, 在实际运用中应注意这一点. 参考文献
[1] 舒勤, 张有正. X (k ) 有零点的卷积反演的DFT 算法. 电子学报,
1990,18(3) :83～89
[2] 舒勤, . 再论X (K ) DFT 算法. 电子
x 3(n ) =[0, 0, 0, 0, 0, 0, 0, 0]y 3(n ) =[0, 0, 0, 0, 0, 0, 0, 0]
从例2可以看出, 由于x 3(n ) =[0, 0, …, 0], 产生所谓病
态问题, 不能得到h (n ) . 在例2中, 由于x 2(0) =4, x 2(4) =
-4, x 23(0) =4, x 23(4) =4, 所以x 3(n ) =x 2(n ) x 23(n ) =[0, 0, …, 0].
解决这类病态问题, 不必进行到x p (n ) 这一步, 由于在x p -1(n ) 中, 出现|x p -1(a 0) |=|x p -1(a N/2) |, 其它元素为0. 而且在y p -1(n ) 中, 出现|y p -1(a 0) |=|y p -1(a N/2) |, |y p -1
(a 1) |=|y p -1(a N/2+1) |, …, |y p -1(a N/2-1) |=|y p -1(a N -1) |,
,1992(12) [3] . . ,1985,13(4) :9～14and G Deconv olution of time domain
presence of noise. NBS T ech. 1047,1981
[5]A. Bennia and S. M. Riad. An optim ization technique for iterative fre 2
quency 2domain deconv olution. IEE T rans Instr. M eas. 1990,I M 239(2) :358～362
[6] T. Dhaene ,L. M artens. Extended Bennia 2Riad Criterion for iterative
frequency 2domain deconv olution. IEEE T rans. Instr. M eas. Apl. 1994, I M42(2) :176～180
[7] 刘明亮, 白瑞林. 测量误差对反卷积影响的分析及其滤波器的
这时将x p -1(n ) 分解为x p -1(n ) =[a 0, 0, …, 0]+[0, …, ,
a N/2, 0, …, 0]=x ′p -1(n ) +x ″p -1(n ) , 将y p -1p -1
(n ) =[a 0, a 1, …, a N/2-, 0, ], a N/21, , a N -1]=y ′p -1(n ) y p -(这样就有h (n ) =y 1x ′p -1a 0) .
, 就是圆周卷积的数据重叠问题. 当两个长为N 的序列线性卷积时, 结果序列长为(2N -1) . 所以用卷积计算反卷积时, 当y (n ) 长为(2N -1) , x (n ) 长为N 时, 那么h (n ) 应为N 长. 例如在例2中, x (n ) 中有4个非零元素, y (n ) 为8元序列, 这说明h (n ) 中应有5个非零元素. 所以用圆周卷积计算就会有一个端数据重叠, 为了避免这一问题, 可以适当补零. 另外考虑到序列总长度为N =2p (为便于FFT 计算) . 所以例2中的x (n ) 、y (n ) 补零后的长度为16. 因此例2的结果如下:y p -1(n ) =[4, 24, -12, -4, 4, 0, 0, 0, 0, 0,
(n ) =y p -1(n ) /x p -1(n ) =[1, 0, 0, 0, 0, 0, 0], x p -1(a 0) =4, h ′
6, -3, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 再将序列恢复到8元, 则h (n ) =[1, 6, -3, -1, 1, 0, 0, 0], 易证y (n ) =x (n )
h (n ) .
设计思想. 电子测量与仪器学报,1994,8(3) :44～
49
刘明亮 北京工业大学计算机学院教授,
1943年生,1967年毕业于大连理工大学无线电技
术专业,1994———1995为日本东京工业大学高级访问学者. 合作出版书六本, 其中一本被评为全国优秀教材; 完成的科研项目多次受到部级奖励, 已发表论文二十多篇, 目前主要研究方向有数字滤波与数字反卷积,DSP 的开发与应用, 数
字信号处理的理论与方法.
蔡永泉 北京工业大学计算机学院副教授、工学博士,1956年生人1992年西北工业大学计算机应用硕士毕业,1997年获得北京农业工程大学电气化与自动化工学博士。现主要从事计算机网络及数字信号处理方面的研究。在国内外发表十几篇论文。
饶 敏 北京工业大学计算机学院副教授,1961年生,1983年毕业于北京师范大学无线电电子学系. 先后完成科研项目八项, 其中包括“七五”国家重点科技攻关项目, 并获“七五”科技攻关重大成果集体荣誉奖, 发表论文9篇. 现从事计算机应用的研究工作.
4 误差估计
测量误差或噪声对各种算法的反卷积的影响都很
大[4～7]. 同样测量误差对卷积算法实现反卷积这一方法的影响也不可忽视. 下面我们来分析这一情况.
如上所述, 由x (n ) 、y (n ) 经过圆周卷积计算出x p (n ) 、y p (n ) 要经过p =log 2N 步. 设x (n ) 的误差为e x , y (n ) 的误差为e y , 则x p (n ) 的误差为2p e x , y p (n ) 的误差为(2p -1) e x +e y , 所以h (n ) =y p (n ) /x p (n ) 的误差为:
e h =[(2p +1-1) e x +e y ]m =(2N -1) |e x |+|e y |
(4)
例3:设例1的测量误差|e x |=|e y |=1%, 即已知x (n ) =[0199, 3103, 1198, 1101, 0, 0101, 0101, 0], y (n ) =[2102, 12187, 30169, 38138, 31131, 14, 3196, 0], 用卷积运算实现反卷积的方法求得h (n ) =y 3(n ) /A =[21038, 6176, 6125, 41026,
刘 第 北京工业大学计算机学院硕士研究生,1974年生人,
1997年毕业于哈尔滨工业大学电子与通信工程系, 现从事数字信号
处理的研究工作.