单个信号、多个信号的窄带波束形成方法-MVDR实现matlab.zip

%% clear all； close all; clc; %% rad2deg=180/pi; deg2rad=pi/180; Tr = 40e-6;%线性调频信号脉宽，发射脉冲时宽20us c = 3e8; %光速 B = 5e6; %调频信号带宽 K = B/Tr; %调频率 f0 = 5*10^6; %载频 B_J = 3e6; %调频信号带宽 K_J = B_J/Tr; %调频率 f0_J = 6*10^6; %载频 fj = f0+B*0.5; %载频 Fs = 10*B; %采样率 Ts = 1/Fs; %采样时间间隔 N = round(Tr/Ts); %采样点数 % t = linspace(-Tr/2,Tr/2,N); %时间序列 t = linspace(0,Tr,N); %时间序列 d = 0.5*c/(f0+B); % 阵元间距 % 宽带信号频域处理 theta_sig = 45*deg2rad; % 指向角 theta_J = -60*deg2rad; %干扰信号角度 M =16; % 阵元数为M % 时延计算 deltaT_sig = [0:M-1]*sin(theta_sig)*d/c; deltaT_J = [0:M-1]*sin(theta_J)*d/c; % 阵列效果仿真 Sig_ary=zeros(M,N); J_ary=zeros(M,N); for m=1:M t_m = t-deltaT_sig(m); Sig_ary(m,:)=exp(1j*pi*K*t_m.^2).*exp(1j*2*pi*f0*t_m); t_m_J = t-deltaT_J(m); J_ary(m,:)=1*exp(1j*2*pi*fj*t_m_J); % 干扰信号 J_ary(m,:)=exp(1j*pi*K_J*t_m_J.^2).*exp(1j*2*pi*f0_J*t_m_J); end %% if(0) st = J_ary(5,:); fai_t = pi*K*t.^2; %线性调频信号相位 ft = K*t; %线性调频信号频率 TBP=K*Tr*Tr; st_fft=fftshift(fft(st)); f=linspace(0,Fs,length(t))-Fs/2; r=-pi*f.^2*K; figure; subplot(221); plot(f,abs(st_fft)); %线性调频信号 xlabel("频率");ylabel("信号幅值"); title("线性调频信号频谱幅值"); grid on;axis tight; %1.画网格，2.使得屏幕适合图像范围 subplot(222); plot(f,real(st_fft)); %线性调频信号 xlabel("频率");ylabel("信号实部"); title("线性调频信号频谱实部"); grid on;axis tight; %1.画网格，2.使得屏幕适合图像范围 subplot(223); plot(f,imag(st_fft)); %线性调频信号 xlabel("频率");ylabel("信号虚部"); title("线性调频信号频谱虚部"); grid on;axis tight; %1.画网格，2.使得屏幕适合图像范围 subplot(224); plot(f,r); %线性调频信号 xlabel("相位");ylabel("频率"); title("线性调频信号相位谱"); grid on;axis tight; %1.画网格，2.使得屏幕适合图像范围 end %% %-----------信号叠加--------------- FFTN = N; % FFT的点数 R = zeros(M,M); % 接收数据协方差矩阵 x = Sig_ary + 1*J_ary; % 信号叠加 % ------------时域-》频域------- spe_x = zeros(M,N); for m=1:M spe_x(m,:) = fft(x(m,:)); end % figure;plot(abs(spe_x(1,:))); abs_spe = abs(spe_x(1,:)); abs_spe=abs_spe./max(abs_spe); figure;plot(10*log10(abs_spe)); %-----------最优化求解--------------- xf= zeros(M,N); f = [0:N-1]*Fs/N; W = zeros(M,N); for m=1:M xf = spe_x(m,:);%第m个接收通道 am = exp(-1j*2*pi*f*sin(theta_sig)*m*d/c).'; Rx = xf*xf'/N;%协方差矩阵 RxR = lsqminnorm(Rx,1);%求逆 w=RxR*am*lsqminnorm(am'*RxR*am,1);%最优估计 W(m,:)=w; end %% %-----------信号逆变换------------ spe_filter = zeros(M,N); for m=1:M spe_filter(m,:)=conj(W(m,:)).*spe_x(m,:); end spe_aft_filter = sum(spe_filter); S = ifft(spe_aft_filter); abs_spe = abs(fft(S)); abs_spe=abs_spe./max(abs_spe); figure;plot(10*log10(abs_spe)); %% scan_v = 90*[-1:0.01:1]*deg2rad; fn = 350; fi = f(fn); for m=1:length(scan_v) a=exp(-1j*2*pi*fi*[0:M-1]'*sin(scan_v(m))*d/c); Wopt=W(:,fn); z(m)=Wopt'*a; end Z=20*log10(abs(z)/max(abs(z))+esp); figure;plot(scan_v*rad2deg,Z);hold on;grid on; axis([-90 90 -50 0]); plot(theta_sig*rad2deg,-30:0,'.'); plot(theta_J*rad2deg,-30:0,'.'); xlabel('\theta/o'); ylabel('Amplitude in dB'); title('宽带-频域-MVDR'); %% % v = exp(-1j*2*pi*fi*[0:M-1].'*sin(scan_v)*d/c); % w=W(:,fn); % B=(w'*v);%方向图计算 % B=20 * log10(abs(B)); % figure;plot(scan_v*rad2deg,B);grid on;title('增益方向图');xlabel('theta(°)');ylabel('Gain(dB)');