clc;
clear;
close all;
SNR=10;%输入信号的信噪比
f0=9000;
f1=10000;%单频信号的频率
f2=11000;
T=0.05;
B=1/T;%脉宽、带宽
fs=50000;
Ts=1/fs;%采样频率
fa=8000;fb=12000;
t = 0:1/fs:(T-1/fs);
si=cos(2*pi*f0*t)+cos(2*pi*f1*t)+cos(2*pi*f2*t);
si1=[zeros(1,2500) si zeros(1,2500)];
N=1:1:length(si1)
D=256;
sigma=sqrt(var(si)/(10^(SNR/10)));
noise0=wgn(1,numel(si1)+D,0);
h0=fir1(512,[fa/(fs/2),fb/(fs/2)]);%生成fa~fb范围的带通滤波器
noise1=filter(h0,1,noise0);
noise2=noise1(D+1:numel(si1)+D)*sigma/std(noise1);
noise3=noise2*sigma/std(noise2);
outD=noise3+si1;
x0c=cos(2*pi*f0*N/fs);
x0s=sin(2*pi*f0*N/fs);
x1c=cos(2*pi*f1*N/fs);
x1s=sin(2*pi*f1*N/fs);
x2c=cos(2*pi*f2*N/fs);
x2s=sin(2*pi*f2*N/fs);
i=2*pi*B/fs;
w0c=zeros(1,length(si1)+1);
w0s=zeros(1,length(si1)+1);
w1c=zeros(1,length(si1)+1);
w1s=zeros(1,length(si1)+1);
w2c=zeros(1,length(si1)+1);
w2s=zeros(1,length(si1)+1);
e=zeros(1,length(si1));
y=zeros(1,length(si1));
y0=zeros(1,length(si1));
y1=zeros(1,length(si1));
y2=zeros(1,length(si1));
for k=1:length(si1)
y0(k)=w0c(k)*x0c(k)+w0s(k)*x0s(k);
y1(k)=w1c(k)*x1c(k)+w1s(k)*x1s(k);
y2(k)=w2c(k)*x2c(k)+w2s(k)*x2s(k);
y(k)=y1(k)+y0(k)+y2(k)
e(k)=outD(k)-y(k);
w0c(k+1)=w0c(k)+2*i*e(k)*x0c(k);
w0s(k+1)=w0s(k)+2*i*e(k)*x0s(k);
w1c(k+1)=w1c(k)+2*i*e(k)*x1c(k);
w1s(k+1)=w1s(k)+2*i*e(k)*x1s(k);
w2c(k+1)=w2c(k)+2*i*e(k)*x2c(k);
w2s(k+1)=w2s(k)+2*i*e(k)*x2s(k);
end
figure(1),plot(N/fs*1000,outD);
xlabel('时间(t)');
ylabel('幅度');
title('输入三个频率CW信号+带限白噪声');
figure(2),plot(N/fs*1000,y);
xlabel('时间(t)');
ylabel('幅度');
title('Notch滤波器输出三个频率CW信号');
figure(3),plot(N/fs*1000,y1+y2);
xlabel('时间(t)');
ylabel('幅度');
title('Notch滤波器输出需要的两个频率CW信号');
outD_fft=abs(fft(outD,fs));
F=0:length(outD_fft)-1;
figure(4),plot(F,outD_fft);
xlabel('频率(f)');
ylabel('幅度');
title('输入CW信号+带限白噪声的幅频特性曲线');
figure(5),plot(F,abs(fft(y,fs)),'b');
xlabel('频率(f)');
ylabel('幅度');
title('Notch滤波器输出三个频率CW信号的幅频特性曲线');
figure(6),plot(F,abs(fft(y1+y2,fs)),'b');
xlabel('频率(f)');
ylabel('幅度');
title('Notch滤波器输出需要的两个频率CW信号的幅频特性曲线');