基于cyclic-MUSIC的DOA估計(jì)Matlab仿真

1、信息與通信工程學(xué)院 陣列信號(hào)處理實(shí)驗(yàn)報(bào)告 （基于cyclic-MUSIC的DOA估計(jì)Matlab仿真）學(xué) 號(hào)：XXXXXX專 業(yè)：XXXXXX學(xué)生姓名：XXX任課教師：XXX2015年X月題目：基于cyclic-MUSIC的DOA估計(jì)Matlab仿真1. 算法簡述：DOA估計(jì)算法一般使用與平穩(wěn)信號(hào)，而生活中測量的信號(hào)往往是非平穩(wěn)的。在非平穩(wěn)的信號(hào)中，有一類信號(hào)屬于周期平穩(wěn)信號(hào)，利用周期平穩(wěn)性得到的算法稱循環(huán)（共軛循環(huán)）平穩(wěn)算法。循環(huán)平穩(wěn)的循環(huán)自相關(guān)函數(shù)表示為其中為循環(huán)角頻率。共軛循環(huán)平穩(wěn)的循環(huán)自相關(guān)函數(shù)表示為 基于天線陣列協(xié)方差矩陣的特征分類DOA估計(jì)算法中，MUSIC算法具有普遍的實(shí)用性。只

2、要已知天線陣的分布形式，無論直線陣還是圓陣，不管陣元是否等間隔分布，都可以得到高分辨率的估計(jì)結(jié)果。2. 實(shí)驗(yàn)內(nèi)容與結(jié)果：實(shí)驗(yàn)先仿真分析BPSK信號(hào)的循環(huán)平穩(wěn)和共軛循環(huán)平穩(wěn)性與循環(huán)頻率和延時(shí)的的關(guān)系，然后使用12陣元均勻線陣，陣元間距為信號(hào)波長的一半，輸入信號(hào)為1個(gè)BPSK信號(hào)、設(shè)置載波頻率10MHz、采樣頻率50MHz、快拍數(shù)1000、信號(hào)方位角-4、信噪比為25dB情況下，得到DOA估計(jì)結(jié)果。圖1 BPSK 的循環(huán)平穩(wěn)性圖2 BPSK共軛循環(huán)平穩(wěn)性圖3 DOA仿真結(jié)果3. 仿真分析由仿真結(jié)果可知，BPSK信號(hào)具有循環(huán)平穩(wěn)特性，當(dāng)循環(huán)頻率為，信號(hào)具有循環(huán)平穩(wěn)性，當(dāng)循環(huán)頻率為，信號(hào)具有共軛循環(huán)平

3、穩(wěn)性。由圖一圖二可得到合適的循環(huán)頻率和延時(shí)，例如圖一，循環(huán)頻率為、延時(shí)為40Ts時(shí)其循環(huán)自相關(guān)函數(shù)值最大。實(shí)驗(yàn)最后針對(duì)cyclic-music算法的DOA估計(jì)，使用共軛循環(huán)平穩(wěn)性設(shè)計(jì)算法，并在給定實(shí)驗(yàn)條件下仿真，有效估計(jì)出了信號(hào)的來向。4. 程序1. 程序1：clear all;close all;clc;%-%線陣 波束形成%-%參數(shù)設(shè)置S_No = 1;sensor_No = 12;azimuth = 40*pi/180 ;Fs = 40e6;Fc = 10e6; %干擾功率RB = 1e6;M = 2;%二進(jìn)制wavelength = 3e8/Fc;d = wavelength/2;K

4、= 1000; %快拍數(shù)%-%dBPs = -0;SNR = 20;M_No = 40; %碼速率等于信息速率Data_No = M_No*Fs/RB;t = 1/Fs:1/Fs:Data_No/Fs;%-%功率轉(zhuǎn)換Ps_l = 10(Ps/10);Pn_l = 10(Ps-SNR)/10);%-%信號(hào)生成bit = randint(1,M_No); %產(chǎn)生信息序列bitstream = ;for i = 1:M_No if bit(i) = 1 bitstream = bitstream, ones(1, Fs/RB); else bitstream = bitstream, -ones(1

5、, Fs/RB); endendCarrier_R = cos(2*pi*Fc*t);S_R = Carrier_R .* bitstream;Carrier_I = sin(2*pi*Fc*t);S_I = Carrier_I .* bitstream;Signal_R = sqrt(Ps_l)*S_R/sqrt(S_R*S_R/length(S_R);Signal_I = sqrt(Ps_l)*S_I/sqrt(S_I*S_I/length(S_I);%生成BPSK復(fù)信號(hào)Signal(1,:) = complex(Signal_R, Signal_I);%-%-%模擬天線接收As = ze

6、ros(sensor_No,S_No);for i = 1:sensor_No for ii = 1:S_No As(i,ii) = exp(-1i*2*pi*sin(azimuth(ii)*d/wavelength*(i-1); endendx = As * Signal;%-%加噪聲noise_R = randn(sensor_No,Data_No)*sqrt(Pn_l);noise_I = randn(sensor_No,Data_No)*sqrt(Pn_l);noise = complex(noise_R, noise_I);s = x + noise;%X = Signal(1,(1

7、:K);Nm = length(X(1,:); % 產(chǎn)生共軛循環(huán)矩陣 flagg=1; Num1= 100;Num2=5;for n = 1:Num1 for Km = 0:39 Fcj=(n-Num1/2)/Num2*1e6;%選擇循環(huán)頻率 x1=X(1:Nm-Km)-mean(X(1:Nm-Km); x2=(X(1+Km:Nm)-mean(X(1+Km:Nm).*exp(-1i*2*pi*Fcj/Fs*(1:Nm-Km); if flagg=1 xx=x1*x2; %循環(huán)平穩(wěn)算法 xn=xx/(Nm-Km);elseif flagg=2 xx=x1*conj(x2); %共軛循環(huán)平算法 x

8、n=xx/(Nm-Km);endphas=exp(1i*2*pi*Fcj/Fs*Km); R(n,Km+1)=xn*phas; endend% Max=0;% for m1=1:Num1% for m2=1:50% if R(m1,m2)=Max;% m3=m1;m4=m2;% Max=R(m1,m2);% end% end% end% m3=m3-Num1/2% m4aaa=-(Num1/2-1)/Num2:1/Num2:Num1/2/Num2;bbb=1:2:80;Rr=abs(R);meshc(bbb,aaa,Rr);ylabel(循環(huán)頻率);xlabel(延時(shí));zlabel(循環(huán)頻率

9、模);2. 程序2：clear all;close all;clc;%-%線陣 music DOA%-%參數(shù)設(shè)置Signal_No = 1;Interference_No = 1;S_No = Signal_No + Interference_No;sensor_No = 8;azimuth = 0*pi/180 5*pi/180;Fs = 30e6;Fc = 10e6;F1 = 8e6; %干擾功率theta = 10 45 52 93 43*pi/180;RB = 2e6;M = 2;%二進(jìn)制wavelength = 3e8/Fc;d = wavelength/2;K = 300; %快拍

10、數(shù)%-%dBPs = -30;SNR = 10;SIR = 00;M_No = 40; %碼速率等于信息速率Data_No = M_No*Fs/RB;t = 1/Fs:1/Fs:Data_No/Fs;%-% for xunhuan=1:1% mean1 = 0;% mean2 = 0;% for SNR =1:20%-%功率轉(zhuǎn)換Ps_l = 10(Ps/10);Pi_l = 10(Ps-SIR)/10);Pn_l = 10(Ps-SNR)/10);%-%信號(hào)生成bit = randint(1,M_No); %產(chǎn)生信息序列bitstream = ;for i = 1:M_No if bit(i)

11、 = 1 bitstream = bitstream, ones(1, Fs/RB); else bitstream = bitstream, -ones(1, Fs/RB); endendCarrier_R = cos(2*pi*Fc*t);S_R = Carrier_R .* bitstream;Carrier_I = sin(2*pi*Fc*t);S_I = Carrier_I .* bitstream;Signal_R = sqrt(Ps_l)*S_R/sqrt(S_R*S_R/length(S_R);Signal_I = sqrt(Ps_l)*S_I/sqrt(S_I*S_I/len

12、gth(S_I);%生成BPSK復(fù)信號(hào)Signal(1,:) = complex(Signal_R, Signal_I);%-%產(chǎn)生干擾%I_R = sqrt(2*Pi_l)*cos(2*pi*F1*t);%I_I = sqrt(2*Pi_l)*sin(2*pi*F1*t);for i = 2:Interference_No + 1 Signal(i,:) = sqrt(2*Pi_l)*complex(sin(2*pi*F1(i-1)*t+theta(i), sin(2*pi*F1(i-1)*t+theta(i);end%-%模擬天線接收As = zeros(sensor_No,S_No);f

13、or i = 1:sensor_No for ii = 1:S_No As(i,ii) = exp(-1i*2*pi*sin(azimuth(ii)*d/wavelength*(i-1); endendx = As * Signal;%-%加噪聲noise_R = randn(sensor_No,Data_No)*sqrt(Pn_l);noise_I = randn(sensor_No,Data_No)*sqrt(Pn_l);noise = complex(noise_R, noise_I);s = x + noise;x1 = s(:,(1:K);R=x1*x1/K;V,D=eig(R);D

14、1=diag(D);Un=V(:,1:sensor_No-S_No);Gn=Un*Un;searching_doa=-90:0.1:90;%線陣的搜索范圍為-9090度 for i=1:length(searching_doa) a_theta=exp(-j*(0:sensor_No-1)*2*pi*d*sin(pi*searching_doa(i)/180)/wavelength); Pmusic(i)=1./abs(a_theta)*Gn*a_theta); endplot(searching_doa,10*log(Pmusic),r);xlabel(入射角/degree);ylabel(

15、空間譜/dB);legend(Music Spectrum);title(經(jīng)典MUSIC估計(jì));grid on;%-% 查最大值 %-% Pmax_db1=0;% for i=1:909% if Pmusic(i)Pmax_db1% Pmax_db1=Pmusic(i);% Pmax_pitching1=i;% end% end% Pmax_db2=0;% for i=910:1801% if Pmusic(i)Pmax_db2% Pmax_db2=Pmusic(i);% Pmax_pitching2=i;% end% end% err1(SNR) = (Pmax_pitching1-1)/10-90)2;% err2(SNR) = (Pmax_pitching2-1)/10-90-5)2;% end% mean1 = mean1 + err1;% mean2 =