作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
在一些实际应用中需要从加性白噪声干扰的宽带有用信号中除去窄带干扰(narrowband interference，NBI)。比如，能穿透地面和植被障碍物的超宽带雷达工作在0.01到1GHz带宽上，用的是脉冲或者线性调频波形。为了得到高的分辨率，这些波形具有很高的带宽，覆盖在0.01到1GHz 范围上的至少100MHz。然而，这些频率范围被TV、FM广播电台、蜂窝电话和其它相对窄(少于1MHz)的射频信号源所引用，因而对超宽带雷达信号产生了窄带干扰。
⛄ 完整代码
close all;clear;clc
%% 参数设定与信号生成
%-宽带信号的生成
tr=0.4;tf=2;arf=2.3;N=4096;
Fs=2;Ts=1/Fs;
point=(-2):Ts:(6);
g=1./(exp(-arf*point/tr)+exp(arf*point/tf));
s=g(2:end)-g(1:end-1);% s为有用宽带信号
plot(s);axis tight;
set(gca,'XTick',1:1.8:16);
set(gca,'XTickLabel',{'-2','-1','0','1','2','3','4','5','6'});
xlabel('时间(ns)');
ylabel('幅度');
title('宽带信号');
Ks=fft(s,1024);
Ks=fftshift(Ks);
AKs=abs(Ks);
PKs=phase(Ks);
figure;
subplot(211);
plot(AKs);axis tight;
set(gca,'XTick',1:102:1024);
set(gca,'XTickLabel',{'-1','-0.8','-0.6','-0.4','-0.2','0','0.2','0.4','0.6','0.8','1'});
xlabel('频率(GHz)');
ylabel('幅度');
title('幅度');
subplot(212);
plot(PKs);axis tight;
set(gca,'XTick',1:102:1024);
set(gca,'XTickLabel',{'-1','-0.8','-0.6','-0.4','-0.2','0','0.2','0.4','0.6','0.8','1'});
xlabel('频率(GHz)');
ylabel('相位(rad)');
title('相位');
%-窄带干扰信号的生成
F=0.1*[0.6,1,1.8,2.1,3,4.8,5.2,5.7,6.1,6.4,6.7,7,7.8,9.3]';
L=length(F);
om=2*pi*F/Fs;
A=[0.5,1,1,0.5,0.1,0.3,0.5,1,1,1,0.5,0.3,1.5,0.5];
rand('seed',1954);
phi=2*pi*rand(L,1);
n=1:N;
y=A*sin(om*n+phi*ones(1,length(n)));
Ky=periodogram(y);
Ky=10*log10(Ky);
figure;
plot(Ky);axis tight;axis([1 N/2 -100 50]);
set(gca,'XTick',1:204:N/2+1);
set(gca,'XTickLabel',{'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'});
xlabel('频率(GHz)');
ylabel('功率(dB)');
title('干扰信号周期图');
%-高斯白噪声生成
v=wgn(1,N,-10);
%-原始信号生成
n0=1000;
x=zeros(1,N);
x(n0-4:n0-4+length(s)-1)=x(n0-4:n0-4+length(s)-1)+30*s; %待传输宽带信号
figure;subplot(221);plot(x(n0-99:n0+100));axis tight;
set(gca,'XTick',1:100:200);
set(gca,'XTickLabel',{'950','1000','1050'});
xlabel('时间(ns)');
ylabel('幅度');
title('宽带信号');
x=x+v; %加入高斯白噪声
subplot(222);plot(x(n0-99:n0+100));axis tight;
set(gca,'XTick',1:100:200);
set(gca,'XTickLabel',{'950','1000','1050'});
xlabel('时间(ns)');
ylabel('幅度');
title('受干扰信号');
x=x+y; %加入窄带干扰信号
subplot(223);plot(x(n0-99:n0+100));axis tight;
set(gca,'XTick',1:100:200);
set(gca,'XTickLabel',{'950','1000','1050'});
xlabel('时间(ns)');
ylabel('幅度');
title('观测信号');
%% 干扰消除
D=600;
M=100;
rela=xcorr(x,M+D);%计算相关系数
r=zeros(M);
for c=1:M %构建输入信号的自相关矩阵R
r(c,:)=rela((M+D+2-c):(2*M+D+1-c));
end
q=rela((M+2*D+1):(2*M+2*D));%计算输入与输出信号的互相关q
h=r\q';
h1=zeros(1,M);
for c=1:M
h1(c)=h(M+1-c);%后向预测器系数
end
res=conv(h,x); %预测滤波器输出响应
res1=conv(h1,x);%后向预测滤波器输出响应
re=x;
re(D+1:end)=re(D+1:end)-res(1:N-D); %NBI后的输出信号
re(1:D)=re(1:D)-res1((D+1):2*D);
subplot(224);plot(re(n0-99:n0+100));axis tight;
set(gca,'XTick',1:100:200);
set(gca,'XTickLabel',{'950','1000','1050'});
xlabel('时间(ns)');
ylabel('幅度');
title('滤波器输出信号');
%% 结果展示
Kx=periodogram(x);
Kx=10*log10(Kx);
hh=zeros(1,N);
hh(1)=1;
hh(D+1:D+length(h))=h*(-1);
Kh=periodogram(hh);
Kh=10*log10(Kh);
Kre=periodogram(re);
Kre=10*log10(Kre);
figure;
subplot(311);
plot(Kx);axis tight;axis([1 N/2 -100 50]);
set(gca,'XTick',1:204:N/2+1);
set(gca,'XTickLabel',{'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'});
xlabel('频率(GHz)');
ylabel('功率(dB)');
title('原信号周期图');
subplot(312);
plot(Kh);axis tight;axis([1 N/2 -100 50]);
set(gca,'XTick',1:204:N/2+1);
set(gca,'XTickLabel',{'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'});
xlabel('频率(GHz)');
ylabel('功率(dB)');
title('滤波器频率响应');
subplot(313);
plot(Kre);axis tight;axis([1 N/2 -100 50]);
set(gca,'XTick',1:204:N/2+1);
set(gca,'XTickLabel',{'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'});
xlabel('频率(GHz)');
ylabel('功率(dB)');
title('输出信号周期图');
figure;spectrogram(x,512,500,4096,2e9);
figure;spectrogram(re,512,500,4096,2e9);
⛄ 运行结果
⛄ 参考文献
[1]孟藏珍, 徐向东, 袁俊泉,等. 一种基于最小二乘的DS—SS系统抗窄带干扰方法[J]. 空军预警学院学报, 2006, 20(3):174-176.
