实验方法及数字信号分析处理
课程大作业
第一次作业:
题目:x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*200*t)+0.1*rand(1,m);
①画原信号的曲线;
②从上述信号中滤出5Hz信号,再画曲线;(分别用卷积和递推低通滤波器)
③从上述信号中滤出50Hz信号,再画曲线;(分别用卷积和递推带通滤波器)
Matlab源代码:
dt=0.002; %dt<1/2*fmax=0.0025,满足采样定理
t=0:dt:1; %采样时间序列
m=length(t); %求取时间采样点数
x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*200*t)+0.1*rand(1,m);
%原始信号
%figure
Subplot(321);plot(t,x); %画出原始信号图形
title('原始信号')
%卷积低通滤波器设计
F=10; %截止频率
N=50; %根据经验取N
for i=-N:N
f(1,N+1+i)=sin(2*pi*5*i*dt)/(pi*i);%根据设计公式求fi
end
f(1,N+1)=F*2*dt;
%求取f与x的卷积
M=size(f,2);NN=size(x,2); %求取矩阵列长
X=zeros(1,M+NN-1);A=rot90(f);
B=[zeros(1,M-1),x,zeros(1,M-1)];
for i=1:M+NN-1
for j=1:M
temp=A(j)*B(j+i-1);
X(i)=X(i)+temp;
end
end
%figure
Subplot(322) ;plot(X) %画出低滤波后频率为5Hz信号图
title('卷积低通滤波器')
%卷积带通滤波器设计
F2=100;F1=10; %取上下截止频率
>> for i=-N:N
ff(1,N+1+i)=2*sin(pi*(F2-F1)*i*dt)*cos(pi*(F2+F1)*i*dt)/(pi*i);
%根据带通滤波器系数公式求ff
end
ff(1,N+1)=2*(F2-F1)*dt;
Y=conv(ff,x);%直接利用卷积公式求ff与信号x的卷积
%M=size(ff,2);NN=size(x,2); %利用上述源程序求卷积
%X=zeros(1,M+NN-1);A=rot90(ff);
%B=[zeros(1,M-1),x,zeros(1,M-1)];
%for i=1:M+NN-1
%for j=1:M
%temp=A(j)*B(j+i-1);
%X(i)=X(i)+temp;
%end
%end
%figure
Subplot(323) ;plot(Y) %画出经带通滤波器滤波后的频率为50Hz信号图形
title('卷积带通滤波器')
%递推低通滤波器设计
w=tan(dt*pi*F); %固有频率计算
f0=w*w/(1+2^0.5*w+w^2); %根据递推低通滤波器系数设计公式计算各系数
f1=2*w*w/(1+2^0.5*w+w^2);
f2=w*w/(1+2^0.5*w+w^2);
g1=-2*(1-w*w)/(1+2^0.5*w+w^2);
g2=(1-2^0.5*w+w^2)/(1+2^0.5*w+w^2);
%对y1,y2初始化
y(1)=f0*x(1);
y(2)=f0*x(2)+f1*x(1)-g1*y(1);
for k=3:length(x) %循环滤波
y(k)=f0*x(k)+f1*x(k-1)+f2*x(k-2)-g1*y(k-1)-g2*y(k-2);
end
%figure
Subplot(324);
plot(t,y) %画出递推低通滤波器滤波后频率为5Hz信号图
title('递推低通滤波器')
[nn,wn]=buttord(0.001,0.2 , 0.3,50); %直接利用MatLAB中ButterWorth公式滤波
[B,A]=butter(nn,wn);
yy=filter(B,A,x);
%figure
Subplot(325);
plot(t,yy