zaboduixiao3.m

这是无源雷达中如何消除直达波的源程序

clc;
clear;
fc=10.7e6;%10.7e6;%载频
fs=0.408e6*2;
df1=1000;%多普勒
df2=2000;%多普勒
df3=3000;%多普勒
m=14.5;td=m/fs; %时延6
N=8192*2;%fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%白噪声过滤低通波器 作为调制信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1 = randn(1,N); %wgn(1,fs,0);

ffty1=fft(y1);
cofft=ffty1.*conj(ffty1);% plot(abs(cofft));

f2=3.5e5;
% fs=0.408e6*2;
b = fir1(127,f2/fs);
 scope=(1:length(b))*fs/length(b);
% figure(gcf+1);
figure(1)
plot(scope,20*log10(abs(fft(b))));
modsig=conv(y1,b);
% plot(abs(fft(modsig)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NN=N/8;%4*1024;%512*2*2*2%0.25s  data length
t=(0:1/fs:(N-1)/fs);%采样时间序列
deviation=4;fm=5000;
% ref=cos(2*pi*fc*t+deviation*sin(2*pi*fm*t));% 参考信号
ref=cos(2*pi*fc*t+deviation*cumsum(modsig(1:length(t))));% 参考信号

ref=conv(ref,b);  %%%%%%%%%%%%%%%%%%%%%
% scope=(-2:0.1:2)  %观察幅值的分布状况
% figure(gcf+1);    %
% hist(ref,scope);  % 察看其柱状图
% grid on;

% scope=(1:NN)*fs/NN;
% figure(1);plot(scope,20*log10(abs(fft(ref(1:NN)))));

% taget=cos(2*pi*(fc+df)*(t-td)+deviation*sin(2*pi*fm*(t-td)));% 目标信号
taget1=0.0001*cos(2*pi*(fc+df1)*(t)+deviation*cumsum(modsig(1:length(t))));% 目标信号
taget2=0.0001*cos(2*pi*(fc+df2)*(t)+deviation*cumsum(modsig(1:length(t))));% 目标信号
taget3=0.0001*cos(2*pi*(fc+df3)*(t)+deviation*cumsum(modsig(1:length(t))));% 目标信号
taget1=conv(taget1,b);
n1=10;n2=20;n3=30;
delay1=[zeros(1,n1) 1];taget1=conv(taget1,delay1);
delay2=[zeros(1,n2) 1];taget2=conv(taget2,delay2);
delay3=[zeros(1,n3) 1];taget3=conv(taget3,delay3);
taget=taget1(1:N);%   +taget2(1:N)+taget3(1:N);

% snr=200;
% taget=awgn(taget,snr,'measured') ;

ref=ref(1:N);
sig=taget(1:N);
% load ref.dat;
% load sig.dat;
len=N;
L=100;
w=zeros(1,L);
% sig=cos(2*pi*20*t);
% ref=sig+0.2*randn(size(sig));
% 
% tem=sig;
% sig=ref;
% ref=tem;

%figure(gcf+1)
sig=sig+circshift(ref,[0 10])+circshift(ref,[0 200])+0.2*circshift(ref,[0 300])+0.05*circshift(ref,[0 400])+0.005*circshift(ref,[0 100])+0.001*circshift(ref,[0 200])++0.0001*circshift(ref,[0 500])+0.00001*circshift(ref,[0 1000]);
figure(2)
subplot(311);plot(sig(100:end));title('回波信号');
subplot(312);plot(ref(100:end));title('参考信号');


% u=1/abs(ref*ref');
u=0.0001
for k=1:(len-L+1)
    outp(k+L-1)=w*ref(k:k+L-1).';
    err(k+L-1)=sig(k+L-1)-outp(k+L-1);
    w=w+2*u*err(k+L-1)*conj(ref(k:k+L-1));     
end;
% % 变步长法 
% beita=0.999
% gama=0.01
%  MUmax=(1/(ref*ref'/N*L))*0.4;
%  mu=MUmax;
%     for k=1:(len-L+1)
%     outp(k+L-1)=w*ref(k:k+L-1).';
%     err(k+L-1)=sig(k+L-1)-outp(k+L-1);
%     w=w+2*mu*err(k+L-1)*conj(ref(k:k+L-1));
%     mu=beita*mu+(1-beita)*gama*MUmax;
%    end;


tem=sig';
tem1=ref'; 
subplot(313);plot(err(100:end));title('自适应对消后的信号');%(fftshift(abs(fft(err(100:end)))));
errsig=err(100:end)+ref(100:end);
%subplot(414);plot(errsig);
figure(5);
scope=(1:length(sig(1000:end)))*fs/length(sig(1000:end))-fs/2;
subplot(121);plot(scope,20*log10(fftshift(abs(fft(sig(1000:end))))));grid on;title('回波信号频谱');
scope=(1:length(err(8000:end)))*fs/length(err(8000:end))-fs/2;
subplot(122);plot(scope,20*log10(fftshift(abs(fft(err(8000:end))))));grid on;title('自适应对消后信号频谱');

scope=(1:length(taget(1000:end)))*fs/length(taget(1000:end))-fs/2;
figure(6);plot(scope,taget(1000:end));grid on;title('理想回波信号');

scope=(1:length(taget(1000:end)))*fs/length(taget(1000:end))-fs/2;
figure(7);plot(scope,20*log10(fftshift(abs(fft(taget(1000:end))))));grid on;title('理想回波信号频谱');