CSI筆記【12】：陣列信號處理及MATLAB實現（第2版）閱讀隨筆（四）

• • Chapter6 陣列信號處理常用MATLAB代碼
  • • • (1).LCMV波束形成MATLAB程序：
      • (2).LMS自適應波束形成MATLAB程序：
      • (3).root-MUSIC算法MATLAB程序：
      • (4).譜峰搜索傳播算子算法MATLAB程序：
  • Reference

Chapter6 陣列信號處理常用MATLAB代碼

(1).LCMV波束形成MATLAB程序：

“線性約束最小方差準則（LCMV）”在對有用信號形式和信號來向完全未知情況下，在某種約束條件下使陣列輸出的方差最小。

"LCMV"方法的代價函數可以表示為 $J(\omega )={\omega }^{H}R\omega$，約束條件為 ${\omega }^{H}a(\theta )=f$。取 $f=1$，得到最佳解為 $\omega =R^{?1}c[c^H{R}^{?1}c{]}^{?1}$。

clc; close all clear all; M=18; % 天線數 L=100; % 快拍數 thetas=10; % 信號入射角度 thetai=[-30 30]; % 干擾入射角度 n=[0:M-1]'; % 構造一個一維列矩陣 vs=exp(-1i*pi*n*sin(thetas/180*pi)); % 信號方向矢量 vi=exp(-1i*pi*n*sin(thetai/180*pi)); % 干擾方向矢量 f=16000; % 載波頻率 t=[0:1:L-1]/200; % 構造時間變量 snr=10; % 信噪比 inr=10; % 干噪比 xs=sqrt(10^(snr/10))*vs*exp(1i*2*pi*f*t); % 構造有用信號 xi=sqrt(10^(inr/10)/2)*vi*[randn(length(thetai),L)+1i*randn(length(thetai),L)]; % 構造干擾信號 noise=[randn(M,L)+1i*randn(M,L)]/sqrt(2); % 噪聲 X=xi+noise; % 含噪信號 R=X*X'/L; % 構造協方差矩陣 wop1=inv(R)*vs/(vs'*inv(R)*vs); % 波束形成 sita=48*[-1:0.001:1]; % 掃描方向范圍 v=exp(-1i*pi*n*sin(sita/180*pi)); % 掃描方向向量 B=abs(wop1'*v); % 求不同角度的增益 plot(sita,20*log10(B/max(B)),'k'); title('波束圖');xlabel('角度/degree');ylabel('波束圖/dB'); grid on axis([-48 48 -50 0]); hold off figure1 LCMV波束形成算法仿真結果

(2).LMS自適應波束形成MATLAB程序：

“最小均方算法（LMS）”采用迭代模式，在每個迭代步驟 $n$ 時刻的權向量加上一個校正量后，即組成 $n+1$ 時刻的權向量，用它逼近最佳權向量。

clear all close all clc M=16; % 天線數 K=2; % 信源數 theta=[0 30]; % 信號入射角度 d=0.3; % 天線間距 N=500; % 采樣數 Meann=0; varn=1; % 噪聲均值、方差 SNR=20; INR=20; pp=zeros(100,N);pp1=zeros(100,N); rvar1=sqrt(varn)*10^(SNR/20); % 信號功率 rvar2=sqrt(varn)*10^(INR/20); % 干擾功率 for q=1:100s=[rvar1*exp(1i*2*pi*(50*0.001*[0:N-1]));rvar2*exp(1i*2*pi*(100*0.001*[0:N-1]+rand))]; % 生成源信號A=exp(-1i*2*pi*d*[0:M-1].'*sin(theta*pi/180)); % 方向向量e=sqrt(varn/2)*(randn(M,N)+1i*randn(M,N)); % 噪聲Y=A*s+e; % 接收信號% LMS算法L=200;de=s(1,:);mu=0.0005;w=zeros(M,1);for k=1:Ny(k)=w'*Y(:,k); % 預測下一個采樣和誤差e(k)=de(k)-y(k); % 誤差w=w+mu*Y(:,k)*conj(e(k)); % 調整權向量end end% 波束形成 beam=zeros(1,L); for i=1:La=exp(-1i*2*pi*d*[0:M-1].'*sin(-pi/2+pi*(i-1)/L));beam(i)=20*log10(abs(w'*a)); end% 作圖 figure angle=-90:180/200:(90-180/200); plot(angle,beam); grid on xlabel('方向角/degree'); ylabel('幅度響應/dB'); figure for k=1:Nen(k)=(abs(e(k))).^2; end semilogy(en);hold on xlabel('迭代次數'); ylabel('MSE');

運行結果：

figure2 LMS自適應波束形成算法仿真結果 figure3 LMS自適應波束形成算法誤差圖

(3).root-MUSIC算法MATLAB程序：

求根"MUSIC"，即"root-MUSIC"算法是"MUSIC"算法的一種多項式求根形式，其基本思想是"Pisarenko"分解。相比于"MUSIC"算法，"root-MUSIC"算法無須譜峰搜索，降低了復雜度。

clear all close all derad=pi/180; % 角度->弧度 radeg=180/pi; % 弧度->角度 twpi=2*pi; kelm=8; % 陣元數 dd=0.5; % 陣元間距 d=0:dd:(kelm-1)*dd; iwave=3; % 信源數 theta=[10 20 30]; % 波達方向 snr=10; % 信噪比 n=200; % 采樣數 A=exp(-1i*twpi*d.'*sin(theta*derad)); % 方向向量 S=randn(iwave,n); % 信源信號 X0=A*S; % 接收信號 X=awgn(X0,snr,'measured'); % 添加噪聲 Rxx=X*X'; % 計算協方差矩陣 InvS=inv(Rxx); [EVx,Dx]=eig(Rxx); % 特征值分解 EVAx=diag(Dx)'; [EVAx,Ix]=sort(EVAx); % 特征值從小到大排序 EVAx=fliplr(EVAx); % 左右翻轉，從大到小排序 EVx=fliplr(EVx(:,Ix));% root-MUSIC Unx=EVx(:,iwave+1:kelm); % 噪聲子空間 syms z; pz=z.^([0:kelm-1]'); pz1=(z^(-1)).^([0:kelm-1]); fz=z.^(kelm-1)*pz1*Unx*Unx'*pz; % 構造多項式 a=sym2poly(fz); % 符號多項式->數值多項式 zx=roots(a); % 求根 rx=zx'; [as,ad]=(sort(abs((abs(rx)-1)))); DOAest=asin(sort(-angle(rx(ad([1,3,5])))/pi))*180/pi

運行結果：

DOAest=9.9925 20.1515 29.9966

(4).譜峰搜索傳播算子算法MATLAB程序：

譜峰搜索傳播算子算法利用方向向量和 $Q$ 矩陣的正交性構造譜函數，通過一維譜峰搜索，得到波達方向估計。相比于“MUSIC”算法，它無須對接數據的協方差矩陣進行特征值分解，降低了復雜度。算法的主要步驟如下：

根據接收信號構造協方差矩陣 $R$。

對協方差矩陣分塊，計算傳播算子 $P$。

構造 $Q$ 矩陣。

使 $\theta$ 變化，按照 $P_{PM}(\theta )=1/[a^H(\theta )Q{Q}^{H}a(\theta )]$ 計算譜函數，通過尋找譜峰來得到波達方向的估計值。

clear all close all derad=pi/180; radeg=180/pi; twpi=2*pi; kelm=16; % 陣元數 dd=0.5; % 陣元間距 d=0:dd:(kelm-1)*dd; iwave=3; % 信源數 theta=[10 20 30]; % DOA pw=[1 0.8 0.7]'; % 信號功率 nv=ones(1,kelm); % 歸一化噪聲方差 snr=20; % 信噪比 snr0=10^(snr/10); n=200; % 樣本數量A=exp(-1i*twpi*d.'*sin(theta*derad)); % 方向向量K=length(d); cr=zeros(K,K); L=length(theta); data=randn(L,n); data=sign(data); twpi=2.0*pi; derad=pi/180.0; s=diag(pw)*data; A1=exp(-1i*twpi*d.'*sin([0:0.2:90]*derad)); received_signal=A*s; cx=received_signal+diag(sqrt(nv/snr0/2))*(randn(K,n)+1i*randn(K,n)); Rxx=cx*cx'/n;% PM算法 G=Rxx(:,1:iwave); H=Rxx(:,iwave+1:end); P=inv(G'*G)*G'*H; % 傳播算子矩陣 Q=[P',-diag(ones(1,kelm-iwave))]; % Q矩陣for iang=1:361angle(iang)=(iang-181)/2;phim=derad*angle(iang);a=exp(-1i*twpi*d*sin(phim)).';SP(iang)=1/(a'*Q'*Q*a); end SP=abs(SP); SPmax=max(SP); SP=10*log10(SP/SPmax);% 作圖 figure h=plot(angle,SP,'-k'); set(h,'Linewidth',2) xlabel('angle(degree)') ylabel('magnitude(dB)') axis([0 60 -60 0]) set(gca,'XTick',[0:10:60]) grid on hold on legend('Propagator Method')

運行結果：

figure4 譜峰搜索傳播算子算法仿真運行結果

Reference

[1] 陣列信號處理及MATLAB實現（第2版）$\text{c?}$ 張小飛 李建峰 徐大專 等 著.

總結

以上是生活随笔為你收集整理的CSI笔记【12】：阵列信号处理及MATLAB实现（第2版）阅读随笔(四)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。