《实验方法信号处理Matlab源代码》由会员分享，可在线阅读，更多相关《实验方法信号处理Matlab源代码（12页珍藏版）》请在金锄头文库上搜索。
1、文档供参考，可复制、编制，期待您的好评与关注！ 实验方法及数字信号分析处理 课程大作业 / 第一次作业：题目：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 for i=-N:Nff(1,N+1+i)=2*sin(pi*(F2-F1)*i*dt)*cos(pi*(F2+F1)*i*dt)/(pi*i);%根据带通滤波器系数
2、公式求ffend 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 %figureSubplot(323) ;plot(Y) %画出经带通滤波器滤波后的频率为50Hz信号图形title(卷积带通滤
3、波器) %递推低通滤波器设计w=tan(dt*pi*F); %固有频率计算 f0=w*w/(1+20.5*w+w2); %根据递推低通滤波器系数设计公式计算各系数 f1=2*w*w/(1+20.5*w+w2); f2=w*w/(1+20.5*w+w2); g1=-2*(1-w*w)/(1+20.5*w+w2); g2=(1-20.5*w+w2)/(1+20.5*w+w2); %对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)
4、-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);%figureSubplot(325);plot(t,yy)title(ButterWorth公式滤波) %递推带通滤波器设计 %固有频率计算 W=tan(dt*pi*(F2-F1); B=cos(dt*pi*(F2
5、+F1)/cos(dt*pi*(F2-F1);%根据递推带通滤波器系数设计公式计算各系数ff0=W*W/(1+20.5*W+W2);ff2=-2*W*W/(1+20.5*W+W2);ff4=W*W/(1+20.5*W+W2);gg1=-(4*B+2(3/2)*B*W)/(1+20.5*W+W2); gg2=(4*B2-2*W2+2)/(1+20.5*W+W2); gg3=(2(3/2)*W*B-4*B)/(1+20.5*W+W2); gg4=(1-20.5*W+W2)/(1+20.5*W+W2); dk(1)=ff0*x(1); %对dk1，dk2，dk3，dk4初始化 dk(2)=ff0*x
6、(2)-gg1*dk(1); dk(3)=ff0*x(3)+ff2*x(1)-gg1*dk(2)-gg2*dk(1); dk(4)=ff0*x(4)+ff2*x(2)-gg1*dk(3)-gg2*dk(2)-gg3*dk(1); for k=5:length(x) %循环滤波dk(k)=ff0*x(k)+ff2*x(k-2)+ff4*x(k-4)-gg1*dk(k-1)-gg2*dk(k-2)-gg3*dk(k-3)-gg4*dk(k-4); end %figure Subplot(326);plot(t,dk) %画出递推带通滤波器滤波后频率为50Hz信号图title(递推带通滤波器)第二次
7、作业：题目：根据已知位移曲线，求速度曲线要求：由数据文件画出位移曲线（ t=0.0005s ）；对位移数据不作处理，算出速度并画出速度曲线；对位移数据进行处理，画出位移曲线，并与原位移曲线对比；画出由处理后的位移数据算出的速度曲线；写出相应的处理过程及分析。Matlab源程序：clear x=load(C:UsersAdministratorDesktopbiancendat2.dat); %加载原始数据x=x; %矩阵倒置 dt=0.0005; %采样时间L=length(x); %计算x的长度t=0:dt:(L-1)*dt;plot(t,x) %画出原始数据图形title(原始位移曲线)x
8、label(时间t/s);ylabel(位移x(t); %求取未经数据处理速度图形 v=zeros(1,L);for i=1:L-1 %利用for循环求未经处理位移的速度 v(i)=(x(i+1)-x(i)/dt;end figure plot(t,v); %画出未经数据处理速度图形 title(数据未处理的速度曲线)xlabel(t/s);ylabel(v(t); %频谱分析X=fft(x); %对原始数据傅里叶变换%PXX=abs(X); PXX=X.*conj(X)/L; %求取功率谱幅值f=(1/dt)*(0:L/2-1)/L;figure(3); plot(f,PXX(1:L/2)
9、%画出功率谱图形% stem(f,PXX(1:L/2);title(频谱分析) %利用第一次作业的低通ButterWorth滤波器对原始信号滤波F=20; %根据上述功率谱图形，取截止频率20Hzw=tan(dt*pi*F); f0=w*w/(1+20.5*w+w2); f1=2*w*w/(1+20.5*w+w2); f2=w*w/(1+20.5*w+w2); g1=-2*(1-w*w)/(1+20.5*w+w2); g2=(1-20.5*w+w2)/(1+20.5*w+w2); y(1,1)=f0*x(1,1);y(1,2)=f0*x(1,2)+f1*x(1,1)-g1*y(1,1); fo
10、r k=3:length(x)y(1,k)=f0*x(1,k)+f1*x(1,k-1)+f2*x(1,k-2)-g1*y(1,k-1)-g2*y(1,k-2);end figure plot(t,y) %画出滤波后位移图形 title(数据处理后的速度曲线)xlabel(t/s);ylabel(y(t); %求取滤波后速度图形vv=zeros(1,L);for i=1:L-1 vv(i)=(y(i+1)-y(i)/dt;end figure; plot(t,vv); %画出滤波后速度图形 title(数据处理后的速度曲线)xlabel(t/s);ylabel(vv(t); 原始数据图原始数据局部放大图原始数据速度图原始数据速度局部放大图频率分析图处理后的位移图处理后的位移局部放大图处理后速度图处理后速度局部放大图
