基于QR分解算法的任意阶复矩阵求逆的DSP实现

2010年第23卷第4期

E l ec tron i c Sci &T ech /Apr 15,

2010

协议 算法及仿真

基于QR 分解算法的任意阶复矩阵求逆的DSP 实现

倪 涛, 丁海锋, 阮黎婷, 张志强, 赵前晟

(西安电子科技大学雷达信号处理国家重点实验室, 陕西西安 710071)

摘 要 常用的矩阵求逆方法不易于工程实现, 计算软件也无法装备到电子系统中去, 文中采用QR 分解算法实现了基于A DSP TS201S DSP 汇编语言的矩阵求逆。该方法克服了高阶矩阵不易求逆, 且效率不高的缺点, 实现了任意阶复矩阵的快速求逆运算。仿真结果表明, 对于16阶的复数矩阵, 该工程方法的计算效率能达到 s 级, 计算精度能达到10-4量级。

关键词 QR 分解算法; 任意阶复矩阵求逆; 汇编语言; TS201S

中图分类号 TP301 6 文献标识码 A 文章编号 1007-7820(2010) 04-099-03

I m ple m ent ation of A rbitrary O rder Co m plexM atrix Inversion Based

on the QR D ecom position A lgorit hm in DSP

N iTao , D i n g H a ifeng , Ruan L iting , Zhang Zhiqiang , Zhao Q iansheng

(State K ey Lab o f R adar S i gnal P r ocessi ng , X idia n Un i versity , X i an 710071, China)

Abstract M atri x i nversion m ethods co mm only used i n eng i neeri ng i s d iffic ult to achieve , and the ca lc ulation soft ware can not be equ i pped i n an electronic syste m. The QR dec o mposition algorith m based on ADSP TS201S DSP

asse mbly language is use d to realize matri x i nversi on i n th i s paper , thus reso l v i ng the proble m of diffi cult ca lc ulation of the h i gh or der matri x i nversion a nd lo w co mputat i onal efficiency ,

and realizi ng an ar b itrary order co m ple x m atri x

-4

i nversi on fast operation . The sm i ulation results sho w that the co m putational efficiency of t he eng i neeri ng approach can ac h i eve the m icrosecond leve l and the calc u l ation accuracy can reach 10order 16.

K ey words QR deco m position al gorit hm; ar bitrary order co mplexm atrix i nversi on ; asse m bl y langua ge ; T S201S

magn itude f or the co m ple x m atri x of

目前计算低阶矩阵逆矩阵方法主要有:定义法、伴随矩阵法、初等行列变换法、等价标准型法、H a m ilton-Ca ley 法、分块矩阵法等。而对于某些特殊矩阵求逆的计算方法则更多, 比如广义范得蒙矩阵求逆快速算法, Toeplitz 矩阵求逆快速算法,

[3]

H anke 矩l 阵求逆快速算法等。但是这些算法只是从数学理论角度提出的, 低阶矩阵求逆容易实现, 阶数过高则难以人工计算出来。对于实际的工程其性能和精度难以实现。

某些计算软件, 如M atlab , 虽然能实现任意阶实数矩阵或者任意阶复数矩阵的求逆运算, 并能达到较高的精度, 但其实时性能却不理想。由于不能实现将计算软件嵌入到实时性要求高的嵌入式电子系统中, 因此而引起的在电子系统中无法计算矩阵求逆或者满足实时性要求的缺陷也越来越明显。文中由此提出如

[2]

[3]

[1]

何在实时性能要求高的电子系统, 如DSP 和FPGA 混合系统中, 实现任意阶矩阵的快速计算, 且还能满足实时性的要求。

工程中常用算法有:全选主元高斯-约当法、QR 分解算法等。这些算法工程实现后的计算效率也有较大差别, 所以对于实时性要求高的场合, 采用一种能实现快速计算的算法, 并将其用执行效率较高的汇编语言实现是重要的。

对于像数字波束形成(DBF) 系统而言, 需要解决的工程问题是如何实时计算出通道快拍数的相关矩阵并求其逆, 此算法实现与否将直接决定整个系统的成功。文中基于AD I 公司的ADSP TS201S 数字信号处理芯片, 实现了任意阶复数矩阵的求逆的快速运算, 在板调试中数值精度和效率完全满足系统性能的要求。

收稿日期:2009 09 07

1 QR 分解算法

前面提到的一些常用计算逆矩阵方法, 在工程实现中并不都是可行的, 文中在简要介绍采用QR 分解

99

作者简介:倪涛(1984-), 男, 硕士研究生。研究方向:雷达信号处理与高速电路设计。

协议 算法及仿真

算法求逆

[4]

倪涛, 等:基于QR 分解算法的任意阶复矩阵求逆的D SP 实现

之后将会详细分析如何用DSP 汇编语言

n ! n

和33 6GB /s的内部存储器宽度。若发挥其单指令多数据的特点, ADSP TS201S 可以提供48亿次40位MAC 运算或12亿次80位MAC 运算。2 2 汇编语言编程

主要计算过程分5个部分:输入数据、计算相关矩阵、QR 分解、计算R 的逆、计算Q 的共轭转置、计算相关矩阵的逆。

程序分为主程序和若干个子程序, 主程序描述整个程序的框架, 各个功能主要由相应的子程序来实现。子程序主要有:(1) 相关矩阵的计算子程序; (2) QR 分解子程序; (3) 计算QR 分解中R 矩阵的逆子程序;

(4) 计算QR 分解中Q 矩阵的共轭转置;

(5) 矩阵相乘子程序。

程序流程图, 如图1

所示。

实现这种算法。

定理1 假设A C 可逆, 则A 可以分解为

H

A =QR, 其中, Q 为酉阵, 即QQ =E, R 是上三角阵;

定理2 假设R =(R ij ) n ! n 是上三角阵, R ij =0, 当i >j 时, 并且R ji ∀0, 1#i #n , 则R 可以通过一下算法得出:

(1) ij =0, 当i >j 时, 即R 也是上三角阵; (2) kk =R kk , 1#k #n; (3) k, k+m=-矩阵, 则Q

-1

-1-1

-1

-1

=( ij ) n ! n

∃R

j=1H

n

k, k+j

k+j , k+m/R kk , 1#k #n -m 。

定理3 假设A 可逆, 且A =QR, 其中Q 是酉

=R Q , Q 是Q 的共轭转置。

H

2 DSP 汇编语言实现求逆算法

文中是基于ADSP TS201S DSP 实现QR 分解求逆算法的, 对于工程实现重要的两个性能指标是:求逆计算的精度和效率。对于DBF 系统, DSP 负责计算权值, 而用于DSP 计算权值的时间不能超过1 s 。计算权值的核心就是计算通道矩阵的逆矩阵, 如何尽可能的缩短逆矩阵的计算时间, 直接影响整个系统的性能。

2 1 ADSP TS201S 的主要性能指标

2003年, AD I 公司发布了T iger SHARC 的新产品

[5]

, 是ADSP TS201S , ADSP TS202S , ADSP

图1 程序流程图

TS203S 。ADSP TS201S 的最高时钟频率高达

600MH z , 而ADSP TS202S 和ADSP TS203S 的最高时钟也能达到500MH z 。本项目采用最高时钟高达600MH z 的ADSP TS201S 芯片。

ADSP TS201S 的性能如下:

(1) 600MH z 的时钟频率, 1 6ns 的指令周期; (2) 双运算模块, 每个计算块包含1个ALU, 1个乘法器, 1个移位器, 1个寄存器组和1个CLU;

(3) 双整数ALU, 提供数据寻址和指针操作功能;

(4) 集成I/O端口, SDRAM 控制器, 可编程标志引脚, 2个定时器和定时器输出引脚。

由此可见, 性能高的静态超标量处理器ADSP TS201S , 将较宽的存储器宽度和双运算模块组合, 使DSP 每周期能够执行多达4条指令, 24个16位定点运算和6个浮点运算。4条相互独立的128位宽度的内部数据总线, 每条总线分别连接6个4MB 内部存储器块中的1个, 提供4字的数据、指令及I /O访问

为了在编程时容易书写, 某些矩阵的字母表示与前文算法中的字母表示有所不同, 现将主要差别列表, 如表1所示。

表1 各字母表示的矩阵名称

程序中的矩阵名表示方式

A Q D E F

D (1):计算最后, D 被A -1覆盖

算法中的矩阵名表示方式

相关矩阵A QR 分解中的Q 矩阵QR 分解中的R 矩阵

中间变量R -1

相关矩阵的逆矩阵A -1

需要注意, 汇编程序中不能用R 作为QR 分解中的矩阵名称的表示字母, 因为R 是汇编中的关键字, 否则会引起冲突。

在用汇编语言实现QR 分解算法过程中, 一定要考虑到具体芯片载体的特定功能。如文中所用ADSP

倪涛, 等:基于QR 分解算法的任意阶复矩阵求逆的D SP 实现

协议 算法及仿真

TS201S 具有4级流水功能; 芯片内部有两个运算块, 可同时进行运算及特殊的单指令多数据操作。所以要综合运用这些特殊的功能, 才能实现程序效率的最大化。

从分析程序中各部分所用时间的比例可以看出, 所有子程序中计算相关矩阵的子程序占用时间最多, 约为34 19%。如图2

所示。

图3 计算精度%绝对误差

3 结束语

在宽带雷达和窄带雷达的应用日益增多的今天, 雷达数字信号处理系统部分需要承担波束合成的任务, 而波束合成的前提就是要根据采样数据进行计算权系数。权系数计算的核心部分就是要计算相关矩阵的逆矩阵。由此可知, 相关矩阵的逆矩阵计算的效率和精度将直接影响整个雷达系统的性能, 甚至决定雷

图2 程序运行时间图

达系统的设计成功与否。

将基于QR 分解算法实现矩阵求逆的DSP 实现应用于常见的DSP 和FPGA 混合电子系统中, 可以实现矩阵求逆程序部分的模块化和通用化。任意阶复数矩阵逆矩阵的工程实现, 将极大的改善宽窄带雷达中数字信号处理单元的效率。

基于QR 分解算法求矩阵逆的工程实现中也有一些不足之处, 比如能否将汇编程序继续优化, DSP 中特殊的多级流水和4字取功能是否充分利用等等。今后的工作就是寻找更加优化的算法, 在此基础上实现

净求逆值74614

所以如何减少占用时间最多地求相关子函数的运行时间, 是实现程序继续优化的关键所在。2 3 程序的效率和计算的精度

程序开始读入的矩阵是16! 64, 首先计算出其相关矩阵, 再求其逆。

现将软仿真运行时时钟计数器指示值列, 表2所示。

表2 软仿真计数器数值起始值计算相关后值求逆后值

核时钟计数

643

39433

114047

DSP 汇编程序的最优化, 如此将会再次提高工程实现时计算的精度和效率, 也必将再次提高整个系统的性能。参考文献

[1] 陈逢明. 逆矩阵的若干求法[J].福建商业高等专科学

校学报,

2006(6) :117-119.

2004,

36

[2] 陆全, 任学明. 一类广义范得蒙矩阵求逆的快速算法

[J].西安建筑科技大学学报:自然科学版, (3):

369-371.

相关矩阵的阶数是16阶, 软仿真运行求相关的时钟周期数为38790, 求相关矩阵逆矩阵的周期数是

74614, 即求一个16阶方阵的逆矩阵所需要的核时钟周期数是74614个。如果ADSP TS201S 运行在核时钟为600MH z 速度时, 只需1 24357e-4s , 即约124 s 。

本程序完全用汇编语言编写, 可实现计算任意阶数矩阵的逆, 并且在16阶下实测的计算时间约为124 s , 完全能满足实时性的需要。

将ADSP TS201S 计算的结果与理想结果进行比较, 计算其绝对误差, 如图3所示。

由图3可以看出, 基于QR 分解算法ADSP TS201S 实现的矩阵求逆, 误差水平小到10量级。

-4

[3] 游兆永, 路浩. T teplitz 矩阵、H ankel 矩阵求逆的快速算法

[J].西安交通大学学报, 1990, 24(6):131-134. [4] 王建锋. 求逆矩阵的快速算法[J].大学数学,

20(1):121-122.

[5] 刘书明, 罗勇江. ADSP TS20XS 系列D SP 原理与应用设

计[M].北京:电子工业出版社, 2007.

2004,

101