mutidetect.m

多目标信号检测（利用并联自适应notch滤波器分离出目标信号）

%mutidetect
% 本程序拟实现在多目标定位的情况下，对独立目标对应信号的包络检测

%并行多(双)通道notch滤波器仿真

%********************************
%x输入信号
%bei=fs/f0采样比
%step步长
%door门限
%jiao相邻正交点间隔
%k为信号前沿点
%**********************************
% function k=ctzsy(x,bei,step,door,jiao)
clc
clear all;
close all;
% close all

%%
f1=14000;
f2=14500;
f3=10000;

fs=200000;
tao=0.016;
N=round(tao*fs);

snr=200;

A1=0;%信号1幅度
A2=1;%信号2幅度
A3=0;%干扰幅度
A0=0;%噪声方差

sig1=round(cos(2*pi*f1*(0:N-1)/fs)*14000);%*32767);
sig2=round(cos(2*pi*f2*(0:N-1)/fs)*14000);%*32767);
sig3=round(cos(2*pi*f3*(0:N-1)/fs)*14000);%*32767);
%sig0=cos(2*pi*(f0-1000)*(0:N-1)/fs);
sig0=A1*sig1+A2*sig2+A3*sig3;

sig1=[zeros(1,N) sig0 zeros(1,N)];%使信号位于中间
 plot(sig1);title('原始信号，未经滤波');
% %对信号带通滤波
b=fir1(128,2*[9000,15000]/fs);   %带通滤波器 算滤波器系数
sig=filter(b,1,sig1);
sig=sig/max(sig)*max(sig1);
figure;plot(sig1,'r');
hold on;
plot(sig,'g');


%窄带噪声
fl=12000;
fh=15000;
noise0=zeros(1,length(sig));
noise1=normrnd(0,1,1,length(sig));
wn=[fl*2/fs fh*2/fs];
hfilter=fir1(256,wn);
%noi=filter(hfilter,1,noise1);
noise2=conv(noise1,hfilter);%卷积
noi=noise2(129:length(noise2)-128);
noi=noi/std(noi); %std 标准差
varn=std(noi);
A=varn*sqrt(2)*10^(snr/20);

%宽带噪声
% noi=normrnd(0,4,1,length(sig));
% varn=std(noi);
% A=varn*sqrt(2)*10^(snr/20);

signoi=sig+A0*noi;
figure,plot(signoi);title('signoi');
figure,plot(abs(fft(signoi,fs)));title('fft signoi');
% axis([0 20000 0 10000])

%%
step=0.00785;

wc1=0;
ws1=0;
wc2=0;
ws2=0;
N=length(signoi);
e=zeros(1,N);
y=zeros(1,N);
for k=1:N
    sinn1=sin(2*pi*k*f1/fs);
    coss1=cos(2*pi*k*f1/fs);
    y1(k)=wc1(k)*coss1+ws1(k)*sinn1;
    
    sinn2=sin(2*pi*k*f2/fs);
    coss2=cos(2*pi*k*f2/fs);
    y2(k)=wc2(k)*coss2+ws2(k)*sinn2;
    
    y(k)=y1(k)+y2(k);
    
    e(k)=signoi(k)-y(k);
%     ws(k-2)=ws(k-3)+step*(e(k)+0.8*e(k-1)+0.5*e(k-2)+0.3*e(k-3))*sinn;
%     wc(k-2)=wc(k-3)+step*(e(k)+0.8*e(k-1)+0.5*e(k-2)+0.3*e(k-3))*coss;
    ws1(k+1)=ws1(k)+step*e(k)*sinn1;
    wc1(k+1)=wc1(k)+step*e(k)*coss1;
    
    targ1(k)=ws1(k+1)*ws1(k+1)+wc1(k+1)*wc1(k+1);
    
    ws2(k+1)=ws2(k)+step*e(k)*sinn2;
    wc2(k+1)=wc2(k)+step*e(k)*coss2;
    
     targ2(k)=ws2(k+1)*ws2(k+1)+wc2(k+1)*wc2(k+1);
end
  %
  figure,plot(sqrt(targ1));title('targ1');grid on;
  figure,plot(sqrt(targ2));title('targ2');grid on;
% %   figure,plot(ws1);title('ws1');
% %   figure,plot(wc1);title('wc1');
% %   figure,plot(ws2);title('ws2');
% %   figure,plot(wc2);title('wc2');
  %
  
figure,plot(y1);title('y1');
figure,plot(abs(fft(y1/1000,fs)));title('fft y1');
figure,plot(y2);title('y2');
figure,plot(abs(fft(y2/1000,fs)));title('fft y2');
% figure,plot((0:6000-1)*1000/fs,y1)
% title('12.5kHz通道输出信号'),xlabel('时间/ ms'),ylabel('信号幅度/ V')
% figure,plot(abs(fft(y1/1000,fs)))
% title('输出信号频谱'),xlabel('频率/ Hz'),ylabel('频谱幅度')
% % axis([0 20000 0 7])
% figure,plot((0:6000-1)*1000/fs,y2)
% title('13.5kHz通道输出信号'),xlabel('时间/ ms'),ylabel('信号幅度/ V')
% figure,plot(abs(fft(y2/1000,fs)))
% title('输出信号频谱'),xlabel('频率/ Hz'),ylabel('频谱幅度')
% % axis([0 20000 0 10])

% figure,plot((0:6000)*1000/fs,ws1)
% title('13.5kHz滤波器权值ws'),xlabel('时间/ ms'),ylabel('ws幅度')
% figure,plot((0:6000)*1000/fs,wc1)
% title('13.5kHz滤波器权值wc'),xlabel('时间/ ms'),ylabel('wc幅度')
% figure,plot((0:6000)*1000/fs,-wc1)
% title('滤波器权值wc'),xlabel('时间/ ms'),ylabel('wc幅度')
%%
% clc
% clear all
% close all
% 
% p=0.055;
% w1=pi/4;
% w2=pi/3;
% w3=2*pi/5;
% W=[w1 w2 w3];
% 
% w=linspace(0,pi,360);
% x=w/pi;
% z=exp(j*w);
% 
% for i=1:3
%     Hc(i,:)=p*(cos(W(i)).*z-1)./(z.^2-2*cos(W(i)).*z+1);
% end
% 
% Hs=Hc(1,:)+Hc(2,:)+Hc(3,:);
% 
% for j=1:3
%     Hby(j,:)=Hc(j,:)./(1+Hs);
% end
% 
% Hbz=1./(1+Hs);
% 
% for i=1:3
%     figure,plot(x,abs(Hby(i,:)))
%     figure,plot(x,angle(Hby(i,:)))
% end
% 
% figure,plot(x,abs(Hbz))
% figure,plot(x,angle(Hbz))