matlab在振动信号处理中的应用代码

程序4-1%最小二乘法消除多项式趋势项%%%%%%%%%%%%%%%%%%%%%%%%clear%清除内存中所有变量和函数clc%清除工作窗口中所显示的内容closeallhidden%关闭所有隐藏的窗口%%%%%%%%%%%%%%%%%%%%%%%%%提示用键盘输入输入数据文件名fni=input('消除多项式趋势项-输入数据文件名:','s');%以只读方式打开数据文件fid=fopen(fni,'r');sf=fscanf(fid,'%f',1);%读入采样频率值m=fscanf(fid,'%d',1);%读入拟合多项式阶数fno=fscanf(fid,'%s',1);%读入输出数据文件名x=fscanf(fid,'%f',inf);%读入时程数据存成列向量%关闭数据文件status=fclose(fid);%取信号数据长度n=length(x);%建立离散时间列向量t=(0:1/sf:(n-1)/sf)';%计算趋势项的多项式待定系数向量aa=polyfit(t,x,m);%用x减去多项式系数a生成的趋势项y=x-polyval(a,t);%将分成2行1列的图形窗口的第1列设为当前绘图区域subplot(2,1,1);%绘制x对于t的时程曲线图形plot(t,x);%在图幅上添加坐标网格gridon;%将分成2行1列的图形窗口的第2列设为当前绘图区域subplot(2,1,2);%绘制y对于t的时程曲线图形plot(t,y);%在图幅上添加坐标网格gridon;%以写的方式打开文件或建立一个新文件fid=fopen(fno,'w');%进行n次循环将计算结果写到输出数据文件中fork=1:n%每行输出两个实型数据，t为时间，y为消除趋势项后的结果fprintf(fid,'%f%f\n',t(k),y(k));%循环体结束语句end%关闭数据文件status=fclose(fid);程序4-2%五点滑动平均法平滑处理%%%%%%%%%%%%%%%%%%%%%%clearclccloseallhidden%%%%%%%%%%%%%%%%%%%%%%fni=input('五点滑动平均法平滑处理-输入数据文件名:','s');fid=fopen(fni,'r');sf=fscanf(fid,'%f',1);%采样频率m=fscanf(fid,'%d',1);%平滑次数fno=fscanf(fid,'%s',1);%输出数据文件名x=fscanf(fid,'%f',inf);%输入数据存成列向量status=fclose(fid);%取信号数据长度n=length(x);%建立离散时间列向量t=(0:1/sf:(n-1)/sf)';%将x赋值给aa=x;%循环m次进行平滑处理计算fork=1:mb(1)=(3*a(1)+2*a(2)+a(3)-a(4))/5;b(2)=(4*a(1)+3*a(2)+2*a(3)+a(4))/10;forj=3:n-2b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;endb(n-1)=(a(n-3)+2*a(n-2)+3*a(n-1)+4*a(n))/10;b(n)=(-a(n-3)+a(n-2)+2*a(n-1)+3*a(n))/5;a=b;end%将a赋值给yy=a;%将分成2行1列的图形窗口的第1列设为当前绘图区域subplot(2,1,1);%绘制平滑前的时程曲线图形plot(t,x);%添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('加速度(g)');%在图幅上添加坐标网格gridon;%将分成2行1列的图形窗口的第2列设为当前绘图区域subplot(2,1,2);%绘制平滑后的时程曲线图形plot(t,y);%添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('加速度(g)');%在图幅上添加坐标网格gridon;%打开文件输出平滑后的数据fid=fopen(fno,'w');fork=1:n%每行写两个实型数据，t为时间，y为平滑后的数据fprintf(fid,'%f%f\n',t(k),y(k));endstatus=fclose(fid);程序4-2%五点滑动平均法平滑处理%%%%%%%%%%%%%%%%%%%%%%clearclccloseallhidden%%%%%%%%%%%%%%%%%%%%%%fni=input('五点滑动平均法平滑处理-输入数据文件名:','s');fid=fopen(fni,'r');sf=fscanf(fid,'%f',1);%采样频率m=fscanf(fid,'%d',1);%平滑次数fno=fscanf(fid,'%s',1);%输出数据文件名x=fscanf(fid,'%f',inf);%输入数据存成列向量status=fclose(fid);%取信号数据长度n=length(x);%建立离散时间列向量t=(0:1/sf:(n-1)/sf)';%将x赋值给aa=x;%循环m次进行平滑处理计算fork=1:mb(1)=(3*a(1)+2*a(2)+a(3)-a(4))/5;b(2)=(4*a(1)+3*a(2)+2*a(3)+a(4))/10;forj=3:n-2b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;endb(n-1)=(a(n-3)+2*a(n-2)+3*a(n-1)+4*a(n))/10;b(n)=(-a(n-3)+a(n-2)+2*a(n-1)+3*a(n))/5;a=b;end%将a赋值给yy=a;%将分成2行1列的图形窗口的第1列设为当前绘图区域subplot(2,1,1);%绘制平滑前的时程曲线图形plot(t,x);%添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('加速度(g)');%在图幅上添加坐标网格gridon;%将分成2行1列的图形窗口的第2列设为当前绘图区域subplot(2,1,2);%绘制平滑后的时程曲线图形plot(t,y);%添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('加速度(g)');%在图幅上添加坐标网格gridon;%打开文件输出平滑后的数据fid=fopen(fno,'w');fork=1:n%每行写两个实型数据，t为时间，y为平滑后的数据fprintf(fid,'%f%f\n',t(k),y(k));endstatus=fclose(fid);程序4-3%滑动平均法消除趋势项%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclccloseallhidden%%%%%%%%%%%%%%%%%%%%%%fni=input('滑动平均法消除趋势项-输入数据文件名:','s');fid=fopen(fni,'r');sf=fscanf(fid,'%f',1);%采样频率l=fscanf(fid,'%d',1);%滑动阶次m=fscanf(fid,'%d',1);%平滑次数fno=fscanf(fid,'%s',1);%输出数据文件名x=fscanf(fid,'%f',[1,inf]);%输入数据存成行向量status=fclose(fid);%取信号数据长度n=length(x);%建立离散时间列向量t=(0:1/sf:(n-1)/sf);%生成一个元素全为1的行向量b=ones(1,l);%信号两端分别向外延伸l个数据a=[b*x(1),x,b*x(n)];b=a;%按平滑次数循环进行滑动平均处理计算趋势项fork=1:mforj=l+1:n-lb(j)=mean(a(j-l:j+l));enda=b;end%用输入信号x减去与平滑趋势项ay=x(1:n)-a(l+1:n+l);%同时绘制x对于t和y对于t的时程曲线plot(t,x,':',t,y,t,a(l+1:n+l),'-.');%添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('位移mm');%在图幅上添加图例legend('输入','输出','趋势');%在图幅上添加坐标网格gridon;%以写的方式打开文件或建立一个新文件fid=fopen(fno,'w');%进行n次循环将结果写到输出数据文件中fork=1:n%每行写两个实型数据，t为时间，y为消除趋势项后的结果fprintf(fid,'%f%f\n',t(k),y(k));endstatus=fclose(fid);程序4-4%五点三次法平滑处理(时域和频域)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclccloseallhidden%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fni=input('五点三次平滑处理-输入数据文件名:','s');fid=fopen(fni,'r');sf=fscanf(fid,'%f',1);%采样频率it=fscanf(fid,'%d',1);%数据类型(1时域,2频域)m=fscanf(fid,'%d',1);%平滑次数fno=fscanf(fid,'%s',1);%输出数据文件名x=fscanf(fid,'%f',[it,inf]);%输入数据存成行向量status=fclose(fid);%取信号数据长度n=length(x(1,:));forl=1:ita=x(l,:);%循环m次进行平滑处理fork=1:mb(1)=(69*a(1)+4*(a(2)+a(4))-6*a(3)-a(5))/70;b(2)=(2*(a(1)+a(5))+27*a(2)+12*a(3)-8*a(4))/35;forj=3:n-2b(j)=(-3*(a(j-2)+a(j+2))+12*(a(j-1)+a(j+1))+17*a(j))/35;endb(n-1)=(2*(a(n)+a(n-4))+27*a(n-1)+12*a(n-2)-8*a(n-3))/35;b(n)=(69*a(n)+4*(a(n-1)+a(n-3))-6*a(n-2)-a(n-4))/70;a=b;endy(l,:)=a;end%绘制平滑前后的曲线图形ifit==1%时域信号%建立离散时间向量nn=1:2000;t=0:1/sf:(n-1)/sf;%同时绘制x对于t和y对于t的时域曲线plot(t(nn),x(nn),':',t(nn),y(nn));xlabel('时间(s)');%添加横向坐标轴的标注ylabel('幅值');%添加纵向坐标轴的标注legend('平滑前','平滑后');%在图幅上添加图例gridon;else%频域信号%建立离散频率向量nn=1:256;f=0:sf/n:(n-1)*sf/n;%同时绘制x的实部对于f和y的实部对于f的频域曲线subplot(2,1,1);plot(f(nn),x(1,nn),':',f(nn),y(1,nn));xlabel('频率(Hz)');%添加横向坐标轴的标注ylabel('实部');%添加纵向坐标轴的标注legend('平滑前','平滑后');%在图幅上添加图例gridon;%同时绘制x的虚部对于f和y的虚部对于f的频域曲线subplot(2,1,2);plot(f(nn),x(2,nn),':',f(nn),y(2,nn));xlabel('频率(Hz)');%添加横向坐标轴的标注ylabel('虚部');%添加纵向坐标轴的标注legend('平滑前','平滑后');%在图幅上添加图例gridon;end%打开文件输出平滑后的数据fid=fopen(fno,'w');fork=1:nifit==1%时域信号输出时间和幅值fprintf(fid,'%f%f\n',t(k),y(k));else%频域信号输出频率、实部和虚部fprintf(fid,'%f%f%f\n',f(k