7-9.m

来自「经典数字信号处理滤波器的源代码 重点是利用巴特沃斯模拟滤波器转而设计其它数字滤」· M 代码 · 共 53 行

%例程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);