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

王济 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;