mymusic.m

很经典的关于计算ＤＯＡ的算法——music算法

clear;
clc


d=1; 
c=342;%声速　　　　　　　　　　
Fs=2000;%采样率
N=256;%帧数
RAND1=2*(rand(10000,1)-0.5);%sound sources

RAND2=2*(rand(10000,1)-0.5);



MainFreqDown=160;%%% narrow band Gauss signal ,from 160Hz to 170Hz%%%
MainFreqUp=170;
SignalFreq=(MainFreqDown+MainFreqUp)/2;%%  Center Freq       %%
NarrowBandFilter= fir1(150,[MainFreqDown*2/Fs,MainFreqUp*2/Fs]);


% Signal1=filter(NarrowBandFilter,1,RAND1);% sound sources
% Signal2=filter(NarrowBandFilter,1,RAND2);
% 

t=0:1/Fs:30;
Signal1=sin(2*pi*160*t)';
Signal2=sin(2*pi*170*t)';






tha1=45/180*pi;%两个声源的方位角
tha2=77/180*pi;


nshiftpoint11=floor(d*cos(tha1)/c*Fs);%第一个声源在当前方位下，在两阵元之间的时延点数
nshiftpoint12=floor(2*d*cos(tha1)/c*Fs);
nshiftpoint13=floor(3*d*cos(tha1)/c*Fs);

nshiftpoint21=floor(d*cos(tha2)/c*Fs);
nshiftpoint22=floor(2*d*cos(tha2)/c*Fs);
nshiftpoint23=floor(3*d*cos(tha2)/c*Fs);




f_nshiftpoint11=d*cos(tha1)/c*Fs-nshiftpoint11;
f_nshiftpoint12=2*d*cos(tha1)/c*Fs-nshiftpoint12;
f_nshiftpoint13=3*d*cos(tha1)/c*Fs-nshiftpoint13;

f_nshiftpoint21=d*cos(tha2)/c*Fs-nshiftpoint21;
f_nshiftpoint22=2*d*cos(tha2)/c*Fs-nshiftpoint22;
f_nshiftpoint23=3*d*cos(tha2)/c*Fs-nshiftpoint23;




Maxnshiftpoint=floor(4*d/c*Fs);%最大时延点数
len=length(Signal1);



x11=Signal1(Maxnshiftpoint:len-Maxnshiftpoint);%起始从Maxnshiftpoint
x12=Signal1(Maxnshiftpoint+nshiftpoint11:len-Maxnshiftpoint+nshiftpoint11);
x13=Signal1(Maxnshiftpoint+nshiftpoint12:len-Maxnshiftpoint+nshiftpoint12);
x14=Signal1(Maxnshiftpoint+nshiftpoint13:len-Maxnshiftpoint+nshiftpoint13);

x21=Signal2(Maxnshiftpoint:len-Maxnshiftpoint);%起始从Maxnshiftpoint
x22=Signal2(Maxnshiftpoint+nshiftpoint21:len-Maxnshiftpoint+nshiftpoint21);
x23=Signal2(Maxnshiftpoint+nshiftpoint22:len-Maxnshiftpoint+nshiftpoint22);
x24=Signal2(Maxnshiftpoint+nshiftpoint23:len-Maxnshiftpoint+nshiftpoint23);



x11=x11;
x12=x12*(1-f_nshiftpoint11)+Signal1(Maxnshiftpoint+nshiftpoint11+1:len-Maxnshiftpoint+nshiftpoint11+1)*f_nshiftpoint11;
x13=x13*(1-f_nshiftpoint12)+Signal1(Maxnshiftpoint+nshiftpoint12+1:len-Maxnshiftpoint+nshiftpoint12+1)*f_nshiftpoint12;
x14=x14*(1-f_nshiftpoint13)+Signal1(Maxnshiftpoint+nshiftpoint13+1:len-Maxnshiftpoint+nshiftpoint13+1)*f_nshiftpoint13;


x21=x21;
x22=x22*(1-f_nshiftpoint21)+Signal2(Maxnshiftpoint+nshiftpoint21+1:len-Maxnshiftpoint+nshiftpoint21+1)*f_nshiftpoint21;
x23=x23*(1-f_nshiftpoint22)+Signal2(Maxnshiftpoint+nshiftpoint22+1:len-Maxnshiftpoint+nshiftpoint22+1)*f_nshiftpoint22;
x24=x24*(1-f_nshiftpoint23)+Signal2(Maxnshiftpoint+nshiftpoint23+1:len-Maxnshiftpoint+nshiftpoint23+1)*f_nshiftpoint23;


x1=x11+x21;
x2=x12+x22;
x3=x13+x23;
x4=x14+x24;


x1=x1(101:356)-mean(x1(101:356));% 只算一次music结果
x2=x2(101:356)-mean(x2(101:356));
x3=x3(101:356)-mean(x3(101:356));
x4=x4(101:356)-mean(x4(101:356));



x=[x1,x2,x3,x4];
x=x-ones(N,1)*mean(x);
x=x ./(ones(N,1)*max(abs(x)));
x=x+(rand(N,4)-0.5)*0.2;%　add noise

x=x+hilbert(x);

Rxx=x'*x;
[V,D,V1]=svd(Rxx);      
EValue=diag(D);
meanValue=mean(EValue);
[K1]=find(EValue>meanValue);
[K2]=find(EValue<meanValue);
S=V(:,K1);%%%%%% Signal subspace%%%%%%%
G=V(:,K2);%%%%%% Noise  subspace%%%%%%%

%谱峰收索
for i=1:1800
    td=[0,d*cos(i*pi*0.1/180)/c,2*d*cos(i*pi*0.1/180)/c,3*d*cos(i*pi*0.1/180)/c]';
    E=exp(-j*2*pi*SignalFreq*td);
    p(i)=1/abs((E'*G*G'*E));
end



tha=0.1:0.1:180;
plot(tha,20*log10(p))