唐毅鋆,陳穎頻,陳 惠,喬嘉琪,陳振雕,李一凡
(1.漳州職業(yè)技術(shù)學(xué)院電子信息學(xué)院,福建 漳州 363000;2.閩南師范大學(xué)物理與信息工程學(xué)院,福建 漳州 363000)
在圖像處理和計算機(jī)視覺的眾多應(yīng)用中廣泛存在具有列循環(huán)結(jié)構(gòu)或行循環(huán)結(jié)構(gòu)的循環(huán)矩陣,例如圖像處理中的差分矩陣[1-2]和相關(guān)跟蹤濾波中由基向量循環(huán)移位得到的預(yù)測樣本矩陣[3-8]。由于循環(huán)矩陣在空間上存在明顯的循環(huán)移位關(guān)系,其信息存在冗余性,因此若直接在空域上處理循環(huán)結(jié)構(gòu)矩陣會導(dǎo)致產(chǎn)生較高的計算復(fù)雜度。為解決這一問題,學(xué)者們將傅里葉變換引入循環(huán)矩陣,并將大型的矩陣相乘和求逆運算轉(zhuǎn)換為向量的點乘與點除運算,其核心原理就是循環(huán)矩陣的對角化性質(zhì)。然而,傳統(tǒng)的循環(huán)矩陣對角化性質(zhì)的證明較難理解和掌握,為此,本文從卷積視角提出一種全新的列循環(huán)矩陣對角化證明方法,進(jìn)而證明行循環(huán)矩陣的對角化性質(zhì)。
循環(huán)矩陣對角化問題描述如下:
(1)
根據(jù)復(fù)數(shù)信號的周期不變性,
ekj2π/N=ekj2π/N·e-2πkj=e-kj2π(N-1)/N.
(2)
同理,
ekj2π(N-1)/N=ekj2π(N-1)/N·e-2πkj=e-kj2π/N.
(3)
因此,
(4)
g(1,k)=x1+x0ekj(N-1)2π/N+…+x2ekj2π/N=ejk2π/Ng(0,k).
(5)
則
(6)
故有
(7)
z=x*y∈N×1,
(8)
將列循環(huán)矩陣C(x)乘以一個信號y,根據(jù)1.2節(jié)中介紹的離散卷積定義,則有
(9)
根據(jù)式(9),在兩邊乘以離散傅里葉變換矩陣FN,則有
FNC(x)y=FN(x*y).
(10)
首先分析式(10)右邊,根據(jù)卷積定理,空域卷積信號的頻譜將等于兩個信號頻譜的點乘[10-11],即
(11)
然后分析式(10)左邊,將C(x)y看作C(x)INy,其中,IN∈N×N表示單位矩陣,可被拆分為則有
(12)
由于式(10)成立,必有
(13)
本節(jié)提供代碼1對列循環(huán)矩陣對角化性質(zhì)加以驗證,代碼如下:
%代碼1
clear all;
N=4;
DFT=zeros(N,N);
n=[0:N-1]; k=[0:N-1];
Wn=exp(-j*2*pi/N);
nk=n’*k;
DFT=Wn.^nk;
C=[-1 1 0 0; 0 -1 1 0; 0 0 -1 1; 1 0 0 -1];
x=[-1,0,0,1].’;
X=fft(x);
D=DFT*C*inv(DFT);
for i=1∶N
Y(i)=D(i,i);
end
Y-X.’
經(jīng)Matlab運行,代碼1運行結(jié)果約等于0。該結(jié)果表明,運用卷積定理進(jìn)行列循環(huán)矩陣對角化與通過傳統(tǒng)數(shù)學(xué)方法計算的誤差極小,可見式(7)成立。
證明 因(7)式恒成立,對該式兩邊同時轉(zhuǎn)置,則有
(14)
(15)
(16)
(17)
本節(jié)提供代碼2對行循環(huán)矩陣對角化性質(zhì)加以驗證,代碼如下:
%代碼2
N=8;
DFT=zeros(N,N);
n=[0∶N-1];
k=[0∶N-1];
Wn=exp(-j*2*pi/N);
nk=n’*k;
DFT=Wn.^nk;
Fn=DFT./sqrt(N);
C=[0 7 6 5 4 3 2 1;
1 0 7 6 5 4 3 2;
2 1 0 7 6 5 4 3;
3 2 1 0 7 6 5 4;
4 3 2 1 0 7 6 5;
5 4 3 2 1 0 7 6;
6 5 4 3 2 1 0 7;
7 6 5 4 3 2 1 0];
D=inv(Fn)*C*(Fn);
for i=1∶N
y(i)=D(i,i);
end
x=[0 7 6 5 4 3 2 1];
FTx=fft(x);
y-FTx
以上代碼中,C以(0, 7, 6, 5, 4, 3, 2, 1)作為第一行,經(jīng)循環(huán)移位擴(kuò)展得到行循環(huán)結(jié)構(gòu)矩陣。代碼2運行結(jié)果趨近于0,該結(jié)果表明,運用卷積定理進(jìn)行行循環(huán)矩陣對角化與通過傳統(tǒng)數(shù)學(xué)方法計算的誤差極小,可驗證式(17)成立。
在相關(guān)濾波跟蹤算法中,常將二維目標(biāo)面片X∈h×w拉成行向量xT∈1×N(N=h×w),并循環(huán)移位擴(kuò)展為行循環(huán)矩陣,以此建模目標(biāo)在空間中移動的各種可能性,如圖1所示。
(a)循環(huán)移位得到的預(yù)測樣本1
(b) 基樣本
(c)循環(huán)移位得到的預(yù)測樣本2
然后以嶺回歸的機(jī)器學(xué)習(xí)建模方法訓(xùn)練濾波器w∈N×1,使得濾波器回歸的響應(yīng)盡可能接近預(yù)先定義的教師信號N×1(該信號為二維高斯窗函數(shù)的拉列信號),如式(18)所示。
(18)
目標(biāo)函數(shù)在空域上對濾波器w進(jìn)行求導(dǎo):
(19)
則w的空域最優(yōu)解為
(20)
若樣本尺寸較大,在空域上直接求解濾波器w的復(fù)雜度很高。行循環(huán)矩陣的對角化性質(zhì)則可較好地解決這一問題。
將式(20)中的C(xT)對角化:
根據(jù)卷積定理可將式(21)轉(zhuǎn)換到頻域中計算,則有
(22)
其中,real表示取復(fù)數(shù)的實部。
圖2展示了利用式(20)和(22)設(shè)計濾波器得到的響應(yīng)圖,從圖2可以看出,兩種算法效果上完全等效。通過上述方法可將濾波器利用快速傅里葉變換加以計算,運算效率會有顯著提高。
(a)利用式(20)設(shè)計濾波器得到的響應(yīng)圖
(b)利用式(22)設(shè)計濾波器得到的響應(yīng)圖
本節(jié)提供代碼3以反映式(20)與(22)計算效率上的差別。代碼如下:
%代碼3
clc
clear all;
N=4096;
X1=[1∶1∶sqrt(N)];
X2=[1∶1∶sqrt(N)];
[X,Y]=meshgrid(X1,X2);
y=exp(-((X-sqrt(N)/2).*(X-sqrt(N)/2)+(Y-sqrt(N)/2).*(Y-sqrt(N)/2))/10);
C=zeros(N,N);
x=linspace(1,N,N);
C(1,:)=x.’;
for i=1∶N-1
C(i+1,:)=circshift(x,i).’;
end
tic
w=inv(C.’*C+0.1)*C.’*y(:);
toc
r1=C*w;
figure(1)
R1=reshape(r1,sqrt(N),sqrt(N));
imagesc(R1);
colorbar;
tic
w2=real(ifft(fft(x.’).*fft(y(:))./(fft(x.’).*conj(fft(x.’))+0.1)));
toc
r2=C*w2;
figure(2)
R2=reshape(r2,sqrt(N),sqrt(N));
imagesc(R2);
colorbar;
r1-r2
通過改變代碼3中N的數(shù)值,對比兩種算法的運行時間,如表1所示。
表1 式(20)與式(22)的耗時對比
從表1可以看出,當(dāng)N值較大時,以式(22)計算濾波器的運算耗時遠(yuǎn)遠(yuǎn)低于以式(20)計算濾波器的運算耗時。這是因為式(20)涉及大型矩陣C(xT)∈N×N的相乘運算和求逆運算,其復(fù)雜度高達(dá)(N3),而式(22)中用到的快速傅里葉變換和快速逆傅里葉變換復(fù)雜度僅為(Nlog2N),且向量的點乘運算復(fù)雜度僅為(N2),故式(22)復(fù)雜度為(N2+2Nlog2N)。因此,利用循環(huán)矩陣對角化將濾波器轉(zhuǎn)換為頻域求解的效率較高。
循環(huán)矩陣對角化性質(zhì)廣泛應(yīng)用于各類圖像處理和計算機(jī)視覺領(lǐng)域。然而,傳統(tǒng)的循環(huán)矩陣對角化證明方法往往晦澀難懂,需要用到大量數(shù)學(xué)技巧。為解決這一問題,本文從列循環(huán)矩陣的性質(zhì)C(x)y=x*y出發(fā),結(jié)合離散卷積定義和卷積定理,從卷積視角巧妙證明了列循環(huán)矩陣的對角化性質(zhì),避免傳統(tǒng)證明方法中涉及的復(fù)雜運算,證明思路易于理解,最后通過行循環(huán)矩陣與列循環(huán)矩陣的轉(zhuǎn)置關(guān)系進(jìn)一步討論行循環(huán)矩陣的對角化性質(zhì),并給出一組實驗以反映行循環(huán)矩陣對角化性質(zhì)在視頻目標(biāo)跟蹤中的重要應(yīng)用。