第1页 / 共79页
第2页 / 共79页
第3页 / 共79页
第4页 / 共79页
第5页 / 共79页
第6页 / 共79页
第7页 / 共79页
第8页 / 共79页
王济-matlab在振动信号处理中的应用代码.docx
% 关闭所有隐藏的窗口
% 清除内存中所有变量和函数
% 清除工作窗口中所显示的内容
程序 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)';
%计算趋势项的多项式待定系数向量 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);
%在图幅上添加坐标网格
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 %f\n',t(k),y(k));
%循环体结束语句
end
%关闭数据文件
status=fclose(fid);
程序 4-2
%五点滑动平均法平滑处理
%%%%%%%%%%%%%%%%%%%%%%
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)';
%将 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 赋值给 y
y=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 %f\n',t(k),y(k));
end
status=fclose(fid);
程序 4-2
%五点滑动平均法平滑处理
%%%%%%%%%%%%%%%%%%%%%%
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)';
%将 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 赋值给 y
y=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 %f\n',t(k),y(k));
end
status=fclose(fid);
%采样频率
%滑动阶次
%平滑次数
程序 4-3
%滑动平均法消除趋势项
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close 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 减去与平滑趋势项 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('输入','输出','趋势');
%在图幅上添加坐标网格
grid on;
%以写的方式打开文件或建立一个新文件
fid=fopen(fno,'w');
%进行 n 次循环将结果写到输出数据文件中
for k=1:n
%每行写两个实型数据，t 为时间，y 为消除趋势项后的结果
fprintf(fid,'%f %f\n',t(k),y(k));
end
status=fclose(fid);
程序 4-4
%五点三次法平滑处理(时域和频域)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('五点三次平滑处理-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1);
it = fscanf(fid,'%d',1);
%平滑次数
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,:);
%采样频率
%数据类型(1 时域,2 频域)
%循环 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 %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
close 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=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);
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
