王济matlab在振动信号处理中的应用代码.docx
- 文档编号:28408940
- 上传时间:2023-07-13
- 格式:DOCX
- 页数:129
- 大小:40.13KB
王济matlab在振动信号处理中的应用代码.docx
《王济matlab在振动信号处理中的应用代码.docx》由会员分享,可在线阅读,更多相关《王济matlab在振动信号处理中的应用代码.docx(129页珍藏版)》请在冰豆网上搜索。
王济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)';
%计算趋势项的多项式待定系数向量a
a=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
%五点滑动平均法平滑处理
%%%%%%%%%%%%%%%%%%%%%%
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)';
%将x赋值给a
a=x;
%循环m次进行平滑处理计算
fork=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;
forj=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赋值给y
y=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));
end
status=fclose(fid);
程序4-2
%五点滑动平均法平滑处理
%%%%%%%%%%%%%%%%%%%%%%
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)';
%将x赋值给a
a=x;
%循环m次进行平滑处理计算
fork=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;
forj=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赋值给y
y=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));
end
status=fclose(fid);
程序4-3
%滑动平均法消除趋势项
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
closeallhidden
%%%%%%%%%%%%%%%%%%%%%%
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:
m
forj=l+1:
n-l
b(j)=mean(a(j-l:
j+l));
end
a=b;
end
%用输入信号x减去与平滑趋势项a
y=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));
end
status=fclose(fid);
程序4-4
%五点三次法平滑处理(时域和频域)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
closeallhidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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:
it
a=x(l,:
);
%循环m次进行平滑处理
fork=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;
forj=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
%绘制平滑前后的曲线图形
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:
n
ifit==1%时域信号输出时间和幅值
fprintf(fid,'%f%f\n',t(k),y(k));
else%频域信号输出频率、实部和虚部
fprintf(fid,'%f%f%f\n',f(k),y(1,k),y(2,k));
end
end
status=fclose(fid);
程序5-1
%频域低通和带通滤波
%%%%%%%%%%%%%%%%%%%%%%
clear
clc
closeallhidden
%%%%%%%%%%%%%%%%%%%%%%
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=2^nextpow2(n);
%四舍五入取整求最小截止频率对应数组元素的下标
ni=round(fmin*nfft/sf+1);
%四舍五入取整求最大截止频率对应数组元素的下标
na=round(fmax*nfft/sf+1);
%进行FFT变换,结果存于y
y=fft(x,nfft);
%建立一个长度为nfft元素全为0的向量
a=zeros(1,nfft);
%将y的正频率带通内的元素赋值给a
a(ni:
na)=y(ni:
na);
%将y的负频率带通内的元素赋值给a
a(nfft-na+1:
nfft-ni+1)=y(nfft-na+1:
nfft-ni+1);
%进行FFT逆变换,结果存于y
y=ifft(a,nfft);
%取逆变换的实部n个元素为滤波结果列向量
y=(real(y(1:
n)))';
%绘制滤波前的时程曲线图形
subplot(2,1,1);
plot(t,x);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
gridon;
%绘制滤波后的时程曲线图形
subplot(2,1,2);
plot(t,y);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
gridon;
%打开文件输出滤波后的数据
fid=fopen(fno,'w');
fork=1:
n
fprintf(fid,'%f%f\n',t(k),y(k));
end
status=fclose(fid);
程序5-2
%频域高通和带阻滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
closeallhidden
%%%%%%%%%%%%%%%%%%%%%%
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=2^nextpow2(n);
%四舍五入取整求最小截止频率对应数组元素的下标
ni=round(fmin*nfft/sf+1);
%四舍五入取整求最大截止频率对应数组元素的下标
na=round(fmax*nfft/sf+1);
%进行FFT变换,结果存于y
y=fft(x,nfft);
%建立一个长度为nfft元素全为0的向量
a=zeros(1,nfft);
%将y的正频率带阻内的元素赋值为0
y(ni:
na)=a(ni:
na);
%将y的负频率带阻内的元素赋值为0
y(nfft-na+1:
nfft-ni+1)=a(nfft-na+1:
nfft-ni+1);
%进行IFFT变换,结果存于a
a=ifft(y,nfft);
%取逆变换的实部n个元素为滤波结果列向量
y=real(a(1:
n));
%绘制滤波前的时程曲线图形
subplot(2,1,1);
plot(t,x);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
gridon;
%绘制滤波后的时程曲线图形
subplot(2,1,2);
plot(t,y);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
gridon;
%打开文件输出滤波后的数据
fid=fopen(fno,'w');
fork=1:
n
fprintf(fid,'%f%f\n',t(k),y(k));
end
status=fclose(fid);
程序5-3
%运用完全设计法的IIR滤波器滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
closeallhidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
iffun<=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)
end
rp=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);
%确定滤波方式
switchfun
case1%低通
ft='low';
case2%高通
ft='high';
case3%带通
ft='bandpass';
case4%带阻
ft='stop';
otherwise
ft='low';
end
%根据滤波器种类进行IIR滤波器设计
switchmod
%巴特沃斯滤波器
case1
[nwn]=buttord(wp/(fs/2),ws/(fs/2),rp,rs);
[ba]=butter(n,wn,ft);
%切比雪夫Ⅰ型滤波器
case2
[nwn]=cheb1ord(wp/(fs/2),ws/(fs/2),rp,rs);
[ba]=cheby1(n,rp,wn,ft);
%切比雪夫Ⅱ型滤波器
case3
[nwn]=cheb2ord(wp/(fs/2),ws/(fs/2),rp,rs);
[ba]=cheby2(n,rs,wn,ft);
%椭圆滤波器
case4
[nwn]=ellipbuttord(wp/(fs/2),ws/(fs/2),rp,rs);
[ba]=ellip(n,rp,rs,wn,ft);
end
%计算滤波器的频率响应
[hw]=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('幅值');
%在图幅上添加坐标网格
gridon;
%绘制滤波前后的信号时程图
subpl
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 王济 matlab 振动 信号 处理 中的 应用 代码
![提示](https://static.bdocx.com/images/bang_tan.gif)