matlab 信号处理

1、一号楼楼板在环境激励下的振动信号处理振动是一种普遍存在的自然现象。为了了解结构对振动的反应，我们对一号楼楼板进行了测验，本文将对所测得的数据进行处理。具体包括振动信号的预处理、时域处理频域处理、数字信号的生成、试验模态参数的频域识别和时域识别以及整体识别。本文将对所测得的加速度数据进行处理。一、 振动信号预处理振动信号预处理是将振动测试中采集到的数据尽可能真实地还原成实际振动状况的最基本的数据加工方法。下面将对其进行消除多项式趋势项处理，采样数据的平滑处理。1、 消除多项式趋势项处理在振动测试中采集到的振动信号数据，由于放大器随温度变化产生的零点漂移、传感器频率范围外低频性能的不稳定以及传感器

2、周围的环境干扰，往往会偏离基线，甚至偏离基线的大小还会随时间变化。偏离基线随时间变化的整个过程被称为信号的趋势项。趋势项直接影响信号的正确性，应将其去除。本文将借助MATLAB进行最小二乘法消除多项式趋势项。a、最小二乘法消除多项式趋势项的程序 %clear % 清除内存中所有变量和函数clc % 清除工作窗口中所显示的内容close all hidden % 关闭所有隐藏的窗口%提示用键盘输入输入数据文件名fni=input('消除多项式趋势项-输入数据文件名:','s'); %以只读方式打开数据文件fid=fopen(fni,'r'); sf

3、 = 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)' %计算趋势项的多项式待定系数向量aa=polyfit(t,x,m); %用x减去多项式系数a

4、生成的趋势项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为消除趋势项

5、后的结果 fprintf(fid,'%f %fn',t(k),y(k); %循环体结束语句 end %关闭数据文件status=fclose(fid); b、处理结果2、 采样数据的平滑处理通过数据采集器采样的道德振动信号数据往往叠加有噪声信号。由于随机干扰信号的平带较宽，有时高频成分所占比例很大，使得采集到离散数据绘成的振动曲线上呈现许多毛刺，很不光滑。为了消除干扰信号的影响，本文进行了五点滑动平均法平滑处理。a、五点滑动平均法平滑处理的程序%clearclcclose all hidden%fni=input('五点滑动平均法平滑处理-输入数据文件名:',&

6、#39;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次进行

7、平滑处理计算 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赋值给yy=a;%将分成2行1列的图形窗口的第1列设为当前绘图区域subplot(2,1,1); %绘制平滑前的时程曲线图形

8、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&

9、#39;); for k=1:n %每行写两个实型数据，t为时间，y为平滑后的数据 fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); b、处理结果二、 振动信号时域处理振动信号时域处理是对振动波形的分析。从记录的时程信号中提取各种有用的信息或将记录的时程信号转换成所需要的形式。通过不同时域处理方法，可以确定实测波形的最大幅值和时间历程，求出相位滞后和波形的时间滞后，有选择地滤除或保留实测波形的某些频率成分，消除波形的畸变状况，再现真实波形面貌。数字滤波是通过数学手段从所采集的离散信号中选取人们所感兴趣的一部分信号的处

10、理方法。它的主要作用有滤除测试信号中的噪声或虚假成分、提高信噪比、平滑分析数据、抑制干扰信号、分离频率分量等。数字进入滤波器后，部分频率可以通过，部分则受阻挡。1、数字滤波（1）、频域低通和带通滤波a、频域低通和带通滤波的程序%clearclcclose all hidden%fni=input('频域带通滤波-输入数据文件名:','s'); fid = fopen(fni,'r'); sf = fscanf(fid,'%f',1); %采样频率fmin = fscanf(fid,'%f',1); %最小截止频率f

11、max = 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)' %取大于并

12、最接近n的2的幂次方为FFT长度nfft=2nextpow2(n); %四舍五入取整求最小截止频率对应数组元素的下标ni=round(fmin*nfft/sf+1); %四舍五入取整求最大截止频率对应数组元素的下标na=round(fmax*nfft/sf+1); %进行FFT变换,结果存于yy=fft(x,nfft); %建立一个长度为nfft元素全为0的向量 a=zeros(1,nfft); %将y的正频率带通内的元素赋值给a a(ni:na)=y(ni:na); %将y的负频率带通内的元素赋值给aa(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);

13、 %进行FFT逆变换,结果存于yy=ifft(a,nfft); %取逆变换的实部n个元素为滤波结果列向量y=(real(y(1:n)' %绘制滤波前的时程曲线图形subplot(2,1,1); plot(t,x); %添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on; %绘制滤波后的时程曲线图形subplot(2,1,2); plot(t,y); %添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on;%打开文件输出滤波后的数据 fid=fopen(fno,'w'

14、);for k=1:n fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); b、处理结果（2）、频域高通和带阻滤波a、频域高通和带阻滤波都的程序%频域高通和带阻滤波%clearclcclose all hidden%fni=input('频域带阻滤波-输入数据文件名:','s'); fid=fopen(fni,'r'); sf = fscanf(fid,'%f',1); %采样频率fmin = fscanf(fid,'%f',1); %最小

15、截止频率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)'

16、%取大于并最接近n的2的幂次方为FFT长度nfft=2nextpow2(n); %四舍五入取整求最小截止频率对应数组元素的下标ni=round(fmin*nfft/sf+1); %四舍五入取整求最大截止频率对应数组元素的下标na=round(fmax*nfft/sf+1); %进行FFT变换,结果存于yy=fft(x,nfft); %建立一个长度为nfft元素全为0的向量 a=zeros(1,nfft); %将y的正频率带阻内的元素赋值为0y(ni:na)=a(ni:na); %将y的负频率带阻内的元素赋值为0y(nfft-na+1:nfft-ni+1)=a(nfft-na+1:nfft-ni

17、+1); %进行IFFT变换,结果存于aa=ifft(y,nfft); %取逆变换的实部n个元素为滤波结果列向量y=real(a(1:n);%绘制滤波前的时程曲线图形subplot(2,1,1);plot(t,x); %添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on;%绘制滤波后的时程曲线图形subplot(2,1,2); plot(t,y);%添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注ylabel(sy); grid on; %打开文件输出滤波后的数据 fid=fopen(fno,'w'); f

18、or k=1:n fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); b、处理结果2、振动信号的积分和微分变换在振动信号测试过程中，由于仪器设备或测试环境的限制，有的物理量往往需要通过对采集到的其他物理量进行转换处理才能得到。频域积分a、频域积分的程序%频域积分%clearclcclose all hidden%fni=input('频域积分-输入数据文件名:','s'); fid=fopen(fni,'r'); sf = fscanf(fid,'%f',

19、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); %纵向坐标

20、轴输出单位的标注 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=2nextpow2(n); %FFT变换y=fft(x,nfft); %计算频率间隔（Hz/s）df=sf/nfft;%计算指定频带对应频率数组的下标ni=round(fmin/df+1); na=round(fmax/df+1);%计算

21、圆频率间隔（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); if it=2 %进行二次积分的相位变换 y=-a; else %进行一次积分的相位变换 real(y)=imag(a); imag(y)=-real(a); en

22、da=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); grid on; %绘制积分后的时程曲线图形 subplot(

23、2,1,2); plot(t,y);%添加横向坐标轴的标注xlabel(sx); %添加纵向坐标轴的标注 ylabel(sy2); grid on; %打开文件输出积分后的数据 fid=fopen(fno,'w'); for k=1:n fprintf(fid,'%f %fn',t(k),y(k); endstatus=fclose(fid);（b）、处理结果三、 振动信号频域处理我们在此对所测数据进行谱分析处理。频域处理也称谱分析，是建立在傅利叶变换基础上的时频变换处理，所得的结果是以频率为变量的函数，称谱函数。下面进行信号谱分析（a）、信号谱分析的程序%随机

24、信号谱分析%clearclcclose all hiddenformat long%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 =

25、 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; %加窗处理switch win case 1 %矩形窗 w=boxcar(nfft); case 2 %汉宁窗 w=hanning(nfft); cas

26、e 3 %海明窗 w=hamming(nfft); case 4 %布莱克曼窗 w=blackman(nfft); case 5 %三角窗 w=triang(nfft); otherwise w=boxcar(nfft);endswitch fun case 1 %计算自谱函数 z=psd(y,nfft,sf,w,nfft/2); case 2 %计算互谱函数 z=csd(x,y,nfft,sf,w,nfft/2); case 3 %计算频响函数 z=tfe(x,y,nfft,sf,w,nfft/2); case 4 %计算相干函数 z=cohere(x,y,nfft,sf,w,nfft/2)

27、; otherwise ;end%绘制幅频曲线图 nn=1:nfft/4;subplot(2,1,1);plot(f(nn),abs(z(nn);xlabel('频率 (Hz)'); ylabel('幅值'); grid on; if fun>1&fun<4 %绘制相频曲线图 subplot(2,1,2); plot(f(nn),angle(z(nn); xlabel('频率 (Hz)'); ylabel('相位'); grid on; end%以写的方式打开文件或建立一个新文件fid=fopen(fno,&#

28、39;w'); %输出自谱和相干函数数据if fun>1&fun<4 for k=1:nfft/2%每行输出2个实型数据，f为频率，z为幅值 fprintf(fid,'%f %fn',f(k),abs(z(k); end%输出互谱和频响函数数据else for k=1:nfft/2%每行输出3个实型数据，f为频率，z的实、虚部 fprintf(fid,'%f %f %fn',f(k),real(z(k),imag(z(k); endendstatus=fclose(fid);（b）、处理结果无3、 生成白噪声信号（a）、程序%生成白噪

29、声信号%clearclcclose all hidden%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);

30、%最大幅值fno = fscanf(fid,'%s',1); %输出数据文件名status=fclose(fid); %计算白噪声信号数据长度nn=round(sf*tl)+1; %建立信号离散时间向量t t=0:1/sf:(n-1)/sf; %大于并最接近n的2的幂次方为FFT长度 nfft=2nextpow2(n); %四舍五入取整求最小截止频率对应数组元素的下标ni=round(fi*nfft/sf+1); %四舍五入取整求最大截止频率对应数组元素的下标na=round(fa*nfft/sf+1); %建立一个长度为nfft/2元素全为0的向量 r=zeros(1,nff

31、t/2); %将r的正频率带通内的元素赋值1，生成幅值谱 r(ni:na)=1; %生成02pi的随机数为随机相位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

32、,boxcar(nft),nft/2); %定义自功率谱的离散频率向量 f=0:sf/nft:(nft/2-1)*sf/nft; %绘制白噪声信号随时间变化的曲线图subplot(2,1,1); plot(t,y);xlabel('时间(s)'); ylabel('幅值'); grid on;%绘制自功率谱曲线图subplot(2,1,2); plot(f(1:nft/4),y1(1:nft/4); xlabel('频率(Hz)'); ylabel('幅值'); grid on;%打开文件输出生成的白噪声信号数据 fid=fope

33、n(fno,'w'); for k=1:n %每一行输出两个实型数据，t为时间，y为白噪声信号值 fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); （b）、运行结果四、 试验模态参数识别的频域识别试验模态参数识别是振动信号处理的一个重要的组成部分，它的主要任务是从测试所得的数据中，确定振动系统的模态参数的估计，其中包括模态固有频率、模态阻尼比、模态振形以及模态质量和模态刚度等。下面进行有理分式多项式法模态参数识别（a） 程序%有理分式多项式法模态参数识别%clearclcclose all hidde

34、nformat long%fni=input('有理分式多项式法模态参数识别-输入数据文件名:','s'); fid=fopen(fni,'r'); mn = fscanf(fid,'%d',1); %模态阶数df = fscanf(fid,'%f',1); %频率间隔fno = fscanf(fid,'%s',1); %输出数据文件名b = fscanf(fid,'%f',2,inf); %实测频响函数实、虚部数据status=fclose(fid);%定义幂多项式的阶次nm =

35、mn*2;%取频响函数的长度n=length(b(1,:);%建立离散频率向量f=0:df:(n-1)*df;%建立离散圆频率向量w=2*pi*f; %建立归一化离散频率向量wi = w/max(w);%建立实测频响函数复数向量H=b(1,:)+i*b(2,:); %计算拟合频响函数的分子和分母系数向量A,B = invfreqs(H,wi,nm,nm,100);%幂多项式方程求根(零点)P = roots(B);%计算模态频率向量 F1 = abs(P)*max(w)/(2*pi);%计算阻尼比向量 D1 = -real(P)./(abs(P);%计算振型系数向量for k=1:nm if

36、k=1 p(1:nm-1)=P(2:nm); else p(1:k-1)=P(1:k-1); p(k:nm-1)=P(k+1:nm); end y=poly(p); S1(k)=polyval(A,P(k)/polyval(y,P(k); end%计算拟合的频响函数H1 = freqs(A,B,wi);%绘制频响函数实部拟合曲线图figure(2)nn=1:n;subplot(2,1,1);plot(f,real(H(nn),':',f,real(H1(nn),'r');xlabel('频率 (Hz)'); ylabel('实部'

37、); legend('实测','拟合');grid on;%绘制频响函数虚部拟合曲线图subplot(2,1,2);plot(f,b(2,:),':',f,imag(H1(nn),'r');xlabel('频率 (Hz)'); ylabel('虚部'); legend('实测','拟合'); grid on;%将模态频率从小到大排序F2,I=sort(F1); %剔除方程解中的非模态项(非共轭根)和共轭项(重复项)m=0;for k=1:nm-1 if F2(k)=F2(k+1) continue; end m=m+1; l=I(k); F(m)=F1(l); %模态频率 D(m)=D1(l); %阻尼比 S(m)=S1(l); %振型系数 end%打开文件输出模态参数数据fid=fopen