王济matlab在振动信号处理中的应用代码.doc

程序4-1%最小二乘法消除多项式趋势项%clear % 清除内存中所有变量和函数clc % 清除工作窗口中所显示的内容close all hidden % 关闭所有隐藏的窗口%提示用键盘输入输入数据文件名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); %在图幅上添加坐标网格grid on; %将分成2行1列的图形窗口的第2列设为当前绘图区域subplot(2,1,2); %绘制y对于t的时程曲线图形 plot(t,y); %在图幅上添加坐标网格grid on; %以写的方式打开文件或建立一个新文件 fid=fopen(fno,w); %进行n次循环将计算结果写到输出数据文件中 for k=1:n %每行输出两个实型数据，t为时间，y为消除趋势项后的结果 fprintf(fid,%f %fn,t(k),y(k); %循环体结束语句 end %关闭数据文件status=fclose(fid); 程序4-2%五点滑动平均法平滑处理%clearclcclose all hidden%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赋值给a a=x; %循环m次进行平滑处理计算 for k=1:m b(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; for j=3:n-2 b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2)/5; end b(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); %在图幅上添加坐标网格grid on; %将分成2行1列的图形窗口的第2列设为当前绘图区域subplot(2,1,2);%绘制平滑后的时程曲线图形 plot(t,y); %添加横向坐标轴的标注xlabel(时间 (s); %添加纵向坐标轴的标注 ylabel(加速度(g); %在图幅上添加坐标网格grid on; %打开文件输出平滑后的数据fid=fopen(fno,w); for k=1:n %每行写两个实型数据，t为时间，y为平滑后的数据 fprintf(fid,%f %fn,t(k),y(k); end status=fclose(fid); 程序4-2%五点滑动平均法平滑处理%clearclcclose all hidden%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赋值给a a=x; %循环m次进行平滑处理计算 for k=1:m b(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; for j=3:n-2 b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2)/5; end b(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); %在图幅上添加坐标网格grid on; %将分成2行1列的图形窗口的第2列设为当前绘图区域subplot(2,1,2);%绘制平滑后的时程曲线图形 plot(t,y); %添加横向坐标轴的标注xlabel(时间 (s); %添加纵向坐标轴的标注 ylabel(加速度(g); %在图幅上添加坐标网格grid on; %打开文件输出平滑后的数据fid=fopen(fno,w); for k=1:n %每行写两个实型数据，t为时间，y为平滑后的数据 fprintf(fid,%f %fn,t(k),y(k); end status=fclose(fid); 程序4-3%滑动平均法消除趋势项%clearclcclose all hidden%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;%按平滑次数循环进行滑动平均处理计算趋势项 for k=1:m for j=l+1:n-l b(j)=mean(a(j-l:j+l); end a=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(输入,输出,趋势); %在图幅上添加坐标网格grid on; %以写的方式打开文件或建立一个新文件fid=fopen(fno,w); %进行n次循环将结果写到输出数据文件中 for k=1:n %每行写两个实型数据，t为时间，y为消除趋势项后的结果 fprintf(fid,%f %fn,t(k),y(k); end status=fclose(fid); 程序4-4%五点三次法平滑处理(时域和频域)%clearclcclose all hidden%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,:); for l=1:it a=x(l,:);%循环m次进行平滑处理 for k=1:m b(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; for j=3:n-2 b(j)=(-3*(a(j-2)+a(j+2)+12*(a(j-1)+a(j+1)+17*a(j)/35; end b(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; end y(l,:)=a;end%绘制平滑前后的曲线图形if it=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(平滑前,平滑后);%在图幅上添加图例 grid on;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(平滑前,平滑后);%在图幅上添加图例 grid on;%同时绘制x的虚部对于f和y的虚部对于f的频域曲线 subplot(2,1,2); plot(f(nn),x(2,nn),:,f(nn),y(2,nn); xlabel(频率 (Hz); %添加横向坐标轴的标注 ylabel(虚部); %添加纵向坐标轴的标注 legend(平滑前,平滑后);%在图幅上添加图例 grid on;end%打开文件输出平滑后的数据fid=fopen(fno,w); for k=1:n if it=1 %时域信号输出时间和幅值 fprintf(fid,%f %fn,t(k),y(k); else %频域信号输出频率、实部和虚部 fprintf(fid,%f %f %fn,f(k),y(1,k),y(2,k); endend status=fclose(fid); 程序5-1%频域低通和带通滤波%clearclcclose all hidden%fni=input(频域带通滤波-输入数据文件名:,s); fid = fopen(fni,r); sf = fscanf(fid,%f,1); %采样频率fmin = fscanf(fid,%f,1); %最小截止频率fmax = fscanf(fid,%f,1); %最大截止频率sx = fscanf(fid,%s,1); %读入横向坐标轴的标注 sy = fscanf(fid,%s,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); %取大于并最接近n的2的幂次方为FFT长度nfft=2nextpow2(n); %四舍五入取整求最小截止频率对应数组元素的下标ni=round(fmin*nfft/sf+1); %四舍五入取整求最大截止频率对应数组元素的下标na=round(fmax*nfft/sf+1); %进行FFT变换,结果存于yy=fft(x,nfft); %建立一个长度为nfft元素全为0的向量 a=zeros(1,nfft); %将y的正频率带通内的元素赋值给a a(ni:na)=y(ni:na); %将y的负频率带通内的元素赋值给aa(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1); %进行FFT逆变换,结果存于yy=ifft(a,nfft); %取逆变换的实部n个元素为滤波结果列向量y=(real(y(1:n); %绘制滤波前的时程曲线图形subplot(2,1,1); plot(t,x); %添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on; %绘制滤波后的时程曲线图形subplot(2,1,2); plot(t,y); %添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on;%打开文件输出滤波后的数据 fid=fopen(fno,w);for k=1:n fprintf(fid,%f %fn,t(k),y(k); end status=fclose(fid); 程序5-2%频域高通和带阻滤波%clearclcclose all hidden%fni=input(频域带阻滤波-输入数据文件名:,s); fid=fopen(fni,r); sf = fscanf(fid,%f,1); %采样频率fmin = fscanf(fid,%f,1); %最小截止频率fmax = fscanf(fid,%f,1); %最大截止频率sx = fscanf(fid,%s,1); %读入横向坐标轴的标注 sy = fscanf(fid,%s,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); %取大于并最接近n的2的幂次方为FFT长度nfft=2nextpow2(n); %四舍五入取整求最小截止频率对应数组元素的下标ni=round(fmin*nfft/sf+1); %四舍五入取整求最大截止频率对应数组元素的下标na=round(fmax*nfft/sf+1); %进行FFT变换,结果存于yy=fft(x,nfft); %建立一个长度为nfft元素全为0的向量 a=zeros(1,nfft); %将y的正频率带阻内的元素赋值为0y(ni:na)=a(ni:na); %将y的负频率带阻内的元素赋值为0y(nfft-na+1:nfft-ni+1)=a(nfft-na+1:nfft-ni+1); %进行IFFT变换,结果存于aa=ifft(y,nfft); %取逆变换的实部n个元素为滤波结果列向量y=real(a(1:n);%绘制滤波前的时程曲线图形subplot(2,1,1);plot(t,x); %添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on;%绘制滤波后的时程曲线图形subplot(2,1,2); plot(t,y);%添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on; %打开文件输出滤波后的数据 fid=fopen(fno,w); for k=1:n fprintf(fid,%f %fn,t(k),y(k); end status=fclose(fid); 程序5-3%运用完全设计法的IIR滤波器滤波%clearclcclose all hidden%fni=input(运用完全设计法的IIR滤波器滤波-输入数据文件名:,s); fid=fopen(fni,r); fs = fscanf(fid,%f,1); %采样频率%滤波方式选择（1低通,2高通,3带通,4带阻）fun = fscanf(fid,%d,1); %滤波器选择（1巴特沃斯,2切比雪夫,3切比雪夫,4椭圆）mod = fscanf(fid,%d,1); if fun=2 wp = fscanf(fid,%f,1);%通带截止频率(Hz) ws = fscanf(fid,%f,1);%阻带截止频率(Hz)else wp = fscanf(fid,%f,2);%通带截止频率(Hz) ws = fscanf(fid,%f,2);%阻带截止频率(Hz)endrp = fscanf(fid,%f,1); %通带波动系数(dB)rs = fscanf(fid,%f,1); %阻带衰减系数(dB)fno = fscanf(fid,%s,1); %输出数据文件名x = fscanf(fid,%f,1,inf);%输入数据存成行向量status=fclose(fid); %确定滤波方式switch fun case 1 %低通 ft=low; case 2 %高通 ft=high; case 3 %带通 ft=bandpass; case 4 %带阻 ft=stop; otherwise ft=low;end %根据滤波器种类进行IIR滤波器设计switch mod %巴特沃斯滤波器 case 1 n wn=buttord(wp/(fs/2),ws/(fs/2),rp,rs); b a= butter(n,wn,ft);%切比雪夫型滤波器 case 2 n wn=cheb1ord(wp/(fs/2),ws/(fs/2),rp,rs); b a= cheby1(n,rp,wn,ft);%切比雪夫型滤波器 case 3 n wn=cheb2ord(wp/(fs/2),ws/(fs/2),rp,rs); b a= cheby2(n,rs,wn,ft);%椭圆滤波器 case 4 n wn=ellipbuttord(wp/(fs/2),ws/(fs/2),rp,rs); b a= ellip(n,rp,rs,wn,ft);end %计算滤波器的频率响应h w=freqz(b,a,1024,fs); m=length(x);t=0:1/fs:(m-1)/fs;%用设计出的滤波器进行滤波y=filter(b,a,x);%绘制滤波器的频率响应图subplot(2,1,1); plot(w,abs(h);%添加横向坐标轴的标注xlabel(频率 (Hz); %添加纵向坐标轴的标注ylabel(幅值); %在图幅上添加坐标网格 grid on; %绘制滤波前后的信号时程图subplot(2,1,2); plot(t,x,:,t,y);%添加横向坐标轴的标注xlabel(时间 (s); %添加纵向坐标轴的标注 ylabel(加速度(g); %在图幅上添加图例 legend(滤波前,滤波后);grid on;%打开文件输出滤波后的数据 fid=fopen(fno,w); for k=1:m fprintf(fid,%f %f n,t(k),y(k); endstatus=fclose(fid);程序5-4%多带FIR滤波器滤波%clearclcclose all %fni=input(多带FIR滤波器滤波-输入数据文件名:,s); fid=fopen(fni,r); fs = fscanf(fid,%f,1); %采样频率%滤波器函数选择（1=fir2,2=firls,3=remez,4=cremez）mod = fscanf(fid,%d,1); n = fscanf(fid,%d,1); %滤波器阶数 k = fscanf(fid,%d,1); %频率点数 f = fscanf(fid,%f,k); %指定频率点(Hz)%通阻状态(0=阻 or 1=通)，对于函数cremez为期望幅值 a = fscanf(fid,%f,k); w = fscanf(fid,%f,k/2); %权向量fno = fscanf(fid,%s,1); %输出数据文件名%输入数据存成行向量（1加噪信号，2原始信号）x = fscanf(fid,%f,2,inf);status=fclose(fid); m=length(x(1,:);t=0:1/fs:(m-1)/fs;%频率归一化处理f=2*f/fs;%多带FIR滤波器设计switch mod case 1 %fir2滤波函数 b= fir2(n,f,a); case 2 %firls滤波函数 b= firls(n,f,a,w); case 3 %remez滤波函数 b= remez(n,f,a,w); case 4 %cremez滤波函数 b= cremez(n,f,a,w,real);end %用多带FIR滤波器进行滤波v=filter(b,1,x(1,:);%消除延时处理y=zeros(1,m); y(1:m-n/2)=v(n/2+1:m);%计算滤波器的频率响应h w=freqz(b,1,1024,fs); %绘制滤波器的频率响应图subplot(2,1,1)plot(w,abs(h);%添加横向坐标轴的标注xlabel(频率 (Hz); %添加纵向坐标轴的标注 ylabel(幅值); grid on; %绘制滤波前后的信号时程比较图subplot(2,1,2)nn=1:200; plot(t(nn),x(1,nn),-.,t(nn),x(2,nn),:,t(nn),y(nn);%添加横向坐标轴的标注xlabel(时间 (s); %添加纵向坐标轴的标注 ylabel(幅值); %在图幅上添加图例legend(加噪信号,原始信号,滤波信号); grid on; %打开文件输出滤波后的数据fid=fopen(fno,w); for k=1:m fprintf(fid,%f %f n,t(k),y(k); endstatus=fclose(fid);程序5-5%频域积分%clearclcclose all hidden%fni=input(频域积分-输入数据文件名:,s); fid=fopen(fni,r); sf = fscanf(fid,%f,1); %采样频率fmin = fscanf(fid,%f,1); %最小截止频率fmax = fscanf(fid,%f,1); %最大截止频率c = fscanf(fid,%f,1); %单位变换系数it = fscanf(fid,%f,1); %积分次数sx = fscanf(fid,%s,1); %横向坐标轴的标注 sy1 = fscanf(fid,%s,1); %纵向坐标轴输入单位的标注 sy2 = fscanf(fid,%s,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; %大于并最接近n的2的幂次方为FFT长度 nfft=2nextpow2(n); %FFT变换y=fft(x,nfft); %计算频率间隔（Hz/s）df=sf/nfft;%计算指定频带对应频率数组的下标ni=round(fmin/df+1); na=round(fmax/df+1);%计算圆频率间隔（rad/s） dw=2*pi*df; %建立正的离散圆频率向量w1=0:dw:2*pi*(0.5*sf-df); %建立负的离散圆频率向量w2=2*pi*(0.5*sf-df):-dw:0; %将正负圆频率向量组合成一个向量w=w1,w2; %以积分次数为指数，建立圆频率变量向量w=w.it; %进行积分的频域变换a=zeros(1,nfft);a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1); if it=2 %进行二次积分的相位变换 y=-a; else %进行一次积分的相位变换 real(y)=imag(a); imag(y)=-real(a); enda=zeros(1,nfft);%消除指定正频带外的频率成分 a(ni:na)=y(ni:na); %消除指定负频带外的频率成分 a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1); %IFFT变换 y=ifft(a,nfft); %取逆变换的实部n个元素并乘以单位变换系数为积分结果y=real(y(1:n)*c; %绘制积分前的时程曲线图形 subplot(2,1,1); plot(t,x); %添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注 ylabel(sy1); grid on; %绘制积分后的时程曲线图形 subplot(2,1,2); plot(t,y);%添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注 ylabel(sy2); grid on; %打开文件输出积分后的数据 fid=fopen(fno,w); for k=1:n fprintf(fid,%f %fn,t(k),y(k); endstatus=fclose(fid);程序5-6%频域微分%clearclcclose all hiddenformat long%fni=input(频域微分-输入数据文件名:,s); fid=fopen(fni,r); sf = fscanf(fid,%f,1); %采样频率fmin = fscanf(fid,%f,1); %最小截止频率fmax = fscanf(fid,%f,1); %最大截止频率c = fscanf(fid,%f,1); %单位变换系数it = fscanf(fid,%f,1); %微分次数sx = fscanf(fid,%s,1); %横向坐标轴的标注 sy1 = fscanf(fid,%s,1); %纵向坐标轴输入单位的标注 sy2 = fscanf(fid,%s,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; %大于并最接近n的2的幂次方为FFT长度 nfft=2nextpow2(n); %FFT变换y=fft(x,nfft); %计算频率间隔（Hz/s）df=sf/nfft;%计算指定频带对应频率数组的下标ni=round(fmin/df+1); na=round(fmax/df+1);%计算圆频率间隔（rad/s） dw=2*pi*df; %建立从负到正的圆频率向量w1=0:dw:2*pi*(0.5*sf-df); w2=2*pi*(0.5*sf-df):-dw:0; w=w1,w2; %以微分次数为指数，建立圆频率变量向量 w=w.it; %进行微分的频域变换a=y.*w; if it=2 %进行二次微分的相位变换 y=-a; else %进行一次微分的相位变换 real(y)=imag(a); imag(y)=-real(a); enda=zero