Matlab振动程序-代码作业
Matlab振动程序-代码作业
一、课题任务要求随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。
本文为机械测试信号的时域和频域分析，其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。
在进行上述分析之前先要对振动信号进行拟合。
机械振动分为确定性振动和随机振动，确定性振动又分为周期振动和非周期振动，周期振动又进一步分为简谐振动和复杂的周期振动。
所以可以根据上述的分类来拟合振动信号。
在设计信号的处理程序时，用MATLAB中的相关函数来对所拟合的振动信号进行时域分析和频域分析，并对绘出的频谱图进行说明。
二、技术路线对机械振动信号的时域和频域采集，根据所拟合的振动信号，选取所需要的时域性能指标和频域分析的性能指标对振动信号进行分析。
其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。
现构造一个振动信号（在该程序中以两个衰减振动分量和一个随机数rand之和来拟合振动信号)，再利用MATLAB中的函数mean ()、min()、max()、std()对离散序列中的平均值、最大值、最小值、标准差等时域性能进行分析，通过调用函数fft(y)； psd(y)； rcep (y)对该振动信号进行频域内的性能分析。
在设计过程中的理论知识有离散傅立叶变换(DFT)、功率谱的概念和意义以及倒谱的概念和意义。
①DFT的定义和意义：DFT的定义式为：X(k)=DFT[x(n)]=t-1 x( n)卬如0WkWN-1n=0 NX(n)= 1 七-1X (k )W - kn =x(n) 0WnWN-1N Nk=0DFT的意义：DFT的意义在于它表示信号中的各个频率的分量的振动幅值的大小，亦即该分量对于振动信号影响的大小。
振动计权加速度级 matlab代码
振动计权加速度级 matlab代码一、前言振动计权加速度级是用来描述振动信号功率谱密度的一种量化参数，是在振动信号分析中经常使用的重要参数。
在振动信号分析中，使用matlab编程可以方便快捷地计算得到振动计权加速度级，本文将介绍如何使用matlab编程实现振动计权加速度级的计算。
二、振动计权加速度级的定义振动计权加速度级是指在一定频率范围内，振动信号功率谱密度在不同频率下面积的平方根。
它是把功率谱密度转化成一个单一的数值，可以代表整个频率范围内振动的能量大小，常用单位是m/s^2。
振动计权加速度级可以通过下式计算得到：```matlabfunction r = rms(x)r = sqrt(mean(x.^2));end```其中，x为振动信号，rms为均方根函数，可以用来计算振动信号的均方根值。
三、matlab代码实现1. 导入振动信号数据```matlabdata = load('vibration_data.txt'); 导入振动信号数据，假设数据保存在vibration_data.txt文件中```2. 计算功率谱密度```matlabFs = 1000; 采样频率为1000HzN = length(data); 数据长度xdft = fft(data); 振动信号的傅立叶变换xdft = xdft(1:N/2+1); 取一半长度psdx = (1/(Fs*N)) * abs(xdft).^2; 计算信号的功率谱密度psdx(2:end-1) = 2*psdx(2:end-1); 修正PSDfreq = 0:Fs/N:Fs/2; 频率范围```3. 计算振动计权加速度级```matlaba_w = sqrt(trapz(freq, psdx)); 计算振动计权加速度级```四、实例下面通过一个具体的实例来演示如何使用上述matlab代码计算振动计权加速度级。
假设导入的振动信号数据如下：```matlabvibration_data = [0.5, 1.2, 3.0, 2.5, 1.8, 1.0, 0.7, 0.3, -0.1, -0.5, -1.0, -1.5, -2.0, -2.5, -3.0];```按照上述matlab代码进行计算，则得到的振动计权加速度级为：```matlaba_w = 2.1; 假设的值，具体数值需要根据实际情况进行计算```五、总结本文介绍了振动计权加速度级的定义以及使用matlab编程实现振动计权加速度级的计算方法，并通过一个具体的实例进行了演示。
Matlab振动程序 代码作业
一、课题任务要求随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。
本文为机械测试信号的时域和频域分析，其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。
在进行上述分析之前先要对振动信号进行拟合。
机械振动分为确定性振动和随机振动，确定性振动又分为周期振动和非周期振动，周期振动又进一步分为简谐振动和复杂的周期振动。
所以可以根据上述的分类来拟合振动信号。
在设计信号的处理程序时，用MATLAB中的相关函数来对所拟合的振动信号进行时域分析和频域分析，并对绘出的频谱图进行说明。
二、技术路线对机械振动信号的时域和频域采集，根据所拟合的振动信号，选取所需要的时域性能指标和频域分析的性能指标对振动信号进行分析。
其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。
现构造一个振动信号（在该程序中以两个衰减振动分量和一个随机数rand 之和来拟合振动信号），再利用MATLAB中的函数mean（）、min（）、max（）、std （）对离散序列中的平均值、最大值、最小值、标准差等时域性能进行分析，通过调用函数fft（y）；psd（y）；rcep（y）对该振动信号进行频域内的性能分析。
在设计过程中的理论知识有离散傅立叶变换（DFT）、功率谱的概念和意义以及倒谱的概念和意义。
①DFT的定义和意义：DFT的定义式为：DFT 的意义：DFT 的意义在于它表示信号中的各个频率的分量的振动幅值的大小，亦即该分量对于振动信号影响的大小。
通过DFT 的快速算法FFT 可以很方便的将振动信号的各个分量的幅值比重计算出来。
即可以把信号的主频分量提取出来。
②功率谱的概念和意义：功率谱的定义式为:若X(Ω)=DFT[x(m)]，x(n)为N 点序列。
王济-matlab在振动信号处理中的应用代码
程序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)';%计算趋势项的多项式待定系数向量aa=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%五点滑动平均法平滑处理%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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赋值给aa=x;%循环m次进行平滑处理计算for k=1:mb(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-2b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;endb(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);%绘制平滑前的时程曲线图形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));endstatus=fclose(fid);程序4-2%五点滑动平均法平滑处理%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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赋值给aa=x;%循环m次进行平滑处理计算for k=1:mb(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-2b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;endb(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);%绘制平滑前的时程曲线图形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));endstatus=fclose(fid);程序4-3%滑动平均法消除趋势项%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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:mfor j=l+1:n-lb(j)=mean(a(j-l:j+l));enda=b;end%用输入信号x减去与平滑趋势项ay=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));endstatus=fclose(fid);程序4-4%五点三次法平滑处理(时域和频域) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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,:));for l=1:ita=x(l,:);%循环m次进行平滑处理for k=1:mb(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-2b(j)=(-3*(a(j-2)+a(j+2))+12*(a(j-1)+a(j+1))+17*a(j))/35;endb(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;endy(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:nif 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));endendstatus=fclose(fid);程序5-1%频域低通和带通滤波%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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变换,结果存于yy=fft(x,nfft);%建立一个长度为nfft元素全为0的向量a=zeros(1,nfft);%将y的正频率带通的元素赋值给aa(ni:na)=y(ni:na);%将y的负频率带通的元素赋值给aa(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);%进行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');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-2%频域高通和带阻滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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变换,结果存于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+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');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-3%运用完全设计法的IIR滤波器滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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);if fun<=2wp = fscanf(fid,'%f',1);%通带截止频率(Hz)ws = fscanf(fid,'%f',1);%阻带截止频率(Hz)elsewp = fscanf(fid,'%f',2);%通带截止频率(Hz)ws = fscanf(fid,'%f',2);%阻带截止频率(Hz)endrp = 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);%确定滤波方式switch funcase 1 %低通ft='low';case 2 %高通ft='high';case 3 %带通ft='bandpass';case 4 %带阻ft='stop';otherwiseft='low';end%根据滤波器种类进行IIR滤波器设计switch mod%巴特沃斯滤波器case 1[n wn]=buttord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= butter(n,wn,ft);%切比雪夫Ⅰ型滤波器case 2[n wn]=cheb1ord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= cheby1(n,rp,wn,ft);%切比雪夫Ⅱ型滤波器case 3[n wn]=cheb2ord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= cheby2(n,rs,wn,ft);%椭圆滤波器case 4[n wn]=ellipbuttord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= ellip(n,rp,rs,wn,ft);end%计算滤波器的频率响应[h w]=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('幅值');%在图幅上添加坐标网格grid on;%绘制滤波前后的信号时程图subplot(2,1,2);plot(t,x,':',t,y);%添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('加速度(g)');%在图幅上添加图例legend('滤波前','滤波后');grid on;%打开文件输出滤波后的数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f \n',t(k),y(k));endstatus=fclose(fid);程序5-4%多带FIR滤波器滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all %%%%%%%%%%%%%%%%%%%%%%%%%%%fni=input('多带FIR滤波器滤波-输入数据文件名:','s');fid=fopen(fni,'r');fs = fscanf(fid,'%f',1); %采样频率%滤波器函数选择（1=fir2,2=firls,3=remez,4=cremez）mod = fscanf(fid,'%d',1);n = fscanf(fid,'%d',1); %滤波器阶数k = fscanf(fid,'%d',1); %频率点数f = fscanf(fid,'%f',k); %指定频率点(Hz)%通阻状态(0=阻or 1=通)，对于函数cremez为期望幅值a = fscanf(fid,'%f',k);w = fscanf(fid,'%f',k/2); %权向量fno = fscanf(fid,'%s',1); %输出数据文件名%输入数据存成行向量（1加噪信号，2原始信号）x = fscanf(fid,'%f',[2,inf]);status=fclose(fid);m=length(x(1,:));t=0:1/fs:(m-1)/fs;%频率归一化处理f=2*f/fs;%多带FIR滤波器设计switch modcase 1 %fir2滤波函数b= fir2(n,f,a);case 2 %firls滤波函数b= firls(n,f,a,w);case 3 %remez滤波函数b= remez(n,f,a,w);case 4 %cremez滤波函数b= cremez(n,f,a,w,'real');end%用多带FIR滤波器进行滤波v=filter(b,1,x(1,:));%消除延时处理y=zeros(1,m);y(1:m-n/2)=v(n/2+1:m);%计算滤波器的频率响应[h w]=freqz(b,1,1024,fs);%绘制滤波器的频率响应图subplot(2,1,1)plot(w,abs(h));%添加横向坐标轴的标注xlabel('频率(Hz)');%添加纵向坐标轴的标注ylabel('幅值');grid on;%绘制滤波前后的信号时程比较图subplot(2,1,2)nn=1:200;plot(t(nn),x(1,nn),'-.',t(nn),x(2,nn),':',t(nn),y(nn)); %添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('幅值');%在图幅上添加图例legend('加噪信号','原始信号','滤波信号');grid on;%打开文件输出滤波后的数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f \n',t(k),y(k));endstatus=fclose(fid);程序5-5%频域积分%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%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);if it==2%进行二次积分的相位变换y=-a;else%进行一次积分的相位变换real(y)=imag(a);imag(y)=-real(a);enda=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(2,1,2);plot(t,y);%添加横向坐标轴的标注xlabel(sx);%添加纵向坐标轴的标注ylabel(sy2);grid on;%打开文件输出积分后的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-6%频域微分%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hiddenformat long %%%%%%%%%%%%%%%%%%%%%%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=y.*w;if it==2%进行二次微分的相位变换y=-a;else%进行一次微分的相位变换real(y)=imag(a);imag(y)=-real(a);enda=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(2,1,2);plot(t,y);%添加横向坐标轴的标注xlabel(sx);%添加纵向坐标轴的标注ylabel(sy2);grid on;%打开文件输出微分后的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-7%相关函数%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%fni=input('相关函数-输入数据文件名:','s');fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率%分析容(1=自相关,2=互相关)fun = fscanf(fid,'%d',1);m = fscanf(fid,'%d',1); %相关数据长度（点数）fno = fscanf(fid,'%s',1); %输出数据文件名x = fscanf(fid,'%f',[fun,inf]);%输入数据存成行向量status=fclose(fid);n=length(x(1,:));%建立离散时间向量t=0:1/sf:(n-1)/sf;if fun==1%计算自相关函数向量a=xcorr(x(1,:),m);else%计算互相关函数向量a=xcorr(x(1,:),x(2,:),m);end%取正频率段的相关函数y=a(m+1:2*m+1);%绘制输入数据的时程曲线图形subplot(2,1,1);nn=1:4000;if fun==1plot(t(nn),x(1,nn));elseplot(t(nn),x(1,nn),':',t(nn),x(2,nn));legend('x','y');endgrid on;%绘制相关函数曲线图形subplot(2,1,2);nn=1:m;plot(t(nn),y(nn));grid on;%打开文件输出相关函数数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序6-1%随机信号谱分析%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 = 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 wincase 1 %矩形窗w=boxcar(nfft);case 2 %汉宁窗w=hanning(nfft);case 3 %海明窗w=hamming(nfft);case 4 %布莱克曼窗w=blackman(nfft);case 5 %三角窗w=triang(nfft);otherwisew=boxcar(nfft);endswitch funcase 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);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,'w');%输出自谱和相干函数数据if fun>1&fun<4for k=1:nfft/2%每行输出2个实型数据，f为频率，z为幅值fprintf(fid,'%f %f\n',f(k),abs(z(k)));end%输出互谱和频响函数数据elsefor k=1:nfft/2%每行输出3个实型数据，f为频率，z的实、虚部fprintf(fid,'%f %f %f\n',f(k),real(z(k)),imag(z(k)));endendstatus=fclose(fid);程序6-2%窗函数（MATLAB中没有这些窗函数）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fni=input('窗函数-输入数据文件名:','s');fid=fopen(fni,'r');n = fscanf(fid,'%d',1); %数据长度%窗函数(1=余弦坡度窗,2=帕曾窗,3=指数窗,4=高斯窗) win = fscanf(fid,'%d',1);%参数(除指数类窗需要输入参数外其余输入1)c = fscanf(fid,'%f',1);fno = fscanf(fid,'%s',1);%输出数据文件名status=fclose(fid);%建立数据长度向量l=0:n-1;%计算窗函数switch wincase 1 %余弦坡度窗m=round(n*4/5);w=[ones(1,m),(1+cos(5*pi*l(m:n-1)/n))/2];case 2 %帕曾窗w=[1-6*(l(1:n/2)/n).^2+6*(l(1:n/2)/n).^3,2*(1-l(n/2+1:n)/n).^3]; case 3 %指数窗w=exp(-c*l);case 4 %高斯窗w=exp(-c*l.^2);end%绘制窗函数曲线图形plot(l/n,w);grid on;%打开文件输出窗函数的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%d %f\n',k,w(k));endstatus=fclose(fid);程序6-3%细化FFT %%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hiddenformat long %%%%%%%%%%%%%%%%%%%%%%fni=input('细化FFT处理-输入数据文件名(20041010):','s');fid=fopen(fni,'r');fs = fscanf(fid,'%f',1); %采样频率fi = fscanf(fid,'%f',1); %最小细化截止频率np = fscanf(fid,'%d',1); %放大倍数nfft = fscanf(fid,'%d',1); %FFT长度fno = fscanf(fid,'%s',1); %输出数据文件名x = fscanf(fid,'%f',[1,inf]); %按行读入原始信号数据status=fclose(fid);%计算输入数据长度nt=length(x);%最大细化截止频率fa=fi+0.5*fs/np;%取大于nt且与nt最接近的2的整数次方为FFT长度nf=2^nextpow2(nt);%确定细化带宽的数据长度na=round(0.5*nf/np+1);%频移%建立一个按1递增的向量n=0:nt-1;%确定单位旋转因子向量b=n*pi*(fi+fa)/fs;%乘以单位旋转因子进行频移y=x.*exp(-i*b);%频移后的低通滤波（频域滤波）%FFT变换b=fft(y,nf);%正频率带通的元素赋值a(1:na)=b(1:na);%负频率带通的元素赋值a(nf-na+1:nf)=b(nf-na+1:nf);%IFFT变换b=ifft(a,nf);%重采样c=b(1:np:nt);%进行细化FFT变换a=fft(c,nfft)*2/nfft;%变换结果重新排序y2=zeros(1,nfft/2);%排列负频率段的数据y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);%排列正频率段的数据y2(nfft/4+1:nfft/2)=a(1:nfft/4);n=0:(nfft/2-1);%定义细化FFT频率向量f2=fi+n*2*(fa-fi)/nfft;%对输入数据做FFT用来变换比较%FFT变换y1=fft(x,nfft)*2/nfft;f1=n*fs/nfft;%定义与细化FFT频率向量相同的频率围ni=round(fi*nfft/fs+1);na=round(fa*nfft/fs+1);%绘制输入时程曲线图形subplot(2,1,1);t=0:1/fs:(nt-1)/fs;plot(t,x);xlabel('时间(s)');ylabel('幅值');grid on;%在相同的频率围下绘制幅频曲线图nn=ni:na;plot(f1(nn),abs(y1(nn)),':',f2,abs(y2));xlabel('频率(Hz)');ylabel('幅值');legend('普通','细化');grid on;%打开文件输出细化后的数据fid=fopen(fno,'w');for k=1:nfft/2%每行输出3个实型数据，f为频率，y2的实和虚部fprintf(fid,'%f %f %f\n',f2(k),real(y2(k)),imag(y2(k))); endstatus=fclose(fid);程序6-4%三分之一倍频程处理%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%fni=input('三分之一倍频程处理-输入数据文件名:','s'); fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率fno = fscanf(fid,'%s',1); %输出数据文件名x = fscanf(fid,'%f',[1,inf]); %按行输入数据status=fclose(fid);%定义三分之一倍频程的中心频率f = [1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.300 8.00]; fc =[f,10*f,100*f,1000*f,10000*f];%中心频率与下限频率的比值oc6=2^(1/6);%取中心频率总的长度nc=length(fc);%输入数据的长度n=length(x);%大于并最接近n的2的幂次方长度nfft=2^nextpow2(n);%FFT变换a=fft(x,nfft);for j=1:nc%下限频率fl=fc(j)/oc6;%上限频率fu=fc(j)*oc6;%下限频率对应的序号nl=round(fl*nfft/sf+1);%上限频率对应的序号nu=round(fu*nfft/sf+1);%如果上限频率大于折叠频率则循环中断if fu > sf/2m=j-1; breakend%以每个中心频率段为通带进行带通频域滤波b=zeros(1,nfft);b(nl:nu)=a(nl:nu);b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);c=ifft(b,nfft);%计算对应每个中心频率段的有效值yc(j)=sqrt(var(real(b(1:n))));end%绘制输入时程曲线图形subplot(2,1,1);t=0:1/sf:(n-1)/sf;plot(t,x);xlabel('时间(s)');ylabel('加速度(g)');grid on;%绘制三分之一倍频程有效值图形subplot(2,1,2);plot(fc(1:m),yc(1:m));xlabel('频率(Hz)');ylabel('有效值');grid on;%打开文件输出三分之一倍频程数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f\n',fc(k),yc(k));endstatus=fclose(fid);程序6-5%倒频谱函数%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%fni=input('倒频谱函数-输入数据文件名:','s'); fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率%分析容(1=实倒频谱,2=复倒频谱)fun = 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;%计算倒频谱函数if fun==1y=rceps(x); %计算实倒频谱,elsey=cceps(x); %计算复倒频谱end%绘制输入时程曲线图形subplot(2,1,1);plot(t,x);xlabel('时间(s)');ylabel('幅值');grid on;%绘制倒频谱时程曲线图形subplot(2,1,2);plot(t,y);xlabel('时间(s)');ylabel('幅值');grid on;%打开文件输出倒频谱的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%d %f\n',t(k),y(k));endstatus=fclose(fid);程序6-6%地震信号加速度反应谱%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hiddenformat long %%%%%%%%%%%%%%%%%%%%%%fni=input('地震信号加速度反应谱-输入数据文件名:','s'); fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率dp = fscanf(fid,'%f',1); %阻尼比hf = fscanf(fid,'%f',1); %计算最高截止频率df = fscanf(fid,'%f',1); %频率间隔fno = fscanf(fid,'%s',1); %输出数据文件名xg = fscanf(fid,'%f',[1,inf]); %输入地震波数据存成行向量status=fclose(fid);%计算输入数据长度n=length(xg);%计算反应谱长度m=round(hf/df)+1;%定义频率向量f=0:df:(m-1)*df;%定义时间向量t=0:1/sf:(n-1)/sf;%计算加速度反应谱for j=0:m-1%计算单位脉冲反应w=2*pi*df*j;wd=w*sqrt(1-dp*dp);e=exp(-t*w*dp);a=t.*wd;s=sin(a).*((1-2*dp*dp)/(1-dp*dp));c=cos(a).*(2*dp/sqrt(1-dp*dp));h=wd*e.*(s+c)/sf;%用卷积计算加速度反应谱值r(j+1)=max(abs(conv(xg,h)));end%绘制输入地震波加速度曲线subplot(2,1,1);plot(t,xg);xlabel('时间(s)');ylabel('加速度(g)');grid on;%绘制地震波加速度反应谱曲线subplot(2,1,2);plot(f,r);xlabel('频率(Hz)');ylabel('加速度(g)');grid on;%打开文件输出地震波加速度反应谱数据fid=fopen(fno,'w');for k=1:m%每行输出两个实型数据，f为频率，r为反应谱值fprintf(fid,'%f %f\n',f(k),r(k));endstatus=fclose(fid);程序7-1%生成扫频信号%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all %%%%%%%%%%%%%%%%%%%%%%fni=input('生成扫频信号-输入数据文件名:','s');fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率f0 = fscanf(fid,'%f',1); %最小截止频率f1 = fscanf(fid,'%f',1); %最大截止频率t0 = fscanf(fid,'%f',1); %开始时刻(s)t1 = fscanf(fid,'%f',1); %终止时刻(s)am = fscanf(fid,'%f',1); %扫频信号幅值md = fscanf(fid,'%d',1); %扫频方式(1=线性,2=对数,3=指数) fno = fscanf(fid,'%s',1); %输出数据文件名status=fclose(fid);%四舍五入取整求扫频信号开始时刻对应数组元素的下标n0=round(sf*t0)+1;%四舍五入取整求扫频信号终止时刻对应数组元素的下标n1=round(sf*t1)+1;%建立离散时间向量tt=0:1/sf:(n1-1)/sf;%按指定的扫频方式计算随时间变化频率向量switch md%计算扫频方式为线性时随时间变化频率向量case 1ft=(t-t0)./(t1-t0);%计算扫频方式为对数时随时间变化频率向量case 2ft=log(t-t0)./log(t1-t0+1);%计算扫频方式为指数时随时间变化频率向量case 3ft=exp((t-t0)./(t1-t0)-1);otherwiseft=(t-t0)./(t1-t0);end%计算随时间变化的园频率向量w=2*pi*(ft.*(f1-f0)+f0);%建立一个长度与时间长度相同元素全为0的向量y=zeros(1,n1);%计算扫频信号y(n0:n1)=am*sin(w(n0:n1).*t(n0:n1));%绘制算扫频信号随时间变化的曲线图plot(t,y);xlabel('时间(s)');ylabel('幅值');grid on;%打开文件输出扫频信号数据fid=fopen(fno,'w');for k=1:n1%每一行输出两个实型数据，t为时间，y为扫频信号值fprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序7-2%生成拍波信号%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%fni=input('生成拍波信号-输入数据文件名:','s');fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率fp = fscanf(fid,'%f',1); %拍波频率np = fscanf(fid,'%d',1); %拍波个数mp = fscanf(fid,'%d',1); %每个拍波的周波数ti = fscanf(fid,'%f',1); %拍波之间的时间间隔am = fscanf(fid,'%f',1); %拍波幅值fno = fscanf(fid,'%s',1) %输出数据文件名status=fclose(fid);%建立完整拍波的离散时间向量ts=0:1/sf:(round(sf*mp/fp)-1)/sf;%建立拍波之间的离散时间间隔向量t0=zeros(1,round(ti*sf)+1);%计算拍波圆频率向量w=2*pi*fp;%生成单个拍波的离散信号y1=am*sin(w*ts).*sin(0.5*w*ts/mp);%建立时间间隔和单个拍波信号的组合向量y=[t0,y1];%通过循环建立完整的拍波信号向量for k=1:np-1y=[y,t0,y1];end%建立与整个拍波信号长度相同的离散时间向量t=0:1/sf:(length(y)-1)/sf;%绘制拍波信号随时间变化的曲线图plot(t,y);xlabel('时间(s)');ylabel('幅值');grid on;%打开文件输出生成的拍波信号数据fid=fopen(fno,'w');for k=1:length(y)%每一行输出两个实型数据，t为时间，y为拍波信号值fprintf(fid,'%f %f\n',t(k),y(k));。
简谐振动合成-Matlab
二、振动的合成实际生活中，一个系统往往会同时参与两个或更多的振动。
例如悬挂在颠簸船舱中的钟摆，两列声波同时传入人耳等。
一般的振动合成显然是比较复杂，下面仅讨论几种间单情况的简谐振动合成。
一、同方向同频率简谐振动的合成若两个同方向的简谐振动，频率都是，它们的运动方程分别为因振动是同方向的，所以这两个谐振动在任意时刻的和位移应在同一直线上，且等于这两个振动位移的代数和，即合位移仍为简谐振动二、两个同方向不同频率简谐振动的合成拍如果两个简谐振动的振动方向相同而频率不同，那么合成后的振动仍与原振动方向相同但不再是简谐振动。
现设两简谐振动的振幅都为A，初相位为零，它们的振动方程分别为合成振动方程为若两个分振动的频率都较大且其差很小时，即，合振动可看作为振幅随时间缓慢变化的近似谐振动，振幅随时间变化且具有周期性，表现出振动或强或弱的现象，称拍，变化的频率称拍频，变化的振幅为变化的频率为三、相互垂直的简谐振动的合成李萨如图如果两个简谐振动分别在x轴和y轴上进行，他们的振动方程分别为合成后，可得质点的轨迹为椭圆方程若两分振动有不同的频率，且两频率之比为有理数时，则合成后的质点运动具有稳定、封闭的轨迹。
称其为李萨如图形。
程序编写我们已经在第一讲中体验了matlab的编程，可是你一定会生出这样的问号，辛辛苦苦在命令窗口写的一大堆代码怎么不保留？不用担心，matlab程序和其他编程工具一样，也有专门的文件格式，称m文件，文件名形式为“文件名.m”。
你可以用matlab自带的编辑器来输入你的程序代码，当然你也可以用其它编辑器或最经济的文本编辑器，不过别忘记添加文件名的后缀“.m”。
下面，请跟我一起用m文件编辑器来编写matlab程序。
例题：两个振动方向相同而频率不同的简谐振动方程分别为合成后的方程是请用matlab程序描述合成波和拍频现象。
编程：第一步：点击matlab图标，打开程序窗口。
第二步：选file—new—m-file,打开编辑器。
面向振动的基于matlab的数据处理编程实现
面向振动的基于matlab的数据处理编程实现振动是物体在力的作用下发生的周期性的来回运动。
在工程领域中，振动的数据处理是非常重要的。
利用振动数据可以分析物体的结构特性、故障诊断以及设计和优化振动控制系统等。
本文将以基于MATLAB 的数据处理编程实现为主题，分为以下步骤进行讨论。
Step 1: 导入振动数据首先，我们需要将振动数据导入到MATLAB 环境中。
可以使用`load` 函数加载预先保存的数据文件，或使用`importdata` 函数读取文本文件、Excel 文件或其他常见的数据格式。
通过在MATLAB 命令窗口中输入相关命令，可以将数据存储在一个变量中以供后续处理使用。
Step 2: 数据预处理在进行振动数据处理之前，通常需要对数据进行预处理。
这包括去除噪声、滤波、数据对齐和裁剪等步骤。
可以使用MATLAB 中丰富的信号处理工具箱来实现这些操作。
例如，使用`butter` 函数可以设计一个巴特沃斯滤波器以去除高频噪声，或使用`medfilt1` 函数进行中值滤波。
此外，还可以使用`resample` 函数对数据进行采样率调整，以适应后续分析的需要。
Step 3: 频域分析频域分析是振动数据处理的重要步骤之一，可以通过它来确定振动信号的主要频率成分。
使用MATLAB 的信号处理工具箱中的傅里叶变换函数（如`fft`）可以将时域振动信号转换为频域。
通过对频域信号进行幅度谱和相位谱分析，可以确定振动信号的频谱和特征频率。
这些特征频率包括共振频率、自然频率、阻尼比等，对于结构特性和故障诊断非常重要。
Step 4: 时域分析时域分析是振动数据处理的另一个重要步骤，主要用于研究振动信号的时变特性。
其中，包络分析是一种常见的时域分析方法。
可以使用MATLAB 的信号处理工具箱中的函数（如`hilbert` 或`envelope`）对振动信号进行包络提取。
包络分析可以揭示振动信号的幅值变化规律，从而实现故障诊断和机械状态监测。
matlab在振动信号处理中的应用代码
程序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)';%计算趋势项的多项式待定系数向量aa=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%五点滑动平均法平滑处理%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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赋值给aa=x;%循环m次进行平滑处理计算for k=1:mb(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-2b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;endb(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);%绘制平滑前的时程曲线图形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));endstatus=fclose(fid);程序4-2%五点滑动平均法平滑处理%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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赋值给aa=x;%循环m次进行平滑处理计算for k=1:mb(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-2b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;endb(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);%绘制平滑前的时程曲线图形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));endstatus=fclose(fid);程序4-3%滑动平均法消除趋势项%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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:mfor j=l+1:n-lb(j)=mean(a(j-l:j+l));enda=b;end%用输入信号x减去与平滑趋势项ay=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));endstatus=fclose(fid);程序4-4%五点三次法平滑处理(时域和频域) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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,:));for l=1:ita=x(l,:);%循环m次进行平滑处理for k=1:mb(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-2b(j)=(-3*(a(j-2)+a(j+2))+12*(a(j-1)+a(j+1))+17*a(j))/35;endb(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;endy(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:nif 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));endendstatus=fclose(fid);程序5-1%频域低通和带通滤波%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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变换,结果存于yy=fft(x,nfft);%建立一个长度为nfft元素全为0的向量a=zeros(1,nfft);%将y的正频率带通内的元素赋值给aa(ni:na)=y(ni:na);%将y的负频率带通内的元素赋值给aa(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);%进行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');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-2%频域高通和带阻滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose 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变换,结果存于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+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');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-3%运用完全设计法的IIR滤波器滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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);if fun<=2wp = fscanf(fid,'%f',1);%通带截止频率(Hz)ws = fscanf(fid,'%f',1);%阻带截止频率(Hz)elsewp = fscanf(fid,'%f',2);%通带截止频率(Hz)ws = fscanf(fid,'%f',2);%阻带截止频率(Hz)endrp = 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);%确定滤波方式switch funcase 1 %低通ft='low';case 2 %高通ft='high';case 3 %带通ft='bandpass';case 4 %带阻ft='stop';otherwiseft='low';end%根据滤波器种类进行IIR滤波器设计switch mod%巴特沃斯滤波器case 1[n wn]=buttord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= butter(n,wn,ft);%切比雪夫Ⅰ型滤波器case 2[n wn]=cheb1ord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= cheby1(n,rp,wn,ft);%切比雪夫Ⅱ型滤波器case 3[n wn]=cheb2ord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= cheby2(n,rs,wn,ft);%椭圆滤波器case 4[n wn]=ellipbuttord(wp/(fs/2),ws/(fs/2),rp,rs);[b a]= ellip(n,rp,rs,wn,ft);end%计算滤波器的频率响应[h w]=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('幅值');%在图幅上添加坐标网格grid on;%绘制滤波前后的信号时程图subplot(2,1,2);plot(t,x,':',t,y);%添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('加速度(g)');%在图幅上添加图例legend('滤波前','滤波后');grid on;%打开文件输出滤波后的数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f \n',t(k),y(k));endstatus=fclose(fid);程序5-4%多带FIR滤波器滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all %%%%%%%%%%%%%%%%%%%%%%%%%%%fni=input('多带FIR滤波器滤波-输入数据文件名:','s');fid=fopen(fni,'r');fs = fscanf(fid,'%f',1); %采样频率%滤波器函数选择（1=fir2,2=firls,3=remez,4=cremez）mod = fscanf(fid,'%d',1);n = fscanf(fid,'%d',1); %滤波器阶数k = fscanf(fid,'%d',1); %频率点数f = fscanf(fid,'%f',k); %指定频率点(Hz)%通阻状态(0=阻or 1=通)，对于函数cremez为期望幅值a = fscanf(fid,'%f',k);w = fscanf(fid,'%f',k/2); %权向量fno = fscanf(fid,'%s',1); %输出数据文件名%输入数据存成行向量（1加噪信号，2原始信号）x = fscanf(fid,'%f',[2,inf]);status=fclose(fid);m=length(x(1,:));t=0:1/fs:(m-1)/fs;%频率归一化处理f=2*f/fs;%多带FIR滤波器设计switch modcase 1 %fir2滤波函数b= fir2(n,f,a);case 2 %firls滤波函数b= firls(n,f,a,w);case 3 %remez滤波函数b= remez(n,f,a,w);case 4 %cremez滤波函数b= cremez(n,f,a,w,'real');end%用多带FIR滤波器进行滤波v=filter(b,1,x(1,:));%消除延时处理y=zeros(1,m);y(1:m-n/2)=v(n/2+1:m);%计算滤波器的频率响应[h w]=freqz(b,1,1024,fs);%绘制滤波器的频率响应图subplot(2,1,1)plot(w,abs(h));%添加横向坐标轴的标注xlabel('频率(Hz)');%添加纵向坐标轴的标注ylabel('幅值');grid on;%绘制滤波前后的信号时程比较图subplot(2,1,2)nn=1:200;plot(t(nn),x(1,nn),'-.',t(nn),x(2,nn),':',t(nn),y(nn)); %添加横向坐标轴的标注xlabel('时间(s)');%添加纵向坐标轴的标注ylabel('幅值');%在图幅上添加图例legend('加噪信号','原始信号','滤波信号');grid on;%打开文件输出滤波后的数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f \n',t(k),y(k));endstatus=fclose(fid);程序5-5%频域积分%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%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);if it==2%进行二次积分的相位变换y=-a;else%进行一次积分的相位变换real(y)=imag(a);imag(y)=-real(a);enda=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(2,1,2);plot(t,y);%添加横向坐标轴的标注xlabel(sx);%添加纵向坐标轴的标注ylabel(sy2);grid on;%打开文件输出积分后的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-6%频域微分%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hiddenformat long %%%%%%%%%%%%%%%%%%%%%%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=y.*w;if it==2%进行二次微分的相位变换y=-a;else%进行一次微分的相位变换real(y)=imag(a);imag(y)=-real(a);enda=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(2,1,2);plot(t,y);%添加横向坐标轴的标注xlabel(sx);%添加纵向坐标轴的标注ylabel(sy2);grid on;%打开文件输出微分后的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序5-7%相关函数%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%fni=input('相关函数-输入数据文件名:','s');fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率%分析内容(1=自相关,2=互相关)fun = fscanf(fid,'%d',1);m = fscanf(fid,'%d',1); %相关数据长度（点数）fno = fscanf(fid,'%s',1); %输出数据文件名x = fscanf(fid,'%f',[fun,inf]);%输入数据存成行向量status=fclose(fid);n=length(x(1,:));%建立离散时间向量t=0:1/sf:(n-1)/sf;if fun==1%计算自相关函数向量a=xcorr(x(1,:),m);else%计算互相关函数向量a=xcorr(x(1,:),x(2,:),m);end%取正频率段的相关函数y=a(m+1:2*m+1);%绘制输入数据的时程曲线图形subplot(2,1,1);nn=1:4000;if fun==1plot(t(nn),x(1,nn));elseplot(t(nn),x(1,nn),':',t(nn),x(2,nn));legend('x','y');endgrid on;%绘制相关函数曲线图形subplot(2,1,2);nn=1:m;plot(t(nn),y(nn));grid on;%打开文件输出相关函数数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序6-1%随机信号谱分析%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 = 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 wincase 1 %矩形窗w=boxcar(nfft);case 2 %汉宁窗w=hanning(nfft);case 3 %海明窗w=hamming(nfft);case 4 %布莱克曼窗w=blackman(nfft);case 5 %三角窗w=triang(nfft);otherwisew=boxcar(nfft);endswitch funcase 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);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,'w');%输出自谱和相干函数数据if fun>1&fun<4for k=1:nfft/2%每行输出2个实型数据，f为频率，z为幅值fprintf(fid,'%f %f\n',f(k),abs(z(k)));end%输出互谱和频响函数数据elsefor k=1:nfft/2%每行输出3个实型数据，f为频率，z的实、虚部fprintf(fid,'%f %f %f\n',f(k),real(z(k)),imag(z(k)));endendstatus=fclose(fid);程序6-2%窗函数（MATLAB中没有这些窗函数）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fni=input('窗函数-输入数据文件名:','s');fid=fopen(fni,'r');n = fscanf(fid,'%d',1); %数据长度%窗函数(1=余弦坡度窗,2=帕曾窗,3=指数窗,4=高斯窗) win = fscanf(fid,'%d',1);%参数(除指数类窗需要输入参数外其余输入1)c = fscanf(fid,'%f',1);fno = fscanf(fid,'%s',1);%输出数据文件名status=fclose(fid);%建立数据长度向量l=0:n-1;%计算窗函数switch wincase 1 %余弦坡度窗m=round(n*4/5);w=[ones(1,m),(1+cos(5*pi*l(m:n-1)/n))/2];case 2 %帕曾窗w=[1-6*(l(1:n/2)/n).^2+6*(l(1:n/2)/n).^3,2*(1-l(n/2+1:n)/n).^3]; case 3 %指数窗w=exp(-c*l);case 4 %高斯窗w=exp(-c*l.^2);end%绘制窗函数曲线图形plot(l/n,w);grid on;%打开文件输出窗函数的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%d %f\n',k,w(k));endstatus=fclose(fid);程序6-3%细化FFT %%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hiddenformat long %%%%%%%%%%%%%%%%%%%%%%fni=input('细化FFT处理-输入数据文件名(20041010):','s');fid=fopen(fni,'r');fs = fscanf(fid,'%f',1); %采样频率fi = fscanf(fid,'%f',1); %最小细化截止频率np = fscanf(fid,'%d',1); %放大倍数nfft = fscanf(fid,'%d',1); %FFT长度fno = fscanf(fid,'%s',1); %输出数据文件名x = fscanf(fid,'%f',[1,inf]); %按行读入原始信号数据status=fclose(fid);%计算输入数据长度nt=length(x);%最大细化截止频率fa=fi+0.5*fs/np;%取大于nt且与nt最接近的2的整数次方为FFT长度nf=2^nextpow2(nt);%确定细化带宽的数据长度na=round(0.5*nf/np+1);%频移%建立一个按1递增的向量n=0:nt-1;%确定单位旋转因子向量b=n*pi*(fi+fa)/fs;%乘以单位旋转因子进行频移y=x.*exp(-i*b);%频移后的低通滤波（频域滤波）%FFT变换b=fft(y,nf);%正频率带通内的元素赋值a(1:na)=b(1:na);%负频率带通内的元素赋值a(nf-na+1:nf)=b(nf-na+1:nf);%IFFT变换b=ifft(a,nf);%重采样c=b(1:np:nt);%进行细化FFT变换a=fft(c,nfft)*2/nfft;%变换结果重新排序y2=zeros(1,nfft/2);%排列负频率段的数据y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);%排列正频率段的数据y2(nfft/4+1:nfft/2)=a(1:nfft/4);n=0:(nfft/2-1);%定义细化FFT频率向量f2=fi+n*2*(fa-fi)/nfft;%对输入数据做FFT用来变换比较%FFT变换y1=fft(x,nfft)*2/nfft;f1=n*fs/nfft;%定义与细化FFT频率向量相同的频率范围ni=round(fi*nfft/fs+1);na=round(fa*nfft/fs+1);%绘制输入时程曲线图形subplot(2,1,1);t=0:1/fs:(nt-1)/fs;plot(t,x);xlabel('时间(s)');ylabel('幅值');grid on;%在相同的频率范围下绘制幅频曲线图nn=ni:na;plot(f1(nn),abs(y1(nn)),':',f2,abs(y2));xlabel('频率(Hz)');ylabel('幅值');legend('普通','细化');grid on;%打开文件输出细化后的数据fid=fopen(fno,'w');for k=1:nfft/2%每行输出3个实型数据，f为频率，y2的实和虚部fprintf(fid,'%f %f %f\n',f2(k),real(y2(k)),imag(y2(k))); endstatus=fclose(fid);程序6-4%三分之一倍频程处理%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%fni=input('三分之一倍频程处理-输入数据文件名:','s'); fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率fno = fscanf(fid,'%s',1); %输出数据文件名x = fscanf(fid,'%f',[1,inf]); %按行输入数据status=fclose(fid);%定义三分之一倍频程的中心频率f = [1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.300 8.00]; fc =[f,10*f,100*f,1000*f,10000*f];%中心频率与下限频率的比值oc6=2^(1/6);%取中心频率总的长度nc=length(fc);%输入数据的长度n=length(x);%大于并最接近n的2的幂次方长度nfft=2^nextpow2(n);%FFT变换a=fft(x,nfft);for j=1:nc%下限频率fl=fc(j)/oc6;%上限频率fu=fc(j)*oc6;%下限频率对应的序号nl=round(fl*nfft/sf+1);%上限频率对应的序号nu=round(fu*nfft/sf+1);%如果上限频率大于折叠频率则循环中断if fu > sf/2m=j-1; breakend%以每个中心频率段为通带进行带通频域滤波b=zeros(1,nfft);b(nl:nu)=a(nl:nu);b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);c=ifft(b,nfft);%计算对应每个中心频率段的有效值yc(j)=sqrt(var(real(b(1:n))));end%绘制输入时程曲线图形subplot(2,1,1);t=0:1/sf:(n-1)/sf;plot(t,x);xlabel('时间(s)');ylabel('加速度(g)');grid on;%绘制三分之一倍频程有效值图形subplot(2,1,2);plot(fc(1:m),yc(1:m));xlabel('频率(Hz)');ylabel('有效值');grid on;%打开文件输出三分之一倍频程数据fid=fopen(fno,'w');for k=1:mfprintf(fid,'%f %f\n',fc(k),yc(k));endstatus=fclose(fid);程序6-5%倒频谱函数%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%fni=input('倒频谱函数-输入数据文件名:','s'); fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率%分析内容(1=实倒频谱,2=复倒频谱)fun = 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;%计算倒频谱函数if fun==1y=rceps(x); %计算实倒频谱,elsey=cceps(x); %计算复倒频谱end%绘制输入时程曲线图形subplot(2,1,1);plot(t,x);xlabel('时间(s)');ylabel('幅值');grid on;%绘制倒频谱时程曲线图形subplot(2,1,2);plot(t,y);xlabel('时间(s)');ylabel('幅值');grid on;%打开文件输出倒频谱的数据fid=fopen(fno,'w');for k=1:nfprintf(fid,'%d %f\n',t(k),y(k));endstatus=fclose(fid);程序6-6%地震信号加速度反应谱%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hiddenformat long %%%%%%%%%%%%%%%%%%%%%%fni=input('地震信号加速度反应谱-输入数据文件名:','s'); fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率dp = fscanf(fid,'%f',1); %阻尼比hf = fscanf(fid,'%f',1); %计算最高截止频率df = fscanf(fid,'%f',1); %频率间隔fno = fscanf(fid,'%s',1); %输出数据文件名xg = fscanf(fid,'%f',[1,inf]); %输入地震波数据存成行向量status=fclose(fid);%计算输入数据长度n=length(xg);%计算反应谱长度m=round(hf/df)+1;%定义频率向量f=0:df:(m-1)*df;%定义时间向量t=0:1/sf:(n-1)/sf;%计算加速度反应谱for j=0:m-1%计算单位脉冲反应w=2*pi*df*j;wd=w*sqrt(1-dp*dp);e=exp(-t*w*dp);a=t.*wd;s=sin(a).*((1-2*dp*dp)/(1-dp*dp));c=cos(a).*(2*dp/sqrt(1-dp*dp));h=wd*e.*(s+c)/sf;%用卷积计算加速度反应谱值r(j+1)=max(abs(conv(xg,h)));end%绘制输入地震波加速度曲线subplot(2,1,1);plot(t,xg);xlabel('时间(s)');ylabel('加速度(g)');grid on;%绘制地震波加速度反应谱曲线subplot(2,1,2);plot(f,r);xlabel('频率(Hz)');ylabel('加速度(g)');grid on;%打开文件输出地震波加速度反应谱数据fid=fopen(fno,'w');for k=1:m%每行输出两个实型数据，f为频率，r为反应谱值fprintf(fid,'%f %f\n',f(k),r(k));endstatus=fclose(fid);程序7-1%生成扫频信号%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcclose all %%%%%%%%%%%%%%%%%%%%%%fni=input('生成扫频信号-输入数据文件名:','s');fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率f0 = fscanf(fid,'%f',1); %最小截止频率f1 = fscanf(fid,'%f',1); %最大截止频率t0 = fscanf(fid,'%f',1); %开始时刻(s)t1 = fscanf(fid,'%f',1); %终止时刻(s)am = fscanf(fid,'%f',1); %扫频信号幅值md = fscanf(fid,'%d',1); %扫频方式(1=线性,2=对数,3=指数) fno = fscanf(fid,'%s',1); %输出数据文件名status=fclose(fid);%四舍五入取整求扫频信号开始时刻对应数组元素的下标n0=round(sf*t0)+1;%四舍五入取整求扫频信号终止时刻对应数组元素的下标n1=round(sf*t1)+1;%建立离散时间向量tt=0:1/sf:(n1-1)/sf;%按指定的扫频方式计算随时间变化频率向量switch md%计算扫频方式为线性时随时间变化频率向量case 1ft=(t-t0)./(t1-t0);%计算扫频方式为对数时随时间变化频率向量case 2ft=log(t-t0)./log(t1-t0+1);%计算扫频方式为指数时随时间变化频率向量case 3ft=exp((t-t0)./(t1-t0)-1);otherwiseft=(t-t0)./(t1-t0);end%计算随时间变化的园频率向量w=2*pi*(ft.*(f1-f0)+f0);%建立一个长度与时间长度相同元素全为0的向量y=zeros(1,n1);%计算扫频信号y(n0:n1)=am*sin(w(n0:n1).*t(n0:n1));%绘制算扫频信号随时间变化的曲线图plot(t,y);xlabel('时间(s)');ylabel('幅值');grid on;%打开文件输出扫频信号数据fid=fopen(fno,'w');for k=1:n1%每一行输出两个实型数据，t为时间，y为扫频信号值fprintf(fid,'%f %f\n',t(k),y(k));endstatus=fclose(fid);程序7-2%生成拍波信号%%%%%%%%%%%%%%%%%%%%%%%%%%%clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%fni=input('生成拍波信号-输入数据文件名:','s');fid=fopen(fni,'r');sf = fscanf(fid,'%f',1); %采样频率fp = fscanf(fid,'%f',1); %拍波频率np = fscanf(fid,'%d',1); %拍波个数mp = fscanf(fid,'%d',1); %每个拍波内的周波数ti = fscanf(fid,'%f',1); %拍波之间的时间间隔am = fscanf(fid,'%f',1); %拍波幅值fno = fscanf(fid,'%s',1) %输出数据文件名status=fclose(fid);%建立完整拍波的离散时间向量ts=0:1/sf:(round(sf*mp/fp)-1)/sf;%建立拍波之间的离散时间间隔向量t0=zeros(1,round(ti*sf)+1);%计算拍波圆频率向量w=2*pi*fp;%生成单个拍波的离散信号y1=am*sin(w*ts).*sin(0.5*w*ts/mp);%建立时间间隔和单个拍波信号的组合向量y=[t0,y1];%通过循环建立完整的拍波信号向量for k=1:np-1y=[y,t0,y1];end%建立与整个拍波信号长度相同的离散时间向量t=0:1/sf:(length(y)-1)/sf;%绘制拍波信号随时间变化的曲线图plot(t,y);xlabel('时间(s)');ylabel('幅值');grid on;%打开文件输出生成的拍波信号数据fid=fopen(fno,'w');for k=1:length(y)%每一行输出两个实型数据，t为时间，y为拍波信号值fprintf(fid,'%f %f\n',t(k),y(k));。
振动测试作业Matlab求解程序
第二题求解系统的模态的matlab程序：m=1200j=1100k1=3000k2=1500l1=2l2=1K=[k1+k2 -k1*l1+k2*l2;k1*l1+k2*l2 -k1*l1^2+k2*l2^2] M=[m 0;0 j][eigve,eigva]=eig(inv(M)*K)运行结果为：eigve =0.8492 0.32360.5281 0.9462eigva =1.4178 00 -7.2133龙格库塔法求解程序：function dy=kjt001(t,y)w1=1;w2=1F=1;T=1a1=1200;b1=1100;a2=4500;b2=7500;a3=-4500;b3=-10500;dy = zeros(4,1);dy(1)=y(2);dy(2)=(F*sin(w1*t)-a2*y(2)-a3*y(3))/a1;dy(3)=y(4);dy(4)=(T*sin(w2*t)-b2*y(2)-b3*y(3))/b1;%[t,y] = ode45(@kjt001,[0 30],[0 0 0 0]);%plot(t,y(:,2),'-',t,y(:,4),'-.')求解结果为：第三题中心差分法求解系统方程，matlab运行程序为：clear%输入参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t=linspace(0,10,290);P=[2*sin(2*t); -2*cos(3*t); sin(3.5*t)];k=50;m=2;time=10; %积分时间dim=3; %维数C=zeros(3,3);M=diag([m m m],0); %质量K=k*[2 -1 0;-1 2 -1;0 -1 1]; %刚度%求解特征值问题%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [eigve,eigva]=eig(inv(M)*K);%%%%%%%%%%%%%%CENTRAL INTEGRATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% calculate the time stepnum=1; %时间步长取临界步长的 1/num倍或者num倍check=1; %求解模式：0-时间步长小于临界时间步长；（给模式求解了num=1,2）%1-时间步长大于临界时间步长 (该模式求解了num=2,5,7)if check==0dt=0.1elsedt=num/((eigva(3,3))^0.5*pi);endnn=floor(time/dt+1)-1;%INITIAL CONDITION AND STARTING PROCEDURE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%a=zeros(3,nn);a(:,1)=0;da0=[0 0 0]';dv0=inv(M)*(P(:,1)-K*a(:,1)-C*a(:,1));a_dt=a(:,1)-da0*dt+0.5*dv0*dt^2;%%%%%%%%%% THE CONST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%c0=1/(dt)^2;c1=1/(2*dt);c2=2*c0;c3=1./c2;EM=c0*M+c1*C; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for i=2:nnif i==2Q=P(:,i)-(K-c2*M)*a(:,i-1)-(c0*M-c1*C)*a_dt;elseQ=P(:,i)-(K-c2*M)*a(:,i-1)-(c0*M-c1*C)*a(:,i-2);enda(:,i)=inv(EM)*Q;endfiguretim=[0:dt:(nn-1)*dt];plot(tim,a(1,:),tim,a(2,:),tim,a(3,:))xlabel('t/s （时间步长为临界步长的num倍）')ylabel('distance/m')title('位移时程响应-central difference')legend('mass 1','mass 2','mass 3')运行结果为：。
matlab计算单自由度振动反应的程序
一、引言在工程领域中，单自由度振动系统是一种常见的动力学模型，其在建筑结构、机械设备等领域都有重要的应用。
而计算单自由度振动系统的反应是工程设计和分析中的重要任务之一，对于系统的稳定性、安全性等具有重要意义。
为了快速而准确地计算单自由度振动系统的反应，工程师和研究人员常常使用MATLAB编写程序来实现计算。
二、MATLAB在单自由度振动反应计算中的优势1. 灵活性：MATLAB是一种功能强大的编程语言和工具，可以实现复杂的算法和数学模型，能够满足工程设计和分析中的各种需求。
2. 可视化：MATLAB具有丰富的绘图和可视化功能，可以直观地展示单自由度振动系统的反应结果，使工程师和研究人员更好地理解系统的运动特性。
3. 高效性：MATLAB提供了丰富的计算和求解工具，可以快速而准确地计算单自由度振动系统的反应，节约了工程师和研究人员的时间和精力。
三、MATLAB编写单自由度振动反应计算程序的基本步骤1. 确定系统参数：首先需要确定单自由度振动系统的质量、刚度、阻尼系数等参数，这些参数将影响系统的振动响应。
2. 构建系统模型：根据系统的参数和运动方程，可以利用MATLAB编写对应的单自由度振动系统模型。
3. 求解运动方程：利用MATLAB提供的求解工具，可以求解单自由度振动系统的运动方程，得到系统的振动响应。
4. 可视化结果：最后可以利用MATLAB的绘图和可视化工具，将系统的振动响应以图表的形式展现出来，便于工程师和研究人员对系统的运动特性进行分析和评估。
四、MATLAB编写单自由度振动反应计算程序的示例代码以下是一个简单的MATLAB示例代码，用于计算单自由度振动系统的阻尼比为0.2时的阻尼比对应的阻尼比比，供读者参考。
```matlab定义系统参数m = 1; 质量k = 10; 刚度zeta = 0.2; 阻尼比omega_n = sqrt(k / m); 自然频率计算阻尼比对应的阻尼比比omega_d = omega_n * sqrt(1 - zeta^2);disp(['阻尼比对应的阻尼比比为：', num2str(omega_d /omega_n)]);```通过上述示例代码，可以看出MATLAB的编写方式简单明了，利用MATLAB可以快速计算单自由度振动系统的各种响应参数。
用MATLAB编写程序对机械振动信号进行剖析2
5，李培芳、孙晖、李江主编，信号与系统分析基础。北京：清华大学出版社， 2006.12
6，王宏主编，MATLAB 6.5 及其在信号处理中的应用。 北京：清华大学出版社，2004.10
指导教师签字
说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。
的各种应用函数和一些相关的用户有好操作界面。而工具箱从深度
和广度上大大扩展了 MATLAB 主包的功能和应用领域。随着自身
的不断完善和发展，MATLAB 功能越来越强大，应用也越来越广泛。
2）信号测试技术与分析
随着机械工业不断向自动化、高精度、智能化等方向的发展，
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关，系电通，力1根保过据护管生高线产中0不工资仅艺料可高试以中卷解资配决料置吊试技顶卷术层要是配求指置，机不对组规电在范气进高设行中备继资进电料行保试空护卷载高问与中题带资2负料2，荷试而下卷且高总可中体保资配障料置2试时32卷，3各调需类控要管试在路验最习；大题对限到设度位备内。进来在行确管调保路整机敷使组设其高过在中程正资1常料中工试，况卷要下安加与全强过，看度并25工且52作尽22下可护都能1关可地于以缩管正小路常故高工障中作高资；中料对资试于料卷继试连电卷接保破管护坏口进范处行围理整，高核或中对者资定对料值某试，些卷审异弯核常扁与高度校中固对资定图料盒纸试位，卷置编工.写况保复进护杂行层设自防备动腐与处跨装理接置，地高尤线中其弯资要曲料避半试免径卷错标调误高试高等方中，案资要，料求编试技5写、卷术重电保交要气护底设设装。备备置管4高调、动线中试电作敷资高气，设料中课并技3试资件且、术卷料中拒管试试调绝路包验卷试动敷含方技作设线案术，技槽以来术、及避管系免架统不等启必多动要项方高方案中式；资，对料为整试解套卷决启突高动然中过停语程机文中。电高因气中此课资，件料电中试力管卷高壁电中薄气资、设料接备试口进卷不行保严调护等试装问工置题作调，并试合且技理进术利行，用过要管关求线运电敷行力设高保技中护术资装。料置线试做缆卷到敷技准设术确原指灵则导活：。。在对对分于于线调差盒试动处过保，程护当中装不高置同中高电资中压料资回试料路卷试交技卷叉术调时问试，题技应，术采作是用为指金调发属试电隔人机板员一进，变行需压隔要器开在组处事在理前发；掌生同握内一图部线纸故槽资障内料时，、，强设需电备要回制进路造行须厂外同家部时出电切具源断高高习中中题资资电料料源试试，卷卷线试切缆验除敷报从设告而完与采毕相用，关高要技中进术资行资料检料试查，卷和并主检且要测了保处解护理现装。场置设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
用MATLAB编写程序对机械振动信号进行分析2
燕山大学课程设计（论文）任务书院（系）：电气工程学院基层教学单位：说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。
年月日二、摘要1）MATLAB的简单介绍MATLAB是美国Mathworks公司开发的新一代科学计算软件：MATLAB是英文MATtrix LABoratory(矩阵实验室)的缩写；MATLAB是一个专门为科学计算而设计的可视化计算器。
利用这个计算器中的简单命令，能快速完成其他高级语言只有通过复杂此案出才能实现的数值计算和图形显示。
MATLAB是一种既可交互使用又能解释执行的计算机编程语言。
所谓交互使用，是指用户输入一条语句后立即就能得到该语句的计算结果，而无需像C语言那样首先编写源程序，然后对之进行编译，连接，才能最终形成可执行文件。
MATLAB语言可以用直观的数学表达式来描述问题，从而避开繁琐的底层编程，因此可大大提高工作效率。
MATLAB是解决工程技术问题的技术平台。
利用它能够轻松完成复杂的数值计算，数据分析，符号计算和数据可视化等任务。
MATLAB软件由主包和各类工具箱构成。
其中，主包基本是一个用C/C++等语言编写成的函数库。
该函数库提供矩阵（或数组）的各种算法以及建立在此基础上的各种应用函数和一些相关的用户有好操作界面。
而工具箱从深度和广度上大大扩展了MATLAB主包的功能和应用领域。
随着自身的不断完善和发展，MATLAB功能越来越强大，应用也越来越广泛。
2）信号测试技术与分析随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。
测试信号通常指的是被测对象的运动或状态信息。
测试信号可以用数学表达式描述，也可以用图形、图表等进行描述。
在工程测试中，有的信号可以用数学公式精确描述，而大量的测试信号却只能用数学公式来近似描述。
Matlab振动程序-代码作业
、课题任务要求随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。
本文为机械测试信号的时域和频域分析，其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。
在进行上述分析之前先要对振动信号进行拟合。
机械振动分为确定性振动和随机振动，确定性振动又分为周期振动和非周期振动，周期振动又进一步分为简谐振动和复杂的周期振动。
所以可以根据上述的分类来拟合振动信号。
在设计信号的处理程序时，用MATLAB^的相关函数来对所拟合的振动信号进行时域分析和频域分析，并对绘出的频谱图进行说明。
二、技术路线对机械振动信号的时域和频域采集，根据所拟合的振动信号，选取所需要的时域性能指标和频域分析的性能指标对振动信号进行分析。
其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT 分析、频谱分析、倒谱分析。
现构造一个振动信号（在该程序中以两个衰减振动分量和一个随机数rand之和来拟合振动信号)，再利用MATLAB^的函数mean()、min ()、max()、std ()对离散序列中的平均值、最大值、最小值、标准差等时域性能进行分析，通过调用函数fft (y)；psd (y)；rcep(y)对该振动信号进行频域内的性能分析。
在设计过程中的理论知识有离散傅立叶变换(DFT、功率谱的概念和意义以及倒谱的概念和意义。
①D FT的定义和意义:DFT的定义式为：1 N 1X( n)二—X(k)W N kn =x (n)N k 0DFT的意义：DFT的意义在于它表示信号中的各个频率的分量的振动幅值的大小，亦即该分量对于振动信号影响的大小。
通过DFT的快速算法FFT可以很方便的将振动信号的各个分量的幅值比重计算出来。
matlab 简谐振动 程序
matlab 简谐振动程序以下是一个简单的MATLAB程序，用于模拟简谐振动： matlab.
% 定义参数。
m = 1; % 质量。
k = 1; % 弹簧常数。
x0 = 1; % 初始位移。
v0 = 0; % 初始速度。
t = 0:0.1:10; % 时间范围。
% 计算简谐振动。
omega = sqrt(k/m); % 角频率。
x = x0 cos(omega t) + (v0/omega) sin(omega t); % 位
置随时间变化的函数。
% 绘制图形。
plot(t, x);
xlabel('时间');
ylabel('位移');
title('简谐振动');
这个程序首先定义了简谐振动的参数，包括质量m、弹簧常数k、初始位移x0和初始速度v0。
然后通过设定时间范围t，在每个时间
点上计算简谐振动的位置x。
最后使用plot函数绘制位移随时间变
化的图形。
当你运行这个程序时，你将会得到一个简谐振动的图形，其中
横轴是时间，纵轴是位置。
这个图形展示了简谐振动的特征，包括
振幅、频率和相位。
你可以根据需要修改参数和时间范围，来观察不同条件下的简谐振动情况。
振动力学ode45matlab程序实例
F
解题思路：
设受力大小为：
F=F0 sin ωt 建立平衡位置并受力分析，可得动力学方程如下：
或可写为： 其中：
mẍ+cẋ+kx=F0 sin ωt
ẍ +2ζω0ẋ +ω02x=
F0 m
sin
ωt
可知:
ω0=√mk
；
ζ= c = c
2√km cc
；
x=Xsin(ωt-ϕ)
X=
F0
; λ= ω
√(k-mω2)2+c2ω2
ωn
方程特解为: 方程通解为:
X2=
F0sin(ωt-ϕ) √(k-mω2)2+c2ω2
ห้องสมุดไป่ตู้
x=Aiζω0tsin(ωt+φ)+ F0sin(ωt-ϕ) √(k-mω2)2+c2ω2
即为:
xk F0
=
√[⊢
(mω )2 ]2
1 +
[2+ζ
(ωωn)2]
=
1 √(1-λ2)2+(2ζλ)2
=|G|
定义函数如下：
function y=vibration(t,x) F=0; c=5; k=32; m=2; y=zeros(2,1); y(1)=x(2); y(2)=F/m-c/m*x(2)-(k/m)^2*x(1); end
ode45 处理如下：
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一、课题任务要求
随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。
本文为机械测试信号的时域和频域分析，其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。
在进行上述分析之前先要对振动信号进行拟合。
机械振动分为确定性振动和随机振动，确定性振动又分为周期振动和非周期振动，周期振动又进一步分为简谐振动和复杂的周期振动。
所以可以根据上述的分类来拟合振动信号。
在设计信号的处理程序时，用MATLAB中的相关函数来对所拟合的振动信号进行时域分析和频域分析，并对绘出的频谱图进行说明。
二、技术路线
对机械振动信号的时域和频域采集，根据所拟合的振动信号，选取所需要的时域性能指标和频域分析的性能指标对振动信号进行分析。
其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。
现构造一个振动信号（在该程序中以两个衰减振动分量和一个随
机数rand 之和来拟合振动信号），再利用MATLAB 中的函数mean （）、min （）、max （）、std （）对离散序列中的平均值、最大值、最小值、标准差等时域性能进行分析，通过调用函数fft （y ）；psd （y ）；rcep （y ）对该振动信号进行频域内的性能分析。
在设计过程中的理论知识有离散傅立叶变换（DFT ）、功率谱的概念和意义以及倒谱的概念和意义。
①DFT 的定义和意义：
DFT 的定义式为：
DFT 的意义：DFT 的意义在于它表示信号中的各个频率的分量的振动幅值的大小，亦即该分量对于振动信号影响的大小。
通过DFT 的快速算法FFT 可以很方便的将振动信号的各个分量的幅值比重计算出来。
即可以把信号的主频分量提取出来。
②功率谱的概念和意义：
功率谱的定义式为:若X (Ω)=DFT[x(m)]，x(n)为N 点序列。
则
X *
（Ω） =DFT[x N (-m)] =0k
从而有 DFT[R(M)]=
N 1 DFT[x(m)] DFT[x N (-m)]
即 S ^N x (Ω)=N 1 X(Ω)X *(Ω)=N 1
|X(Ω)|^2
综上所述，先用FFT 求出随机离散序列的DFT ，再计算幅频特性的平方，再除以N ，即得到该随机信号的功率谱估计。
功率谱的意义在于它可以对信号中的周期成分进行分析。
③倒谱的概念与意义：
倒谱可以分析复杂频谱图上的周期成分，分离和提取在密集泛频信号中的成分。
倒频谱对于整个谱的形状不敏感。
三、实现程序
A1=3;A2=4;
f1=100;f2=250;fs=1000;
t=0:1/fs:2;
N=length(t);
X1=A1*exp*t).*sin(2*pi*f1*t);
X2=A2*exp*t).*sin(2*pi*f2*t);
R=rand(1,N);
Y=X1+X2+R;
figure(1);
plot(t,Y);
title('振动信号的波形'); xlabel('时间/秒');
ylabel('幅度');
grid; hold on;%时域分析m= mean(Y);
disp (m);
mi = min(Y); disp (mi); mx = max(Y); disp (mx); st = std(Y);
disp (st);%频域分析
l=length(Y);
r=fft(Y)/l;r=fftshift(r);
f=linspace(-fs/2,fs/2,l);
figure(2);
plot(f,abs(r));
grid; hold on;
figure(3);
psd(Y,2048,1000,kaiser(512,5),0,; figure(4);
yc=rceps(Y);
plot(yc);
四、运行结果
(1)时域分析结果：
序列的平均值为
序列的最小值为
序列的最大值为
序列的标准差为
(2)频域分析结果:
这是所拟合得到的振动信号的图像。
上图为FFT频谱图，从该频谱中可以看到有三个主要高峰值，即在0Hz，100Hz，250Hz处。
在功率谱中可以很明显的看到振动信号中有100Hz和250Hz两个主要的频率。
表明信号中含有这两个频率的周期成分。
五、总结
通过使用MATLAB中的相关函数编写的程序对这一所拟合的振动信号进行了时域分析和频域的分析，得到了关于该振动信号的时域分析结果并绘制出了频域分析图谱。
通过这学期课程的学习，我认识到MATLAB是一款非常强大的软件，应用在各个行业，包括航空航天、人工智能、机械、控制、科学运算等各个领域。
目前我所学习和应用的只是它强大的功能非常小的一部分，之后还要学习MATLAB更多的知识，还需更加努力！。