1、基本原理
阵列信号处理及matlab实现——张小飞。
P36:阵列的基本概念
P41:均匀圆阵的方向矢量
P92:MUSIC算法原理与设计步骤
music算法步骤归纳:
下面所示代码实现是基于四元麦克风阵列,课本代码是基于线阵的麦克风阵列
clear all
close all
[x1,fs]=audioread('D:\EEE.wav');
[x2,fs]=audioread('D:\EEE.wav');
[x3,fs]=audioread('D:\EEE.wav');
[x4,fs]=audioread('D:\EEE.wav');
x11=x1(1:441000);
x12=x2(1:441000);
x13=x3(1:441000);
x14=x4(1:441000); %截取2s信号
X1=[x11'
x12'
x13'
x14'];
derad = pi/180; % deg -> rad
radeg = 180/pi;
twpi = 2*pi;
kelm = 4; % 阵列数量
dd = 0.058; % space
d=0:dd:(kelm-1)*dd; %
iwave = 1; % number of DOA
snr = 10; % input SNR (dB)
n = 441000; %
ll=340/2500; %2500Hz是信号频率 波长ll=c/f
Rxx=X1*X1'/n;
InvS=inv(Rxx); %%%%求逆
[EV,D]=eig(Rxx); %%%% (eig:返回矩阵的特征向量和特征值)
EVA=diag(D)';
[EVA,I]=sort(EVA);
EVA=fliplr(EVA);%(fliplr:左右翻转)
EV=fliplr(EV(:,I));
% MUSIC
for iang = 1:361
angle(iang)=(iang-181)/2;
phim=derad*angle(iang);
%四元麦克风阵列的波程时延差,根据均匀圆阵查找对应的角度分布设计信号矩阵
d=[dd*sin(phim) dd*(sin(phim)-cos(phim)) -dd*cos(phim) 0];
a=exp(-j*twpi*d/ll).'; %twpi/ll为单位向量波数,课本36页
L=iwave;
Un=EV(:,L+1:kelm); %噪声子空间
SP(iang)=1/(a'*Un*Un'*a);
end
%
SP=abs(SP);
SPmax=max(SP);
SP=10*log10(SP/SPmax);
h=plot(angle,SP);
图1为gps的伪随机码(c/a码)下的music算法角度
图2为随机码下的music算法角度
图3为语音信号下的music算法角度
区别:以上三幅图唯一的区别在于语音信号的功率谱估计更大
以上三幅图不同信号对比的代码如下
2、算法
% DOA estimation by MUSIC
% Developed by xiaofei zhang (南京航空航天大学 电子工程系 张小飞)
% EMAIL:zhangxiaofei@nuaa.edu.cn
clear all
close all
derad = pi/180; % deg -> rad 0.0175
radeg = 180/pi;
twpi = 2*pi;
kelm = 8; % 阵列数量
dd = 0.5; % space
d=0:dd:(kelm-1)*dd; % 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000
iwave = 3; % number of DOA
theta = [10 30 60]; % 角度
snr = 10; % input SNR (dB)
n = 44100;
A=exp(-j*twpi*d.'*sin(theta*derad));%%%% direction matrix,其中A为8x3矩阵
%3个不同声音来自3个方向的角度
x13= ca_code_generator(n);
[x1,fs]=audioread('C:\EEE.wav'); %3x500随机数
x11=x1';
x11=x11(1:n);
x12=randn(1,n);
S=[x11
x12
x13
]; %分别采用不用的x12,x13可得到所需对比的信号
X=A*S; %8x500信号
X1=awgn(X,snr,'measured');%将白色高斯噪声添加到信号中,按照每采样点信噪比为10的方式添加噪声
Rxx=X1*X1'/n; %求协方差矩阵,课本p88,公式4.2.5
InvS=inv(Rxx); %%%%求逆
[EV,D]=eig(Rxx); %%%% (eig:返回矩阵的特征向量和特征值)Rxx的全部特征值,构成对角阵D,Rxx的特征向量构成EV的列向量
EVA=diag(D)'; %0.2353 0.2581 0.2648 0.2739 0.3036 5.4490 6.9303 10.3758
[EVA,I]=sort(EVA); %%%从小到大排序,I表示序号
EVA=fliplr(EVA); %(fliplr:数据左右翻转),此处翻转未做后续变量调用
EV=fliplr(EV(:,I)); %左右翻转,变成从大到小排序,便于for循环中读取变量的噪声子空间
% MUSIC
for iang = 1:361
angle(iang)=(iang-181)/2;%-90至90分成360等分,每0.5一份
phim=derad*angle(iang);%derad = pi/180= 0.0175; deg -> rad。将-90-90变成弧度
a=exp(-j*twpi*d*sin(phim)).';
L=iwave;
Un=EV(:,L+1:kelm); %噪声子空间,kelm个阵元,iwave个信号
SP(iang)=1/(a'*Un*Un'*a);%P93MUSIC算法步骤归纳,计算谱函数
end
%
SP=abs(SP);%虚部基本为0,默认取得实部数据
SPmax=max(SP);
SP=10*log10(SP/SPmax);%dB = 10log (P1/P2)P2表示参考标准
h=plot(angle,SP);
set(h,'Linewidth',2)
xlabel('angle (degree)')
ylabel('magnitude (dB)')
axis([-90 90 -60 0])
set(gca, 'XTick',[-90:30:90])
grid on
以下代码所示为GPS的c/a码发生器所设计的代码
%卫星1-C/A码发生器-定义随机码产生函数
function[c] = ca_code_generator(ca_length)
%ca_length = 1000;
%C/A两个10级寄存器初始化
reg1 = ones(1,10);
reg2 = ones(1,10);
%中间寄存器初始化
reg1_t = ones(1,10);
reg2_t = ones(1,10);
%初始时刻,第一级寄存器置“1”
reg1(1) = 1;
reg2(1) = 1;
%产生ca_length个C/A码h
for i=1:ca_length
%输出第i个C/A码,每个时刻输出一个随机码
%f= mod(reg2(2)+reg2(6),2);
%c(i) = mod(reg1(10)+f,2);
c(i) = mod(reg1(10)+reg2(2)+reg2(6),2);
for j=1:9
reg1_t(11-j) = reg1(10-j);
reg2_t(11-j) = reg2(10-j);
end
%根据C/A码产生图,给第一个寄存器赋值
reg1_t(1) = mod(reg1(3)+reg1(10),2);%reg1_t(1) = mod(reg1(4)+reg1(10),2);
reg2_t(1) = mod(reg2(2)+reg2(3)+reg2(6)+reg2(8)+reg2(9)+reg2(10),2);
%将循环移位的结果保存,进行下次随机码产生
reg1 = reg1_t;
reg2 = reg2_t;
end
for k = 1:ca_length
if(c(k) == 1)
c(k) = -1;
else
c(k) = 1;
end
end