matlab信号处理.docx
《matlab信号处理.docx》由会员分享，可在线阅读，更多相关《matlab信号处理.docx（29页珍藏版）》请在冰点文库上搜索。
matlab信号处理
一号楼楼板在环境激励下的振动信号处理
振动是一种普遍存在的自然现象。
为了了解结构对振动的反应，我们对一号楼楼板进行了测验，本文将对所测得的数据进行处理。
具体包括振动信号的预处理、时域处理频域处理、数字信号的生成、试验模态参数的频域识别和时域识别以及整体识别。
本文将对所测得的加速度数据进行处理。
一、振动信号预处理
振动信号预处理是将振动测试中采集到的数据尽可能真实地还原成实际振动状况的最基本的数据加工方法。
下面将对其进行消除多项式趋势项处理，采样数据的平滑处理。
1、消除多项式趋势项处理
在振动测试中采集到的振动信号数据，由于放大器随温度变化产生的零点漂移、传感器频率范围外低频性能的不稳定以及传感器周围的环境干扰，往往会偏离基线，甚至偏离基线的大小还会随时间变化。
偏离基线随时间变化的整个过程被称为信号的趋势项。
趋势项直接影响信号的正确性，应将其去除。
本文将借助MATLAB进行最小二乘法消除多项式趋势项。
a、最小二乘法消除多项式趋势项的程序
%%%%%%%%%%%%%%%%%%%%%%%%
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）;
b、处理结果
2、采样数据的平滑处理
通过数据采集器采样的道德振动信号数据往往叠加有噪声信号。
由于随机干扰信号的平带较宽，有时高频成分所占比例很大，使得采集到离散数据绘成的振动曲线上呈现许多毛刺，很不光滑。
为了消除干扰信号的影响，本文进行了五点滑动平均法平滑处理。
a、五点滑动平均法平滑处理的程序
%%%%%%%%%%%%%%%%%%%%%%
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）;
b、处理结果
二、振动信号时域处理
振动信号时域处理是对振动波形的分析。
从记录的时程信号中提取各种有用的信息或将记录的时程信号转换成所需要的形式。
通过不同时域处理方法，可以确定实测波形的最大幅值和时间历程，求出相位滞后和波形的时间滞后，有选择地滤除或保留实测波形的某些频率成分，消除波形的畸变状况，再现真实波形面貌。
数字滤波是通过数学手段从所采集的离散信号中选取人们所感兴趣的一部分信号的处理方法。
它的主要作用有滤除测试信号中的噪声或虚假成分、提高信噪比、平滑分析数据、抑制干扰信号、分离频率分量等。
数字进入滤波器后，部分频率可以通过，部分则受阻挡。
1、数字滤波
（1）、频域低通和带通滤波
a、频域低通和带通滤波的程序
%%%%%%%%%%%%%%%%%%%%%%
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）;
b、处理结果
（2）、频域高通和带阻滤波
a、频域高通和带阻滤波都的程序
%频域高通和带阻滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%
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）;
b、处理结果
2、振动信号的积分和微分变换
在振动信号测试过程中，由于仪器设备或测试环境的限制，有的物理量往往需要通过对采集到的其他物理量进行转换处理才能得到。
频域积分
a、频域积分的程序
%频域积分
%%%%%%%%%%%%%%%%%%%%%%%%%%%
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）;%最大截止频率
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=2^nextpow2（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）;
ifit==2
%进行二次积分的相位变换
y=-a;
else
%进行一次积分的相位变换
real（y）=imag（a）;
imag（y）=-real（a）;
end
a=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）;
gridon;
%绘制积分后的时程曲线图形
subplot（2,1,2）;
plot（t,y）;
%添加横向坐标轴的标注
xlabel（sx）;
%添加纵向坐标轴的标注
ylabel（sy2）;
gridon;
%打开文件输出积分后的数据
fid=fopen（fno,'w'）;
fork=1:
n
fprintf（fid,'%f%f\n',t（k）,y（k））;
end
status=fclose（fid）;
（b）、处理结果
三、振动信号频域处理
我们在此对所测数据进行谱分析处理。
频域处理也称谱分析，是建立在傅利叶变换基础上的时频变换处理，所得的结果是以频率为变量的函数，称谱函数。
下面进行信号谱分析
（a）、信号谱分析的程序
%随机信号谱分析
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
closeallhidden
formatlong
%%%%%%%%%%%%%%%%%%%%%%
fni=input（'随机信号谱分析-输入数据文件名:
','s'）;
fid=fopen（fni,'r'）;
%分析内容（1=自谱,2=互谱,3=频响函数,4=相干函数）
fun=fscanf（fid,'%d',1）;
sf=fscanf（fid,'%f',1）;%采样频率
nfft=fscanf（fid,'%d',1）;%FFT长度
%窗函数（1=矩形,2=汉宁,3=海明,4=布莱克曼,5=三角）
win=fscanf（fid,'%d',1）;
fno=fscanf（fid,'%s',1）;%输出数据文件名
%按行读入数据,第1行激励,第2行响应
a=fscanf（fid,'%f',[2,inf]）;
status=fclose（fid）;
%取激励信号向量
x=a（1,:
）;
%取相对的响应信号向量
y=a（2,:
）-a（1,:
）;
%建立频率向量
f=0:
sf/nfft:
sf/2-sf/nfft;
%加窗处理
switchwin
case1%矩形窗
w=boxcar（nfft）;
case2%汉宁窗
w=hanning（nfft）;
case3%海明窗
w=hamming（nfft）;
case4%布莱克曼窗
w=blackman（nfft）;
case5%三角窗
w=triang（nfft）;
otherwise
w=boxcar（nfft）;
end
switchfun
case1%计算自谱函数
z=psd（y,nfft,sf,w,nfft/2）;
case2%计算互谱函数
z=csd（x,y,nfft,sf,w,nfft/2）;
case3%计算频响函数
z=tfe（x,y,nfft,sf,w,nfft/2）;
case4%计算相干函数
z=cohere（x,y,nfft,sf,w,nfft/2）;
otherwise
;
end
%绘制幅频曲线图
nn=1:
nfft/4;
subplot（2,1,1）;
plot（f（nn）,abs（z（nn）））;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
gridon;
iffun>1&fun<4
%绘制相频曲线图
subplot（2,1,2）;
plot（f（nn）,angle（z（nn）））;
xlabel（'频率（Hz）'）;
ylabel（'相位'）;
gridon;
end
%以写的方式打开文件或建立一个新文件
fid=fopen（fno,'w'）;
%输出自谱和相干函数数据
iffun>1&fun<4
fork=1:
nfft/2
%每行输出2个实型数据，f为频率，z为幅值
fprintf（fid,'%f%f\n',f（k）,abs（z（k）））;
end
%输出互谱和频响函数数据
else
fork=1:
nfft/2
%每行输出3个实型数据，f为频率，z的实、虚部
fprintf（fid,'%f%f%f\n',f（k）,real（z（k））,imag（z（k）））;
end
end
status=fclose（fid）;
（b）、处理结果
无
3、生成白噪声信号
（a）、程序
%生成白噪声信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
closeallhidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input（'生成白噪声信号-输入数据文件名:
','s'）;
fid=fopen（fni,'r'）;
sf=fscanf（fid,'%f',1）;%采样频率
fi=fscanf（fid,'%f',1）;%最小截止频率
fa=fscanf（fid,'%f',1）;%最大截止频率
tl=fscanf（fid,'%f',1）;%时间长度（秒）
am=fscanf（fid,'%f',1）;%最大幅值
fno=fscanf（fid,'%s',1）;%输出数据文件名
status=fclose（fid）;
%计算白噪声信号数据长度n
n=round（sf*tl）+1;
%建立信号离散时间向量t
t=0:
1/sf:
（n-1）/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2（n）;
%四舍五入取整求最小截止频率对应数组元素的下标
ni=round（fi*nfft/sf+1）;
%四舍五入取整求最大截止频率对应数组元素的下标
na=round（fa*nfft/sf+1）;
%建立一个长度为nfft/2元素全为0的向量
r=zeros（1,nfft/2）;
%将r的正频率带通内的元素赋值1，生成幅值谱
r（ni:
na）=1;
%生成0~2pi的随机数为随机相位
a=2*pi*rand（1,nfft/2）;
%将幅值谱和相位谱转换成实部和虚部
b=r.*exp（i*a）*nfft/2;
%将正负圆频率向量组合成一个向量
a=[b,b（nfft/2:
-1:
1）];
%IFFT变换，并取实部为生成白噪声信号
x=real（ifft（a,nfft））;
%取指定长度的数据并使数据的最大值为指定的幅值
y=am*x（1:
n）/max（abs（x（1:
n）））;
%定义自功率谱的FFT长度
nft=1024;
%计算生成白噪声信号自功率谱
y1=psd（y,nft,sf,boxcar（nft）,nft/2）;
%定义自功率谱的离散频率向量
f=0:
sf/nft:
（nft/2-1）*sf/nft;
升级会员 