物理应用基于matlab麦克风阵列近场波束形成的典型方法仿真【含Matlab源码 2196期】.zip

%%CSM算法，本质与ISM一样，只不过加了一个聚焦矩阵 %%画波束图，画不同点频率波束图，显示最低频率、中间频率、最高频率 clc;clear all;%close all %%step 1 信号源和阵列参数设置********************************************************************************* f0 = 1500; %%信号中心频率(Hz) B = 200; %%带宽(Hz) fs = 2*B; %%采样频率(Hz),带通采样定理 fl = f0-B/2; %%信号起始频率(Hz) fh = f0+B/2; %%信号终止频率(Hz) T = 1; %%信号持续时间(s) Tr = 5; %%工作周期为5s snr = 0; %%信噪比(dB) angle = 0; %%信号到达方向(度数) M = 16; %%麦克风阵列个数 c = 340; %%声速(m/s) d = 0.5*c/fh; %%阵元间距按最高频率半波长来计算，当然也可以自己设置，注意要满足空间采样定理 FFTLEN = fs; %%FFT的点数 Nm = 64; %%重复次数 N = round(T*fs); %%信号的采样点数 Nr = round(Tr*fs); %%一个周期的点数 NN = ceil(Nr*Nm/FFTLEN); %%频域快拍数 %%step 2产生接收处的源信号********************************************************************************************************** %%期望信号 t = 0:1/fs:T-1/fs; %%有期望信号的时间长度 len = length(t); s = zeros(M,len); for i = 1 : len for n = 1:M s(n,i) = 10^(snr/20)*exp(1j*(2*pi*fl*(t(i)-(n-1)*d/c*sin(angle/180*pi))+ pi * B/T *(t(i)-(n-1)*d/c*sin(angle/180*pi)).^2));%%如果相搞近场，这边得进行修改 end end s_all = zeros(M,round(Tr*fs)); %%一个工作周期分配内存 s_all(:,1:len) = s; %%将期望信号嵌进去 s_all = repmat(s_all,1,Nm); %%Nm个周期； noise_sig = randn(size(s_all))+1i*randn(size(s_all)); %%噪声 array_sig = s_all + noise_sig; %%接收处信号 clear s_all s noise_sig; %%节省内存，删除不必要的变量 %%step 3转换为FFT****************************************************************************************************************** XX = zeros(B,M,NN); %第一维频率点数，第二维是M个阵元,第三维快拍数 for n=1:NN for m=1:M if(NN==n) tmp = fft(array_sig(m,(n-1)*FFTLEN+1:end),FFTLEN); else tmp = fft(array_sig(m,(n-1)*FFTLEN+1:n*FFTLEN),FFTLEN); end XX(:,m,n) = tmp(B+1:end); end end %%step 4计算平均自相关矩阵********************************************************************************************************** R = zeros(M,M,B);%第一、二维为自相关矩阵,第三维为频率数 for i = 1 : B for j = 1 : NN tmp = XX(i,:,j); tmp = tmp(:); %%转为列向理 R(:,:,i) = R(:,:,i) + tmp * tmp'; end R(:,:,i) = R(:,:,i)/NN; end clear X; %%节省内存，删除不必要的变量 %%step 5计算聚焦矩阵和中心频率下的权向量***************************************************************************************************** fre = [fl:fh-1]; p_angle = -90:0.5:90; %%观察区间 b_pattern = zeros(B,length(p_angle)); %%波束形状 a_c = exp(-1j*2*pi*fre(B/2)*d*[0:M-1]'*sin(angle/180*pi)/c); %%中心频率点方向向量 R_c = zeros(M,M); %%中心频率自相关矩阵 for i = 1:B a = exp(-1j*2*pi*fre(i)*d*[0:M-1]'*sin(angle/180*pi)/c); %%当前频率点方向向量 E = a_c * a'; [U,S,V]=svd(E); T_f(:,:,i)=V*U'; %%转换矩阵 R_c=R_c + T_f(:,:,i)*R(:,:,i)*T_f(:,:,i)'; end R_c = R_c/B; tmp = inv(R_c) ; w_c = tmp *a_c/(a_c'*tmp*a_c); %%step 6计算所有频率下的权向量***************************************************************************************************** for i = 1:B w(:,i) = T_f(:,:,i)' * w_c; clear b_a; for m=1:M b_a(m,:) = exp(-1j*2*pi*fre(i)*d*(m-1)*sin(p_angle/180*pi)/c); % 用于方向搜索的方向向量 end b_pattern(i,:) = w(:,i)'*b_a; b_pattern(i,:) = b_pattern(i,:)/max(abs(b_pattern(i,:))); %%保存了波束图 end %%step 7画图******************************************************************************************************************************** %%b_pattern变量大小为[M,N] 其中第一维M对应的是频率fre = [fl:fh-1],间隔为1 Hz 第二维N为度数，即之前的变量p_angle figure; plot(p_angle,20*log10(abs(b_pattern(1,:)))) hold on plot(p_angle,20*log10(abs(b_pattern(B/2,:))),'r'); plot(p_angle,20*log10(abs(b_pattern(B,:))),'k'); legend(['最低频率',num2str(fre(1))],['中心频率',num2str(fre(B/2))],['最高频率',num2str(fre(B))]) xlabel('方位(\circ)') ylabel('幅度(dB)') title(['CSM\_MVDR, 方向角为',num2str(angle),'\circ']) axis([-90,90,-50,10]); grid on