数字信号处理_使用经典滤波器从噪声中提取信号

%用经典滤波器从噪声中提取信号
clear all;close all;clc
N=2000;T=0.001;Tp=N*T;
t=0:T:(N-1)*T;
k=0:N-1;
f=k/Tp;     %%%f=k/N*fs
fsa=1/T;fc=fsa/10;f0=fc/10;
mt=cos(2*pi*f0*t);  
ct=cos(2*pi*fc*t);  
xt=mt.*ct;          
nt=2*rand(1,N)-1;   
figure;
plot(linspace(0,fsa,N),abs(fft(xt)))
%===============================
fp=150; fs=300; 		% 输入给定指标
Rp=1;As=60;Fs=fsa;
fb=[fp,fs];m=[0,1]; 		% 计算remezord函数所需参数f,m,dev，m=[0,1]为设计高通滤波器
dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];
[n,fo,mo,W]=remezord(fb,m,dev,Fs);	% 确定remez函数所需参数
hn=remez(n,fo,mo,W);	% 调用remez函数进行设计
figure;
plot(linspace(0,2*pi,1000),abs(fft(hn,1000)))
title('高通滤波器幅频特性')
xlabel('数字角频率')
ylabel('幅度')
xlim([0,2*pi])
set(gca,'XTick',[0:pi:2*pi])
set(gca,'xtickLabel',{'0','π','2π'})
grid on
%================================
yt=filter(hn,1,10*nt);      %产生高通噪声
%============================
st=xt+yt;           %噪声加信号
fst=fft(st,N);
figure;
subplot(2,1,1)
plot(t,st);grid;xlabel('t/s');ylabel('x(t)');
axis([0,Tp/8,min(st),max(st)]);title('(a) 信号加噪声波形')
subplot(2,1,2)
plot(f,abs(fst)/max(abs(fst)));grid;title('(b) 信号加噪声的频谱')
axis([0,Fs/2,0,1.2]);
xlabel('f/Hz');ylabel('幅度')
%==============================
fp=150; fs=200; 		% 输入给定指标
Rp=0.2;As=60;Fs=fsa;
fb=[fp,fs];m=[1,0]; 		% 计算remezord函数所需参数f,m,dev，m=[1,0]为设计低通滤波器
dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];
[n,fo,mo,W]=remezord(fb,m,dev,Fs);	% 确定remez函数所需参数
hn=remez(n,fo,mo,W);	% 调用remez函数进行设计
fhw=fft(hn,N);		% 求设计出的滤波器频率特性
figure
subplot(2,1,1)
plot(f,abs(fhw));grid;title('(c) 低通滤波器幅频特性')
axis([0,Fs/2,0,max(abs(fhw))+0.5]);
xlabel('f/Hz');ylabel('幅度')
%===============================
zt=filter(hn,1,st);
subplot(2,1,2)
plot(t,zt);grid;
axis([0,Tp/8,-1,1]);xlabel('t/s');ylabel('s(t)');
title('(d) 滤除噪声后的信号波形')