您的位置：360文档中心› 王济-matlab在振动信号处理中的应用代码

王济-matlab在振动信号处理中的应用代码

合集下载

振动计权加速度级 matlab代码

振动计权加速度级 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振动程序-代码作业

Matlab振动程序-代码作业

一、课题任务要求随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。

本文为机械测试信号的时域和频域分析，其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。

在进行上述分析之前先要对振动信号进行拟合。

机械振动分为确定性振动和随机振动，确定性振动又分为周期振动和非周期振动，周期振动又进一步分为简谐振动和复杂的周期振动。

所以可以根据上述的分类来拟合振动信号。

在设计信号的处理程序时，用MATLAB中的相关函数来对所拟合的振动信号进行时域分析和频域分析，并对绘出的频谱图进行说明。

二、技术路线对机械振动信号的时域和频域采集，根据所拟合的振动信号，选取所需要的时域性能指标和频域分析的性能指标对振动信号进行分析。

其中时域分析包括对信号最大值、最小值、中值、方差的分析，频域分析包括FFT分析、频谱分析、倒谱分析。

现构造一个振动信号（在该程序中以两个衰减振动分量和一个随机数rand 之和来拟合振动信号），再利用MATLAB 中的函数mean （）、min （）、max （）、std （）对离散序列中的平均值、最大值、最小值、标准差等时域性能进行分析，通过调用函数fft （y ）；psd （y ）；rcep （y ）对该振动信号进行频域内的性能分析。

在设计过程中的理论知识有离散傅立叶变换（DFT ）、功率谱的概念和意义以及倒谱的概念和意义。

①DFT 的定义和意义：DFT 的定义式为：DFT 的意义：DFT 的意义在于它表示信号中的各个频率的分量的振动幅值的大小，亦即该分量对于振动信号影响的大小。

通过DFT 的快速算法FFT 可以很方便的将振动信号的各个分量的幅值比重计算出来。

即可以把信号的主频分量提取出来。

两自由度振动系统matlab代码

两自由度振动系统matlab代码

两自由度振动系统是一个经典的物理学问题，它描述了两个质点在受到某种力的作用下进行复杂振动的情况。

在工程和物理学领域，研究两自由度振动系统可以帮助我们更好地理解和控制实际系统中的振动现象，比如建筑结构、机械系统和电路等。

而利用Matlab来模拟和分析两自由度振动系统，可以更直观地展现系统的运动规律和特性。

让我们来定义一个经典的两自由度振动系统。

假设有两个质点$m_1$和$m_2$，它们分别位于$x_1$和$x_2$的位置上，通过弹簧和阻尼器相连接。

系统中存在的力包括弹簧力、阻尼力和外力，而质点的运动受到这些力的综合作用而产生振动。

我们可以利用Newton力学原理建立方程来描述系统的运动。

对于一个给定的两自由度振动系统，我们可以通过编写Matlab代码来模拟其振动过程。

在Matlab中，我们可以利用ODE求解器来求解系统的微分方程，从而得到系统的位置随时间的变化。

通过改变系统的参数和初始条件，我们可以观察到系统的不同运动模式和特性。

另外，Matlab还提供了丰富的绘图函数，可以直观地展现系统振动的过程和特性。

当我们编写Matlab代码来模拟两自由度振动系统时，我们需要注意几个关键的步骤。

我们需要建立系统的微分方程模型，这涉及到物理建模和方程的推导。

我们需要选择合适的数值方法和求解器来求解微分方程，比如常见的四阶Runge-Kutta法。

我们可以用Matlab来进行数值模拟，并利用绘图函数来展现系统的振动规律。

我们可以对模拟结果进行分析，比如计算系统的频率响应、模态分析和阻尼比等特性。

对于我个人的观点和理解，两自由度振动系统的分析和模拟在工程和科学研究中具有重要意义。

通过对系统振动特性的研究，我们可以更好地设计和控制实际系统，从而提高系统的性能和稳定性。

而利用Matlab来模拟系统的振动行为，则可以帮助我们更直观地理解系统的动态特性，为工程实践提供有力的支持。

两自由度振动系统是一个具有重要应用价值的物理问题，通过Matlab 代码的编写和模拟分析，我们可以更好地理解和掌握系统的振动特性。

振动信号变分模态分解matlab

振动信号变分模态分解matlab

振动信号变分模态分解matlab振动信号变分模态分解(Matlab)在工程和物理学领域中，振动信号是一种广泛存在的信号类型。

通过对振动信号进行分解可以更好地理解信号，分析信号，并进一步实现信号的相关处理。

其中，变分模态分解（VMD）是一种先进的信号分解技术，可用于将振动信号分解成多个子信号并分析每个子信号。

本文将介绍如何使用Matlab进行VMD分解。

1. 安装VMD软件包VMD分解技术已经被打包为Matlab的一个工具箱，可以通过下载和安装该工具箱来安装VMD软件包。

要开始安装，请浏览VMD软件包的官方网站，下载适用于您的Matlab版本的该软件包。

下载完成后，运行Matlab并将软件包添加到路径中。

接下来，您可以使用Matlab中的VMD函数进行信号分解。

2. VMD分解信号首先，您需要准备一个要分解的信号。

一个简单的方法就是使用Matlab中的随机函数生成一些样本数据进行测试。

示例代码：% Generate sample signalx = randn(1, 1000);% Perform VMD decompositionalpha = 2000;tau = 0;K = 10;DC = 1;init = 1;tol = 1e-7;omega = 0;verbose = 1;[u, ~] = VMD(x, alpha, tau, K, DC, init, tol, omega, verbose);上述代码使用了VMD函数对信号进行分解，其中参数alpha表示VMD算法的控制参数，tau表示信号的初态时间，K表示分解生成的子信号数，DC表示是否保留分解后信号的直流分量，init表示算法的初始条件，tol表示解决非线性优化的容差，omega表示算法的微调选项， verbose表示是否输出调试信息。

3. 分析VMD分解结果一旦分解完成，您可以分析每个子信号以获取信号的相关特征。

例如，您可以使用Matlab中的快速傅里叶变换（FFT）来计算每个子信号的频谱。

matlab振动算法

matlab振动算法

matlab振动算法
MATLAB是一种用于数学建模、仿真和数据分析的强大工具，它提供了许多用于处理振动问题的算法和工具。

在MATLAB中，振动问题通常涉及到求解微分方程、频率分析、模态分析等内容。

以下是一些在MATLAB中处理振动问题常用的算法和工具：
1. 求解微分方程，MATLAB提供了强大的微分方程求解器，如ode45、ode23等，用于求解振动系统的运动方程。

用户可以通过编写自定义的微分方程函数来描述振动系统的运动规律，并利用求解器得到系统的解析解或数值解。

2. 频率分析，MATLAB中的信号处理工具箱提供了丰富的频谱分析函数，如fft、pwelch等，用于分析振动信号的频谱特性。

用户可以通过这些函数对振动信号进行频谱分析，了解系统的频率响应特性。

3. 模态分析，MATLAB中的模态分析工具箱提供了用于计算结构模态参数的函数，如modeShape、naturalFrequency等。

用户可以利用这些函数计算振动系统的模态形状和固有频率，从而了解系统的振动特性。

4. 有限元分析，MATLAB中的有限元分析工具箱提供了用于建立和求解有限元模型的函数，用户可以利用这些函数对复杂结构的振动特性进行分析和预测。

总之，MATLAB提供了丰富的算法和工具，用于处理各种振动问题，用户可以根据具体的振动分析需求选择合适的算法和工具进行使用。

希望以上信息能够帮助到你。

matlab在振动信号处理中的应用代码

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在机械振动信号中的应用

MATLAB在机械振动信号中的应用申振（山东理工大学交通与车辆工程学院）摘要：综述了现代信号分析处理理论、方法如时域分析（包括时域参数识别、相关分析等）、频域分析（包括傅立叶变换、功率谱分解等),并结合MATLAB中的相关函数来对所拟合的振动信号进行时域分析和频域分析,并对绘出的频谱图进行说明.关键词:时域分析频域分析 MATLAB信号是信息的载体，采用合适的信号分析处理方法以获取隐藏于传感观测信号中的重要信息(包括时域与频域信息等），对于许多工程应用领域均具有重要意义。

对获取振动噪声信号的分析处理，是进行状态监测、故障诊断、质量检查、源识别、机器产品的动态性能测试与优化设计等工作的重要环节，它可以预先发现机械部件的磨损和缺陷等故障，从而可以提高产品的质量，降低维护费用。

随着测试技术的迅速发展,各种信号分析方法也随之涌现,并广泛应用在各个领域［1]。

时域描述简单直观，只能反映信号的幅值随时间的变化,而不能明确的揭示信号随时间的变化关系。

为了研究信号的频率组成和各频率成分的幅值大小、相位关系，应对信号进行频谱分析,即把时域信号通过适当的数学方法处理变成频率f(或角频率 ）为独立变量，相应的幅值或相位为因变量的频域描述。

频域分析法将时域分析法中的微分或差分方程转换为代数方程，有利于问题的分析[2］.MATLAB是MathWorks公司于1982年推出的一种功能强大、效率高、交互性好的数值计算和可视化计算机高级语言,它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体,形成了一个极其方便、用户界面良好的操作环境。

随着其自身版本的不断提高，MATLAB的功能越来越强大,应用范围也越来越广，如广泛应用于信号处理、数字图像处理、仿真、自动化控制、小波分析及神经网络等领域[3].本文主要运用了MATLAB R2014a对机械振动信号进行分析.分析过程包括时域分析和频域分析两大部分，时域分析的指标包括随机信号的均值、方差以及均方值。

振动测试作业Matlab求解程序

振动测试作业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进行地震信号处理和振动分析引言地震信号处理和振动分析是地球科学中非常重要的研究领域。

随着计算机技术的发展，利用计算机编程语言进行数据处理和分析已成为地震学和工程地震学的常用方法。

在本文中，将介绍如何使用Matlab进行地震信号处理和振动分析。

一、Matlab简介Matlab是一种强大的科学计算软件，广泛应用于各个领域，包括地震学。

它具有丰富的函数库和图形化界面，提供了各种数据处理和分析工具，非常适合用于地震信号处理和振动分析。

二、地震信号处理在地震学中，地震信号通常是通过地震仪器记录的地震波形数据。

地震信号处理的目标是从原始数据中提取地震波形特征，如到达时间、波形振幅、频率等。

Matlab提供了多种处理方法和函数，方便地进行地震信号的滤波、增益校正、相位校正等操作。

1. 地震信号滤波地震信号通常包含各种频率分量，包括低频、中频和高频分量。

为了分析和识别地震事件，需要对地震信号进行滤波，去除干扰信号并突出地震信号的特征。

Matlab提供了多种滤波函数，如低通滤波、高通滤波、带通滤波等，可以根据需求选择适合的滤波方法。

2. 特征提取地震波形中的各种特征包含了地震事件的重要信息，如震源距离、震级、震中位置等。

Matlab提供了多种特征提取方法和函数，可以从地震波形数据中提取到达时间、波形振幅、频率等特征，并帮助地震学家进行地震事件的分析和研究。

三、振动分析振动分析是工程地震学中的一项重要任务，旨在研究结构在地震或其他振动作用下的响应和受力。

通过对结构振动的分析，可以评估结构的安全性并制定相应的安全标准。

Matlab提供了多种振动分析方法和函数，方便地进行结构的模态分析、响应谱分析等。

1. 结构模态分析结构的模态分析是指在预定边界条件下，确定结构的固有频率、振型和振动模态。

利用Matlab可以进行结构的模态分析，并绘制模态图，有助于工程师评估结构的动力性能和稳定性。

2. 结构响应谱分析结构响应谱分析是指通过计算结构在地震作用下的响应谱，来评估结构的受力性能和安全性。

matlab在振动信号处理中的应用

matlab在振动信号处理中的应用

matlab在振动信号处理中的应用
MATLAB在振动信号处理中得到了广泛的应用。

首先，可以利用它对振动信号进行快速可视化分析，比如选择正确的显示器，看看振动信号的谱图和功率谱等。

其次，它可以为振动信号建立预测模型，采用最小二乘法、贝叶斯估计等技术，从而更好地了解振动的实时行为及其原因。

此外，MATLAB还可以用来分析振动信号，识别各种异常振动，特别是在维护和诊断振动系统中更为重要。

另外，MATLAB中也有相关的机器学习算法，如聚类、回归分析等，可以用来进一步处理振动信号数据，从而更好地理解振动的动态行为及其特征。

如何使用MATLAB进行信号处理的基本操作

如何使用MATLAB进行信号处理的基本操作

如何使用MATLAB进行信号处理的基本操作MATLAB是一种功能强大的数学计算软件，用于信号处理的基本操作。

信号处理是一种涉及测量、分析和操纵信号的技术，广泛应用于通信、音频处理、医学成像等领域。

本文将介绍如何使用MATLAB进行信号处理的基本操作，包括信号生成、采样和重构、频谱分析、滤波和相关性分析。

第一章：信号生成信号生成是信号处理的首要步骤，涉及到产生原始信号以及添加噪声等。

在MATLAB中，可以使用函数来生成各种类型的信号，如正弦信号、方波信号、三角波信号等。

例如，可以使用sin函数生成一个正弦信号：```matlabt = 0:0.1:10; % 时间从0到10，步长为0.1f = 1; % 频率为1Hzx = sin(2*pi*f*t); % 生成正弦信号```此外，可以使用randn函数生成高斯白噪声信号：```matlabn = length(t); % 信号长度noise = 0.1*randn(1,n); % 生成标准差为0.1的高斯白噪声信号```第二章：采样和重构采样是将连续时间信号转换为离散时间信号的过程，重构则是将离散时间信号再转换为连续时间信号。

在MATLAB中，可以使用采样函数进行采样和重构操作。

采样函数包括：A/D（模拟到数字）和D/A（数字到模拟）转换。

例如，可以使用函数`resample`进行信号的采样和重构：```matlabFs = 100; % 采样频率为100HzTs = 1/Fs; % 采样时间间隔t = 0:Ts:1; % 采样时间段为1秒x = cos(2*pi*10*t); % 原始信号，频率为10Hz的余弦信号y = resample(x, 2, 1); % 按2倍采样重构信号```可以通过观察原始信号和重构信号的波形来验证采样和重构的效果。

第三章：频谱分析频谱分析是信号处理中重要的环节，可以用于分析信号的频率成分。

在MATLAB中，可以使用快速傅里叶变换（FFT）函数进行频谱分析。

试析MATLAB在汽车振动分析及控制中的运用

试析MATLAB在汽车振动分析及控制中的运用

试析MATLAB在汽车振动分析及控制中的运用摘要：随着汽车振动会对汽车造成一定的危害，影响汽车的正常运行，因此要对汽车振动的原因进行分析，并提出合理的解决对策。

而MATLAB在汽车有巨大的作用，既能够应用与汽车的振动分析和控制中，还能够预测到汽车的反映，在实际中有较高的使用价值。

关键词：汽车振动；MATLAB；控制；应用在实际中汽车振动具有巨大的危害，因此一方面为了减少振动对汽车的零部件和使用性能造成的损害，另一方面还要针对汽车振动进行研究，尽量把汽车振动运用到汽车设计的服务中来，同时还要尽量利用振动的原理来制造推动机，减少汽车部件工作的强度，从而来提高工作效率。

1.汽车振动产生的原因分析机械的振动实际上是一种形式特殊的运动，构成振动系统的主要因素有阻尼、弹性、质量以及激励。

汽车本身就是具有阻尼、弹簧以及质量的系统，在汽车内部各部分具有不同的频率，因而汽车在平时的行驶中就往往会因为路面不平坦、运动方向以及车速的改变、发动机车轮与传送系统不平衡、汽车齿轮之间的冲击等各种内部与外部的激振作用造成汽车的局部或者是整体发生强烈的振动。

这种振动现象就会使汽车在运动时由于动力性能没有得到充分的发挥就容易出现一些问题。

另外，振动还会对汽车操作的平顺性与稳定性以及汽车的通过性造成不利的影响，还会使乘员在乘坐时产生不舒服的感觉，严重情况下还会对汽车的零部件造成损坏，缩短汽车正常的使用寿命。

2. 汽车振动的理论分析2.1.汽车振动的模型在实际中当一个振动系统比较复杂的时候，而建立的相应的模型也比较复杂，就越接近于真实的情况，相应的模拟情况就越真实，但是这却增加了对系统分析的难度，所以在建立振动的力学模型时就需要在逼真模拟与系统分析之间找出一个平衡点。

在进行汽车振动分析时需要把握的一个关键因素就是在实际中要根据研究的要求和研究内容，把研究的对象和外界的作用力进行简化，同时还要保证简化的模型能够与原来的系统模型在动态分析方面进行等效对比。

加速度转换成位移的matlab代码及说明

加速度转换成位移的matlab代码及说明

加速度转换成位移的matlab代码及说明由测量的加速度离散数据数据转化成位移数据一般不直接在时域进行积分处理，而是由时域转换成频域在频域中进行二次积分再转化到时域中得到位移结果。

相关matlab处理程序方法参考王济老师的matlab在振动信号处理中的应用中的程序如下：%频域积分%%%%%%%%%%%%%%%%%%%%%%%%%%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);%最大截止频率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 == 2y=-a; %进行二次积分的相位变换elsea1=imag(a); a2=real(a); y=a1-a2*i; %进行一次积分的相位变换enda=zeros(1,nfft);%消除指定正频带外的频率成分a(ni:na)=y(ni:na);%消除指定负频带外的频率成分a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);y=ifft(a,nfft); %IFFT变换%取逆变换的实部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:n, fprintf(fid,'%f \n',y(k)); endstatus=fclose(fid);程序使用说明：运行该程序会显示频域积分-输入数据文件名：在之后写出输入数据文件的完整文件名（包括扩展名，例如data.txt）要求数据文件必须在matlab工作空间路径下，该数据文件必须满足如下格式要求：采样频率下限截止频率上限截止频率单位变换系数积分阶数（加速度变位移为2，变速度为1）时间（s）加速度（g）位移（mm）out.txt具体加速度数据（回车或空格隔开）下图为书上给出的例子：满足这样的格式的数据文件才能正常处理，注意单位变换系数和坐标有关，如果按书上的例子加速度以g为单位位移以mm为单位，单位变换系数应该为9.8*1000=9800；而用米每二次方秒为加速度单位以mm为位移单位的话就是1*1000=1000；如果m为位移单位的话那单位变换系数就应该为1了。

Matlab在信号处理中的应用范例

Matlab在信号处理中的应用范例
X1
1.2 MATLAB特色举例
· >> [x,y] = meshgrid(-3:1/8:3); z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2)- 10*(x/5 - x.^3 -
y.^5)... .*exp(-x.^2-y.^2)- 1/3*exp(-(x+1).^2 - y.^2); surf(x,y,z), shading interp; colorbar
X1
1.6 MATLAB 函数编写方法与应用
· 所 谓 MATLAB 程序，大致分为两类：M 脚本文件(M-script) 和 M 函数(M-function), 它们均是普通的ASCII 码构成的文件。 M 脚本文件中包含一族由MATLAB 语言所支持的语句，它类似 于 DOS 下的批处理文件，它的执行方式很简单，用户只需在 MATLAB 的提示符>> 下键入该M 文件的文件名，这样 MATLAB 就会自动执行该M 文件中的各条语句，并将结果直接返回到 MATLAB 的工作空间。M 函数格式是MATLAB 程序设 计的主流，一般情况下， 不建议您使用M 脚本文件格式编程。 · MATLAB 的M 函数是由function 语句引导的，其基本格式 如下： · function [返回变量列表] = 函数名(输入变量列表) 注释说明语句段, 由% 引导 输入、返回变量格式的检测 函数体语句
MATLAB 还支持其他运算，如取整、求余数等。可以使用rond（）, fix(), rem() 等来实现。
X1
1.5 MATLAB 的语句流程与控制
· 作为一种常用的编程语言，MATLAB 支持各种流程控制结构，如循环结构、条件 转移结构、客观结构等另外MATLAB 还支持一种新的结构--- 试探结构。 · 循环语句有两种结构： for ... end 结构和while ... end 结构。 这两种语句结构不 完全相同，各有各的特色。for ... end 语句通常的调用格式为： ·for 循环变量=s1:s3:s2循

第八讲 MATLAB在信号处理中的应用(一)

第八讲 MATLAB在信号处理中的应用(一)

1, n n0 u (n n0 ) 0, n n0
利用MATLAB函数文件表示为： function [x,n]=stepseq(n0,n1,n2) n=n1:n2; x=[(n-n0)>=0]; >> [x,n]=stepseq(3,0,5) >>stem(n,x)

其中 x ( n ) 看成 x(n) 以N为周期的周期延拓，即

~ x(n),0 n N 1 x ( n) 0, 其他n ~ 另一种表示为 x(n) x((n)) N 其中 x((n)) N 表示为 (n模N

即n对N取余数)例如N=16,n=25时 因为n=25=1*16+9，故

例如要实现在区间
上的样本序列，可表示为 1, n n0 ( n n0 ) 0, n n0
n1 n0 n2
用函数文件描述样本序列： function[x,n]=imseq(n0,n1,n2) n=n1:n2; x=[n-n0==0]； >>[x,n]=imseq(3,0,5) >>stem(n,x)
2.信号相乘 function [y,k]=sigmult(f1,k1,f2,k2) k=min(min(k1),min(k2)):max(max(k1),max(k2)); y1=zeros(1,length(k)); y2=y1; y1(find((k>=min(k1))&(k<=max(k1))==1))=f1; %find函数表示寻找非零元素的下标 y2(find((k>=min(k2))&(k<=max(k2))==1))=f2; y=y1.*y2;

王济matlab在振动信号处理中的应用代码

王济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**sf-df);%建立负的离散圆频率向量w2=2*pi**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**sf-df);w2=2*pi**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 %%%%%%%%%%%%%%%%%%%%%%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+*fs/np;%取大于nt且与nt最接近的2的整数次方为FFT长度nf=2^nextpow2(nt);%确定细化带宽的数据长度na=round*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 = [ ];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%地震信号加速度反应谱%%%%%%%%%%%%%%%%%%%%%%%%%%% clearclcformat 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*w*ts/mp);%建立时间间隔和单个拍波信号的组合向量y=[t0,y1];%通过循环建立完整的拍波信号向量for k=1:np-1y=[y,t0,y1];。

MATLAB在振动台试验数据处理中的应用

MATLAB在振动台试验数据处理中的应用
ｍ ｅｈ ｄ ｉ ｖ ｒ ａ ｉｆｃｏ ｙ，ａ ｄ ｉ ｃ ｎ ｂ ｎ ｒ ｄ ｃ ｄ ｔ ｔ ｅ ａ ａｐ ｏ ｅｓｎ ｏ ｘｐ ｒｍ ｅ ｔ ． ｔ ｏ ｓ ｅ ｙ ｓ ｔｓａ ｔ ｒ ｎ ｔ ａ ｅｉｔｏ ｕ ｅ ｏ ｏ ｈ ｒｄ ｔ ｒ ｃ ｓｉｇ ｆ ｒｅ ｅ ｉ ｎ ｓ
摘 要 ：利 用 ＭＡＴ Ａ Ｌ Ｂ语 言编 写 了 用于 振 动 台试 验 数据 处理 的程 序 ， 据模 型试 验 所 得 的 加 速 度 波 形 根
可 计 算 得 出相 应 的积 分 速 度 及 位 移 。 通 过 对 某 换 流 站 大 跨 度 结 构 振 动 台试 验 数 据 的 处 理 表 明 ， 方 法 该 所 得 结 果 正 确 且 通 用性 强 ， 可推 广应 用到 其 它试 验 数 据 处 理 中。 关键词 ： ＭＡ Ｌ Ｂ； 动 台试 验 ： 据 处 理 Ｔ Ａ 振 数 中 图分 类 号 ： ３ ５ ８ Ｐ １ ． 文 献标 识码 ： Ａ 文章 编 号 ： ６ ２ １４ （０ ８０ — ０ ２ — ０ １７ — １４ ２ ０ ）１ １ １ ２
果。
２ 用 ＭＡ Ａ ＴＬ Ｂ求 积分位移
在振 动 台试 验 中 ， 测量 加速 度 比 测 量 位 移 相 对 来 说 方 便 些 ， 据 模 型 加速 度 波 形 ， 过 数 值 积 分 求 得 速 度 和 位 移 波 根 通
用 由拟 建 厂 房 当 地 地 震 局 为 该 工 程 场 地 提 供 的 场 地 波 。 Ｅ ｌ
—
批 功 能 强 大 的核 心 内部 函数 和工 具 箱 函 数 ， 需 要 高 深 的 不
Ｃ ｎｒ 和 Ｔ ｆ 波所 含 频 率 成 分 丰 富 ， 地 条 件接 近 本 工 ｅｔ ｏ波 ａ ｔ 场

MATLAB在车辆振动分析中的应用研究

MATLAB在车辆振动分析中的应用研究

MATLAB在车辆振动分析中的应用研究MATLAB是一种强大的工具，用于进行车辆振动分析。

车辆振动分析是机械工程中的一个重要领域，也是汽车工程中的重要分支，对汽车的性能和安全性有很大的影响。

在这篇文章中，我们将探讨如何在MATLAB中进行车辆振动分析以及一些应用案例。

MATLAB是一种非常流行的数学软件，可以用于处理各种数学和工程问题。

MATLAB提供了许多功能，用于实现车辆振动分析。

以下是其中一些常用功能：1. 振动信号处理：MATLAB提供了许多振动信号处理工具，用于对车辆振动信号进行分析和处理。

这些工具包括傅里叶变换、小波变换、滤波器等。

2. 模态分析：MATLAB可以用于进行车辆结构的模态分析，以确定车辆的固有频率和振型。

这对于设计和优化车辆的结构非常重要。

3. 车辆动力学分析：MATLAB可以用于对车辆的动力学参数进行分析和计算，如车辆的加速度、速度、位移、角度等参数。

4. 有限元分析：MATLAB可以与有限元分析软件配合使用，进行车辆结构的有限元分析。

这有助于确定车辆的结构是否足够强壮，以应对各种振动载荷。

1. 车辆悬挂系统分析：车辆悬挂系统对车辆的振动特性有很大的影响。

使用MATLAB 可以对车辆悬挂系统进行动态分析，以计算车辆的自然频率、模态形状和阻尼比等参数。

这些参数可以用于优化车辆悬挂系统的设计，提高车辆的舒适性和稳定性。

3. 车辆噪声控制分析：车辆噪声是一种常见的底盘振动问题。

MATLAB可以用于对车辆噪声进行预测和控制。

使用MATLAB可以对车辆噪声进行频谱分析，以确定噪声的频率特性。

然后，可以使用MATLAB中的振动控制算法来减小车辆噪声。

4. 车辆结构损伤诊断：车辆经过长期使用后，车辆结构上可能出现疲劳、裂纹和变形等损伤。

MATLAB可以用于对车辆结构损伤进行诊断和预测。

利用MATLAB可以对车辆振动信号进行分析和处理，以检测车辆结构的损伤状态。

相关主题
  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

程序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));。

相关文档
最新文档