%*************************************************************************
% 一个非常经典的LS滤波用于雷达干扰自治应副瓣对消的程序
% 多通道数据条件下的LS滤波
% 基于正规方程、QR分解、奇异值分解求解LS滤波系数
%*************************************************************************
clc;close all;clear;
%% 定义频率、波长和波数
f=600e6;
c=3e8;
lamda=c/f;
k=2*pi/lamda;
d=lamda/2;
%% 定义主天线和辅助天线的个数和位置
mainN=1;
mainPos=5; %主天线的位置
assistantPos=[1 3 7 9]; %辅助天线的位置
assistantN=length(assistantPos); %辅助天线的个数
dx=(assistantPos-mainPos)*d; %辅助天线至主天线的距离
%% 干扰角度、高度和数量
jamerAzimuth=[-100 5 40 65]*pi/180; %干扰的方位
jamerElevation=[5 10 50 30]*pi/180; %干扰的高度?
jamerNum=length(jamerAzimuth); %干扰的数量
%% LFM参数
tao=20e-6; %LFM信号时间20us
B=5e6; %LFM信号带宽
mu=B/tao; %LFM信号调频斜率
fs=2*B; %LFM信号的计算机采样
t=(-tao/2+1/fs:1/fs:tao/2); %离散采样时间点
PRT=150e-6; %脉冲持续时间150us
sample=fix(fs*tao); %LFM信号采样点数
totalSample=fix(fs*PRT); %总的采样点数
echo=zeros(1,totalSample); %存储回波数据
%% 定义噪声、信噪比等
noiseAmp=1;
SNR=20;
signalAmp=10^(SNR/20);
noise=1/sqrt(2)*(randn(1,totalSample)+1i*randn(1,totalSample));%高斯白噪声
signal=zeros(1,totalSample);
signal(101:101+sample-1)=signalAmp*exp(1i*pi*mu*t.^2); %在信号中指定含LFM信号的数据
Sig=signal+noise;
figure;
subplot(2,2,1);
plot(real(signal));
title('chirp signal without noise');
subplot(2,2,2);
plot(real(noise));
title('Gauss noise')
subplot(2,2,3);
plot(real(Sig));
title('Chirp signal with noise');
subplot(2,2,4);
plot(abs(fftshift(fft(Sig))));
title('FFT for Chirp signal with noise');
%% 定义干扰信号比,干扰幅度,
JSR=30;
jamerAmp=signalAmp*10^(JSR/20)/16; %20lg(N^2*JamerAmp/SignalAmp)=JSR,所以JamerAmp=SignalAmp*10^(JSR/20)/(N^2)
gaussNoise1=1/sqrt(2)*(randn(1,totalSample)+1i*randn(1,totalSample));
aa1=jamerAmp*gaussNoise1;
bb=fir1(63,0.5); %带通滤波
cc1=conv(aa1,bb); %滑动卷积
J1=cc1(64:end); %截取数据?
gaussNoise2=1/sqrt(2)*(randn(1,totalSample)+1i*randn(1,totalSample));
aa2=jamerAmp*gaussNoise2;
bb=fir1(63,0.5); %带通滤波
cc2=conv(aa2,bb); %滑动卷积
J2=cc2(64:end); %截取数据?
gaussNoise3=1/sqrt(2)*(randn(1,totalSample)+1i*randn(1,totalSample));
aa3=jamerAmp*gaussNoise3;
bb=fir1(63,0.5); %带通滤波
cc3=conv(aa3,bb); %滑动卷积
J3=cc3(64:end); %截取数据?
gaussNoise4=1/sqrt(2)*(randn(1,totalSample)+1i*randn(1,totalSample));
aa4=jamerAmp*gaussNoise1;
bb=fir1(63,0.5); %带通滤波
cc4=conv(aa4,bb); %滑动卷积
J4=cc4(64:end); %截取数据?
J=[J1;J2;J3;J4];
%% 信号模型
Y_main=signal+noise+J1+J2+J3+J4; %主天线接收到的信号
n=1/sqrt(2)*(randn(assistantN,totalSample)+1i*randn(assistantN,totalSample));
for i=1:assistantN
X_assistant(i,:)=[exp(1i*k*dx(i)*cos(jamerElevation).*sin(jamerAzimuth))]*J+n(i,:); %辅助天线接收到的信号
end
snapNum=64; %不要超过前面的101,因为要选择无目标区的距离单元。在该区域内y(n)=i(n)+v(n)作为期望信号
delayNum=1;
Y=(Y_main(delayNum+1:delayNum+snapNum));%期望信号
X=(X_assistant(:,delayNum+1:delayNum+snapNum));
%% 基于正规方程求解LS滤波系数
Rxx=(X*X'/snapNum); %注意此处X为行向量,书上为列向量
Rxy=X*Y'/snapNum;
Wopt=inv(Rxx)*(Rxy);
Y_res=Y_main-Wopt'*X_assistant; %对消输出
CR=db(mean(abs(Y_main(1:1:end)))/mean(abs(Y_res(1:1:end))));%抑制比
% 观察LS滤波前和LS滤波后的结果
figure;
subplot(2,1,1);
plot(real(Y_main));title('Received Chip signal with interference and noise');
subplot(2,1,2);
plot(real(Y_res));title('The signal after LS filtering');
%% QR分解求解LS滤波系数
[Q2,R2]=qr(X');% 对X进行QR分解
Wopt2=inv(R2(1:4,:))*Q2(:,1:4)'*Y';
Y_res2=Y_main-Wopt2'*X_assistant;
CR2=db(mean(abs(Y_main(1:1:end)))/mean(abs(Y_res2(1:1:end))));%抑制比
figure;
subplot(2,1,1);
plot(real(Y_main));title('Received Chip signal with interference and noise');
subplot(2,1,2);
plot(real(Y_res2));title('The signal after QR-LS filtering');
%% 奇异值分解求解LS滤波系数
[U,S,V]=svd(X');%对X进行奇异值分解
U1=U(:,1:4);%取出U的分块矩阵
S1=S(1:4,:);%取奇异值
V1=V;%取出V的分块矩阵
Wopt3=V1*inv(S1)*U1'*Y';%计算LS滤波系数
Y_res3=Y_main-Wopt3'*X_assistant;
CR3=db(mean(abs(Y_main(1:1:end)))/mean(abs(Y_res3(1:1:end))));%抑制比
figure;
subplot(2,1,1);
plot(real(Y_main));title('Received Chip signal with interference and noise');
subplot(2,1,2);
plot(real(Y_res3));title('The signal after SVD-LS filtering');