实验方法及数字信号分析处理

课程大作业

第一次作业：

题目：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