数字信号处理实验九

三、实验内容
1.对周期方波信号进行滤波。
（1）生成一个基频为10Hz的周期方波信号x(t)。
（2）设计一数字滤波器，滤去该周期信号中40Hz以外的频率成分，观察滤波前后信号时域波形及频谱。
（3）若该信号x(t)淹没在噪声中（随机噪声用randn((1,N)生成），试用filter函数滤去噪声。
滤波前的时域波形
ylabel('幅度')
HC1=filter(b,a,HC);
figure
plot(f,abs(HC1));
title('滤除噪声之后')
xlabel('频率f')
ylabel('幅度')
2.若原始落信号由5Hz、15Hz、30Hz三个幅度相等的正弦信号构成。分别设计一个FIR和IIR数字滤波器滤除5Hz和30Hz频率成分。
X=fft(x0);
figure
plot(abs(X));
title('原始信号频谱');
[N,Wn] = buttord([0.11,0.15],[0.05,0.4],1,40)
[b,a] = butter(N,Wn);
x1=filter(b,a,x);
figure
plot(x1);axis([0,200,-0.5,0.5])
[p,q] = butter(M,Rp,Wn,'stop');
figure
freqz(p,q,256,1000)
z=filter(p,q,y);
figure
stem(z)
figure
plot(f,abs(fftshift(fft(z,512))));
低通滤波器频谱
通过低通滤波器的心电图频谱
带阻滤波器频谱
x0=cos(2*pi*5*k*0.01)+cos(2*pi*15*k*0.01)+cos(2*pi*30*k*0.01);
x=cos(2*pi*5*t)+cos(2*pi*15*t)+cos(2*pi*30*t);
figure
plot(x);axis([0,1000,-3,3])
title('原始信号时域波形');
滤波后的频谱
N=512;
f=(-N/2)*Fs/N:Fs/N:(N/2-1)*Fs/N;
fy=fftshift(fft(y,N))/N;
figure
plot(f,abs(fy))
title('滤波后信号的频谱')
xlabel('频率')
ylabel('函数值fx')
figure
freqz(b,a);
提示：心电信号通常均分布在200Hz范围内。
在实验中，以x(n)作为输入序列，滤除其中的干扰成分。{x(n)}={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0}
title('滤波后时域波形');
figure
freqz(b,a,256,1000)
m=freqz(b,a,256,1000);
title('iir带通滤波器');
%axis([0,500,-80,0]);
grid on
y=X.*m';
figureplot(abs(y));
title('滤波后频谱');
思考题
直接运行程序，结果输出滤波器幅频特性曲线图，把有噪声的心电图采集信号波形图和经过三级二阶滤波器滤波后的心电图信号波形图进行对比，总结滤波效果。
x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,...
0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,...
fx=fftshift(fft(x,N))/N;
figure
plot(f,abs(fx));
title('原心电图的频谱')
xlabel('频率')
ylabel('函数值fx')
[N,Wc]=buttord(150/500,250/500,3,60);
[b,a]=butter(N,Wc);
figure
Rs=0.01;
f=[0.05,0.11,0.15,0.4];
a=[0,1,0];
dev=Rs*ones(1,length(a));
[M,Wc,beta,ftype]=kaiserord(f,a,dev);
h=fir1(M,Wc,ftype,kaiser(M+1,beta));
omega=linspace(0,pi,256)
figure
stem(x)
X=fft(x);
figure
plot(abs(X))
FIR滤波器
k=0:255;
x=cos(2*pi*5*k*0.01)+cos(2*pi*15*k*0.01)+cos(2*pi*30*k*0.01);
X=fft(x);
plot(abs(X));title('原始信号频谱');
4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
figuபைடு நூலகம்e
stem(x);
title('原心电图的时域序列')
Fs=1000; %采样频率
T=1/10; %采样长度
N=T*Fs; %采样点数
f=(-N/2)*Fs/N:Fs/N:((N/2-1))*Fs/N;
title('低通滤波器的频率特性')
滤除噪声
noise=randn(1,100);
HC=x1+noise;
f=(-N/2)*Fs/N:Fs/N:(N/2-1)*Fs/N;
HC=fftshift(fft(HC,N))/N;
figure
plot(f,abs(HC));
title('滤除噪声之前')
xlabel('频率f')
freqz(b,a,256,1000)
for i=1:3
y=filter(b,a,x);
end
figure
plot(f,abs(fftshift(fft(y,512))));
Wp =[45 55]/500; Ws =[45-5 55+5]/500;
Rp = 3; Rs = 40;
[M,Wn] = buttord(Wp,Ws,Rp,Rs)
mag=freqz(h,[1],omega);
figure
plot(omega/pi,20*log10(abs(mag)));title('fir带通滤波器');
Y=X.*mag;
figure
plot(abs(Y));title('滤波后频谱');
%IIR
k=0:255;
t=0:.001:2.55;
（1）绘制原始信号时域波形和幅度频谱。
（2）分别设计FIR和IIR数字滤波器，滤绘制其幅频特性。
（3）利用设计的滤波器对信号进行滤波，绘制输出信号的时域波形和幅度频谱。
时域波形
（1）k=0:255;
x=cos(2*pi*5*k*0.01)+cos(2*pi*15*k*0.01)+cos(2*pi*30*k*0.01);
二、实验原理
首先对待滤波的信号进行频谱分析，观察信号频率分布的规律，从而确定数字滤波器的类型（FIR滤波器、IIR滤波器、自适应滤波器、小波滤波器等）。
在加性噪声的情况下，若信号的频谱与噪声的频谱基本不重叠，可以采用频率选择滤波器（FIR滤波器、IIR滤波器）。
若信号的频谱与噪声的频谱重叠较多，可以采用自适应滤波、小波滤波等。
2.如何根据含有噪声信号的频谱特性选择滤波器的类型和设计指标？
答：根据采集到的信号获得频谱图，由时域频域的对应关系，确定需要滤除噪声的特性，然后从以下角度确定所需要的滤波器：（1）频响特性角度：IIR滤波器设计时不考虑相位特性，且通常相位都是非线性的，而FIR滤波器在满足幅频特性要求的同时，还能获得比较严格的线性相位特性，利用窗函数或者其他算法可以逼近更加任意的频响特性，因此性能优越，使用范围更广；（2）稳定性问题：IIR滤波器设计时，极点必须在单位圆之内；而FIR滤波器极点在单位圆内，因此始终稳定；（3）滤波器结构的影响：IIR滤波器一般采用递归结构，存在有输出对输入的反，而IIR滤波器阶次相对较低，运算次数少，存储单元少，FIR滤波器正好相反；（4）设计工作量：FIR无表可查，需要用到迭代法，计算量较大；而IIR滤波器相对简单，有现成的计公式和数据表格可用。
若为乘性噪声，可以根据同态滤波的原理对信号进行预处理，然后再按照加性噪声的情况处理。
在确定了数字滤波器类型后，还需要根据信号时域特性、频域特性、或时频特性确定滤波器的设计参数，设计出相应的数字滤波器。
最后，利用该数字滤波器对信号进行滤波，在时域和频域观察信号滤波的主观及客观效果。若主观及客观效果满足要求，说明分析过程和滤波方法正确有效，若不满足要求，需要重新分析和设计。
x1=square(2*pi*10*t1);
f=(-N/2)*Fs/N:Fs/N:((N/2-1))*Fs/N;
fx=fftshift(fft(x,N))/N;