%例程7－9  自適應頻域窄帶干擾抑制
%生成調制信號，并加入加性噪聲和窄帶干擾

k=0:1023;
x1=20*sin(70*2*pi*k/40);
x2=20*cos(70*2*pi*k/40);
x=x1.*data+x2.*datb+randn(1,1024)+15*sin(25*2*pi*k/40)+20*sin(35*2*pi*k/40); 
 
w=blackman(1024);                  %通過blackman窗，減小頻譜泄漏 
Y=fft(w'.*x,1024);                   %對截取的數據進行快速傅立葉變換  R=abs(real(Y)); 
I=abs(imag(Y)); 
Ampl=abs(Y);                       %計算頻譜幅度
subplot(2,1,1);
plot(Ampl);
axis([0 1000 0 5000]);

u=0;
for i=1:1024
    u=u+Ampl(i);
end

mean1=u./1024;                     %計算頻譜幅度的均值，便于分析頻譜特性
Ampldb=10.*log10(Ampl);            %為防止溢出將幅度取分貝
u=0;
m=0;
for i=1:1024
    u=u+Ampldb(i);
end
mean=u./1024;                     %取分貝后的幅度的均值
for i=1:1024
    m=m+(Ampldb(i)).^2;
end

div=sqrt(m./1024-mean.^2);          %計算標準差
M=0.8;                           %可以通過改變此值的大小來改變門限值
Th=mean+M.*div;                  %計算出噪聲門限
a=0.06;                           %設置衰減系數
for i=1:1024
    if Ampldb(i)>Th
       X(i)=a.*Y(i);               %系數衰減法
        %X(i)=0;                 %門限置零法
        %X(i)=mean1;             %置均值法
    else
        X(i)=Y(i);                %未超過門限，保持不變
    end
end

Amplout=abs(X);
out=ifft(X,1024);                  %傅立葉反變換，恢復時域信號
subplot(2,1,2);
plot(Amplout);