王济 matlab在振动信号处理中地应用代码.docx
《王济 matlab在振动信号处理中地应用代码.docx》由会员分享，可在线阅读，更多相关《王济 matlab在振动信号处理中地应用代码.docx（112页珍藏版）》请在冰豆网上搜索。
王济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;
升级会员 