精华内容
下载资源
问答
  • 思路 粘弹性边界因为能够考虑地基辐射阻尼而使得结构抗震的计算结果更趋于合理,所以在需要考虑结构地基...粘弹性边界是通过在有限元模型的地基边界节点上施加弹簧阻尼器实现的,在ABAQUS的实现有以下几种方法...

     思路

    粘弹性边界因为能够考虑地基辐射阻尼而使得结构抗震的计算结果更趋于合理,所以在需要考虑结构地基相互作用的结构抗震计算时,是较为常用的地基边界处理和地震荷载施加方法。而ABAQUS软件是经常用来进行结构响应分析的有限元软件。下面介绍一种在ABAQUS中实现粘弹性边界及地震荷载施加的方法。

    粘弹性边界是通过在有限元模型的地基边界节点上施加弹簧阻尼器实现的,在ABAQUS中的实现有以下几种方法:第一种,通过ABAQUS自有的弹簧单元spring单元和阻尼单元dashpot实现,具体的单元参数可以参考文献[1],这种较为精确;第二种是通过ABAQUS的UEL子程序实现,可以看下文献[2];还有一种是等效单元替代的方法,就是在地基周围加一层单元,然后设置近似的材料参数,参考文献[3],这一种精度较差,但实现起来较为简单。我采用的是第一种方法,但操作起来较为繁琐,具体程序及过程后面介绍。

    采用粘弹性边界,其配套的地震荷载输入方法就是在已知输入地震位移和速度的情况下,计算各个时刻地基边界各个结点上应当施加的集中力荷载,然后施加荷载,一步一步的进行计算。地震荷载的施加在ABAQUS中也有两种不同的思路,文献[2]中的方法是通过ABAQUS的DLOAD和UTRACLOAD两个子程序实现。DLOAD子程序用于施加边界面的法向荷载,UTRACLOAD用于施加边界面的切向荷载。而文献[1]中则是将边界上每一个节点每个时刻的力都计算出来,然后导入到ABAQUS中作为幅值数据,对每个对应节点施加。

    我最初的想法是两篇文章的思路各取一半,用文献[1]的方法实现粘弹性边界,用文献[2]的方法施加地震荷载。然而尝试了很久,发现这样做的效果并不是太好,可能我编的程序哪儿还是有问题吧。最后放弃了,统一采用文献[1]的方法实现,具体实现采用MATLAB语言生成ABAQUS的input文件,然后将生成的input文件在模型文件的指定位置插入,用ABAQUS运行即可。

    MATLAB程序与运行前的准备工作

    首先需要准备一些必要的数据文件(上图中红色框内的文件),其余黄色框内为模型文件,蓝色框内的文件为程序运行后的生成文件,

    boundary1~5.rpt是从ABAQUS反力文件中提取的反力文件,其值代表某一节点的控制面积,可以通过在地基边界(5个面)施加值为1的压力荷载,即可提取得到这些反力文件;

    coord_point.rpt为5个边界面上节点的坐标文件,提取方法可以在很容易的百度到;

    DIS.txt是三个方向地震波的位移文件;

    VEL.txt是三个方向地震波的速度文件;

    job-996.inp是模型的inp文件;

    Amplitude.inp里面是计算过程中边界节点上随时间要施加的所有集中力荷载,文件较大;

    load.inp是将Amplitude.inp里面的幅值施加到对应节点的荷载命令;

    springs&dashpot.inp是模型地基边界施加的弹簧阻尼器文件;

     

    三个input文件在模型inp文件中的插入位置:

    记住springs&dashpot.inp在Assembly部分,所以搜索到关键字*End Assembly,把springs&dashpot.inp放在*End Assembly之前,Amplitude.inp放在*End Assembly之后,load.inp放在load部分即可,如下所示

    ·································

    ·································

    *Include, Input= springs&dashpot.inp

    *End Assembly

    *Include, Input=Amplitude.inp

    ·································

    ·································

    ** LOADS

    **

    *Include, Input=load.inp

    **

    ** OUTPUT REQUESTS

     

     

    以下为MATLAB程序,记得依据模型修改标红的参数准备并必要的文件

    %%%%

    %%%%-----------------------------说明------------------------------------

    %%%%

    %         1.本程序用于ABAQUS隐式计算的粘弹性边界的inp文件及荷载输入文件的生成

    %         2.需要准备5个边界的节点反力文件(用于节点提取控制面积)、地震动位移和速度文件以及边界节点坐标文件

    %         3.输出为三个文件,springsanddashpot.inp用于施加粘弹性边界;Amplitude.inp为各节点幅值;load.inp为荷载文件

    %         4.首先生成Boundray_point,保存节点编号、节点坐标、控制面积

    %         5.再生成弹簧阻尼器的input文件

    %         6.再生成荷载文件等

    %                                                      by w_tao

    %                                                      2018/10/01

    %%%%---------------------------------------------------------------------

     

    %%%%

    %%%%---------------地基及相关计算参数------------------------------------

    %%%%

    d_time=0.01;

    %地震波的时间间隔,也将是计算步长

    Z0=-600

    %地基底面坐标,Z方向为垂直方向

    Z2=0

    %地基地面坐标

    H=Z2-Z0

    %地基高度

    X0=0

    %地基水平X向坐标

    X1=400

    %地基水平X向坐标

    LLLX=X1-X0

    %地基X向宽度

    Y0=0

    %地基水平Y向坐标

    Y1=400

    %地基水平Y向坐标

    LLLY=Y1-Y0

    %地基Y向宽度

    EEE=4.88e9;

    %地基弹模

    psb=0.22;

    %地基泊松比

    GGG=EEE/2/(1+psb);

    %地基剪切模量

    DENSITY=2000;

    %地基密度

    CP=sqrt((1-psb)*EEE/(1+psb)/DENSITY/(1-2*psb));

    %计算地基中纵波波速

    CS=sqrt(EEE/2/(1+psb)/DENSITY);

    %计算地基中横波波速

    lmt=(CP*CP-2*CS*CS)*DENSITY;

    %拉梅常数中的lmt

    alpha_N=1.33;

    %弹簧阻尼器系数

    alpha_T=0.67;

    %弹簧阻尼器系数

    RRR(1)=LLLX/2;

    RRR(2)=LLLY/2;

    RRR(3)=H;

    %%%%---------------------------------------------------------------------

     

    %%%%

    %%%%--------------------整理节点数据-------------------------------------

    %%%%

    %%%%加载文件数据

    coord_point=load('coord_point.rpt')

    boundary_Z0=load('boundary1.rpt')

    boundary_X1=load('boundary2.rpt')

    boundary_X0=load('boundary3.rpt')

    boundary_Y1=load('boundary4.rpt')

    boundary_Y0=load('boundary5.rpt')

    %  加载地震波数据,默认位移和速度数据是一样长的

    dis=load('DIS.txt');

    vel=load('VEL.txt');

    m=length(dis);

    n(1)=length(boundary_Z0);%底面点数

    n(2)=length(boundary_X1);%X正向面的点数

    n(3)=length(boundary_X0);%X负向面的点数

    n(4)=length(boundary_Y1);%y正向面的点数

    n(5)=length(boundary_Y0);%y负向面的点数

    n(6)=length(coord_point);%总节点数,但不等于5个面的和,因为有重复的点

    error_point_num=0;

     

    %coord_point第1列节点号,第2-4列XYZ坐标,第5列归属的边界面号(1/2/3/4/5),第6列R值,第7列面积

    %循环依据坐标确定归属号和R值

    for i=1:n(6)

    if(coord_point(i,4)==Z0)%底面点

    coord_point(i,5)=1;

    coord_point(i,6)=RRR(3);

    elseif(coord_point(i,2)==X1)%X正向面的点

    coord_point(i,5)=2;

    coord_point(i,6)=RRR(1);

    elseif(coord_point(i,2)==X0)%X负向面的点

    coord_point(i,5)=3;

    coord_point(i,6)=RRR(1);

    elseif(coord_point(i,3)==Y1)%y正向面的点

    coord_point(i,5)=4;

    coord_point(i,6)=RRR(2);

    elseif(coord_point(i,3)==Y0)%y负向面的点

    coord_point(i,5)=5;

    coord_point(i,6)=RRR(2);

    else%理论上没有其他的点了,如果error_point_num非0,说明节点坐标文件不正确

    coord_point(i,5)=0;

    coord_point(i,6)=0;

    error_point_num=error_point_num+1;

    end

    end

    %循环确定节点控制面积

    for i=1:n(6)

    if(coord_point(i,5)==1)

    for j=1:n(1)

    if(coord_point(i,1)==boundary_Z0(j,1))

    coord_point(i,7)=abs(boundary_Z0(j,4));

    end

    end

    elseif (coord_point(i,5)==2)

    for j=1:n(2)

    if(coord_point(i,1)==boundary_X1(j,1))

    coord_point(i,7)=abs(boundary_X1(j,2));

    end

    end

    elseif (coord_point(i,5)==3)

    for j=1:n(3)

    if(coord_point(i,1)==boundary_X0(j,1))

    coord_point(i,7)=abs(boundary_X0(j,2));

    end

    end

    elseif (coord_point(i,5)==4)

    for j=1:n(4)

    if(coord_point(i,1)==boundary_Y1(j,1))

    coord_point(i,7)=abs(boundary_Y1(j,3));

    end

    end

    elseif (coord_point(i,5)==5)

    for j=1:n(5)

    if(coord_point(i,1)==boundary_Y0(j,1))

    coord_point(i,7)=abs(boundary_Y0(j,3));

    end

    end

    else

    error_point_num=error_point_num+1;

    end

    end

    %%%%---------------------------------------------------------------------

     

    %%%%

    %%%%--------------------生成弹簧阻尼器和地震荷载的input文件--------------------------

    %%%%

    fid=fopen('springs&dashpot.inp','w')

    fin_amp=fopen('Amplitude.inp','w')

    fin_load=fopen('load.inp','w')

    k_sum=0;

    m_sum=0;

     

    % 方向向量

    for i=1:3

    dri(i)=i;

    end

    %%%%

    %循环每个节点

    for i=1:n(6)

    %依据节点属于哪个面,确定alpha和C_vel的顺序

    if (coord_point(i,5)==1)%如果是底面(法向方向平行于Z面)的话

    alpha(1)=alpha_T;

    alpha(2)=alpha_T;

    alpha(3)=alpha_N;

    C_vel(1)=CS;

    C_vel(2)=CS;

    C_vel(3)=CP;

    R=RRR(3);

    elseif (coord_point(i,5)==2)%如果法向方向平行于X方向的话

    alpha(1)=alpha_N;

    alpha(2)=alpha_T;

    alpha(3)=alpha_T;

    C_vel(1)=CP;

    C_vel(2)=CS;

    C_vel(3)=CS;

    R=RRR(1);

    elseif (coord_point(i,5)==3)%如果法向方向平行于X方向的话

    alpha(1)=alpha_N;

    alpha(2)=alpha_T;

    alpha(3)=alpha_T;

    C_vel(1)=CP;

    C_vel(2)=CS;

    C_vel(3)=CS;

    R=RRR(1)

    elseif (coord_point(i,5)==4)%如果法向方向平行于Y方向的话

    alpha(1)=alpha_T;

    alpha(2)=alpha_N;

    alpha(3)=alpha_T;

    C_vel(1)=CS;

    C_vel(2)=CP;

    C_vel(3)=CS;

    R=RRR(2);

    elseif (coord_point(i,5)==5)%如果法向方向平行于Y方向的话

    alpha(1)=alpha_T;

    alpha(2)=alpha_N;

    alpha(3)=alpha_T;

    C_vel(1)=CS;

    C_vel(2)=CP;

    C_vel(3)=CS;

    R=RRR(2);

    else

    error_point_num=error_point_num+1;

    end

    %%%%         计算每个节点的弹簧阻尼器参数,并写入文件中

    %%%%         X方向弹簧阻尼器

    k_sum=k_sum+1;

    fprintf(fid,'%s%d\r\n','*Element, type=Spring1, elset=Springs-',k_sum)

    m_sum=m_sum+1;

    fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))

    fprintf(fid,'%s%d\r\n','*Element, type=Dashpot1, elset=Dashpots-',k_sum)

    m_sum=m_sum+1;

    fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))

    fprintf(fid,'%s%d\r\n','*Spring, elset=Springs-',k_sum)

    fprintf(fid,'%d\r\n',dri(1))

    KKK(1)=alpha(1)*GGG/R*coord_point(i,7);

    fprintf(fid,'%d\r\n',KKK(1))

    fprintf(fid,'%s%d\r\n','*Dashpot, elset=Dashpots-',k_sum)

    fprintf(fid,'%d\r\n',dri(1))

    CCC(1)=DENSITY*C_vel(1)*coord_point(i,7);

    fprintf(fid,'%d\r\n',CCC(1))

    %%%%        Y方向弹簧阻尼器

    k_sum=k_sum+1;

    fprintf(fid,'%s%d\r\n','*Element, type=Spring1, elset=Springs-',k_sum)

    m_sum=m_sum+1;

    fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))

    fprintf(fid,'%s%d\r\n','*Element, type=Dashpot1, elset=Dashpots-',k_sum)

    m_sum=m_sum+1;

    fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))

    fprintf(fid,'%s%d\r\n','*Spring, elset=Springs-',k_sum)

    fprintf(fid,'%d\r\n',dri(2))

    KKK(2)=alpha(2)*GGG/R*coord_point(i,7);

    fprintf(fid,'%d\r\n',KKK(2))

    fprintf(fid,'%s%d\r\n','*Dashpot, elset=Dashpots-',k_sum)

    fprintf(fid,'%d\r\n',dri(2))

    CCC(2)=DENSITY*C_vel(2)*coord_point(i,7);

    fprintf(fid,'%d\r\n',CCC(2))

    %%%%        Z方向弹簧阻尼器

    k_sum=k_sum+1;

    fprintf(fid,'%s%d\r\n','*Element, type=Spring1, elset=Springs-',k_sum)

    m_sum=m_sum+1;

    fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))

    fprintf(fid,'%s%d\r\n','*Element, type=Dashpot1, elset=Dashpots-',k_sum)

    m_sum=m_sum+1;

    fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))

    fprintf(fid,'%s%d\r\n','*Spring, elset=Springs-',k_sum)

    fprintf(fid,'%d\r\n',dri(3))

    KKK(3)=alpha(3)*GGG/R*coord_point(i,7);

    fprintf(fid,'%d\r\n',KKK(3))

    fprintf(fid,'%s%d\r\n','*Dashpot, elset=Dashpots-',k_sum)

    fprintf(fid,'%d\r\n',dri(3))

    CCC(3)=DENSITY*C_vel(3)*coord_point(i,7);

    fprintf(fid,'%d\r\n',CCC(3))

     

    %%%%         计算每个节点的集中力时程并写入文件中

    Z1=coord_point(i,4)-Z0;

    %%%%        时间循环

    for j=1:m

    time=(j-1)*d_time; %第一个时刻为0

    UNUM(1)=(time-Z1/CS)/d_time;

    UNUM(2)=(time-(2*H-Z1)/CS)/d_time;

    UNUM(3)=(time-Z1/CS)/d_time;

    UNUM(4)=(time-(2*H-Z1)/CS)/d_time;

    UNUM(5)=(time-Z1/CP)/d_time;

    UNUM(6)=(time-(2*H-Z1)/CP)/d_time;

     

    if(UNUM(1)>0)

    K_NUM(1)=round(UNUM(1))+1;

    U(1)=dis(K_NUM(1),1);

    V(1)=vel(K_NUM(1),1);

    else

    K_NUM(1)=0;

    U(1)=0;

    V(1)=0;

    end

    if(UNUM(2)>0)

    K_NUM(2)=round(UNUM(2))+1;

    U(2)=dis(K_NUM(2),1);

    V(2)=vel(K_NUM(2),1);

    else

    K_NUM(2)=0;

    U(2)=0;

    V(2)=0;

    end

    if(UNUM(3)>0)

    K_NUM(3)=round(UNUM(3))+1;

    U(3)=dis(K_NUM(3),2);

    V(3)=vel(K_NUM(3),2);

    else

    K_NUM(3)=0;

    U(3)=0;

    V(3)=0;

    end

    if(UNUM(4)>0)

    K_NUM(4)=round(UNUM(4))+1;

    U(4)=dis(K_NUM(4),2);

    V(4)=vel(K_NUM(4),2);

    else

    K_NUM(4)=0;

    U(4)=0;

    V(4)=0;

    end

    if(UNUM(5)>0)

    K_NUM(5)=round(UNUM(5))+1;

    U(5)=dis(K_NUM(5),3);

    V(5)=vel(K_NUM(5),3);

    else

    K_NUM(5)=0;

    U(5)=0;

    V(5)=0;

    end

    if(UNUM(6)>0)

    K_NUM(6)=round(UNUM(6))+1;

    U(6)=dis(K_NUM(6),3);

    V(6)=vel(K_NUM(6),3);

    else

    K_NUM(6)=0;

    U(6)=0;

    V(6)=0;

    end

     

    %依据节点属于哪个面,确定集中力荷载中的自由场应力项sigma

    if (coord_point(i,5)==1)%如果是底面(法向方向平行于Z面)的话

    sigma(1)=DENSITY*CS*(V(1)-V(2))*coord_point(i,7);

    sigma(2)=DENSITY*CS*(V(3)-V(4))*coord_point(i,7);

    sigma(3)=DENSITY*CP*(V(5)-V(6))*coord_point(i,7);

    elseif (coord_point(i,5)==2)%如果法向方向平行于X方向的话

    sigma(1)=-lmt/CP*(V(5)-V(6))*coord_point(i,7);

    sigma(2)=0;

    sigma(3)=-DENSITY*CS*(V(1)-V(2))*coord_point(i,7);

    elseif (coord_point(i,5)==3)%如果法向方向平行于X方向的话

    sigma(1)=lmt/CP*(V(5)-V(6))*coord_point(i,7);

    sigma(2)=0;

    sigma(3)=DENSITY*CS*(V(1)-V(2))*coord_point(i,7);

    elseif (coord_point(i,5)==4)%如果法向方向平行于Y方向的话

    sigma(1)=0;

    sigma(2)=-lmt/CP*(V(5)-V(6))*coord_point(i,7);

    sigma(3)=-DENSITY*CS*(V(3)-V(4))*coord_point(i,7);

    elseif (coord_point(i,5)==5)%如果法向方向平行于Y方向的话

    sigma(1)=0;

    sigma(2)=lmt/CP*(V(5)-V(6))*coord_point(i,7);

    sigma(3)=DENSITY*CS*(V(3)-V(4))*coord_point(i,7);

    else

    error_point_num=error_point_num+1;

    end

     

    %节点3个方向的集中力时程

    amp(j,1)=KKK(1)*(U(1)+U(2))+CCC(1)*(V(1)+V(2))+sigma(1);

    amp(j,2)=KKK(2)*(U(3)+U(4))+CCC(2)*(V(3)+V(4))+sigma(2);

    amp(j,3)=KKK(3)*(U(5)+U(6))+CCC(3)*(V(5)+V(6))+sigma(3);

    end

     

    for j=1:3

    %将节点3个方向的集中力时程写入文件中

    fprintf(fin_amp,'%s%d%s%d\r\n','*Amplitude, name=Amp-',coord_point(i,1),'-',dri(j))

    for k=1:m

    fprintf(fin_amp,'%f%s%f\r\n',k*d_time,',   ',amp(k,j))

    end

    end

     

    for j=1:3

    %生成集中力荷载施加的文件

    fprintf(fin_load,'%s%d%s%d\r\n','*Cload, amplitude=Amp-',coord_point(i,1),'-',dri(j))

    fprintf(fin_load,'%s%d%s%d%s\r\n','PART-1-1.',coord_point(i,1),', ',dri(j),', 1.')

    end

     

    end

    status=fclose(fid);

    status=fclose(fin_amp);

    status=fclose(fin_load);

    算例验证

    以文献[2]中的模型作为验证模型,具体参数见下图:

     

    模型如下图所示:

    验证模型400m×400m×600m,紫色为弹簧阻尼器

    计算结果如下所示

    Z方向地基地面点和地基上部点的位移对比

    X方向地基地面点和地基上部点的位移对比

    Y方向地基地面点和地基上部点的位移对比

     

    入射波幅值为0.5m,地面放大2倍接近1m,反射波0.5m

     

    参考文献

    [1]何建涛, 马怀发, 张伯艳,等. 黏弹性人工边界地震动输入方法及实现[J]. 水利学报, 2010, 39(8):960-969.

    [2]苑举卫, 杜成斌, 陈灯红. 基于ABAQUS的三维粘弹性边界单元及地震动输入方法研究[J]. 三峡大学学报(自然科学版), 2010, 32(3):9-13.

    [3]谷音, 刘晶波, 杜义欣. 三维一致粘弹性人工边界及等效粘弹性边界单元[J]. 工程力学, 2007, 24(12):31-37.

    [4]潘坚文. ABAQUS水利工程应用实例教程[M]. 中国建筑工业出版社, 2015.

    转载于:https://www.cnblogs.com/w-tao13614/p/10498697.html

    展开全文
  • 1,将波形图片(截图)保存为test.png或test.jpg,并将图片放于matlab工作目录,如下图示例所指定的目录: 2,新建文件,输入如下程序代码,将文件保存为jpg2data.m(名字可以随便取): 代...

    有如下的波形图,如何从中精确提取出全部的数据:

               42455e886df6e4e2e3cd523cf032b962.png            

    1,将波形图片(截图)保存为test.png或test.jpg,并将图片放于matlab工作目录中,如下图示例所指定的目录中:

               8109fc7e789db760890c19f7f5cc18ae.png            

    2,新建文件,输入如下程序代码,将文件保存为jpg2data.m(名字可以随便取):           572f6f6f76db9d2ed085a719aae4b12a.png      代码:

    % 提取图片中的曲线数据

    clear,clc,close all

    %% 图片与曲线间的定标

    im=imread('test.jpg');%读入图片(替换成需要提取曲线的图片)

    im=rgb2gray(im);%灰度变化

    thresh = graythresh(im);%二值化阈值

    im=im2bw(im,thresh);%二值化

    set(0,'defaultfigurecolor','w')

    imshow(im)%显示图片

    [y,x]=find(im==0);%找出图形中的“黑点”的坐标。该坐标是一维数据。

    y=max(y)-y;%将屏幕坐标转换为右手系笛卡尔坐标

    y=fliplr(y);%fliplr()——左右翻转数组

    plot(x,y,'r.','Markersize', 2);

    disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');

    [Xx,Yy]=ginput(2);%Xx,Yy——指实际坐标框的两个顶点

    min_x=input('最小的x值');%输入x轴最小值,已知图片,可直接赋确定值如min_x=-30;

    max_x=input('最大的x值');%输入x轴最大值,可直接赋确定值如max_x=130;

    min_y=input('最小的y值');%输入y轴最小值,可直接赋确定值如min_y=-2;

    max_y=input('最大的y值');%输入y轴最大值,可直接赋确定值如max_y=5;

    x=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;

    y=(y-Yy(2))*(min_y-max_y)/(Yy(1)-Yy(2))+max_y;

    plot(x,y,'r.','Markersize', 2);

    axis([min_x,max_x,min_y,max_y])%根据输入设置坐标范围

    title('由原图片得到的未处理散点图')

    %% 将散点转换为可用的曲线

    %需处理的问题与解决思路

    %(1)散点图中可能一个x对应好几个y 保留mean()-std()到mean()+std()之间的y值 并取平均处理

    %(2)曲线的最前端和最后段干扰较大 去掉曲线整体的前(如5%)和后5%

    %(3)曲线的最顶端和最底段干扰较大 去掉曲线整体的上10%和下10%

    %参数预设

    rate_x=0.08;%曲线的最前端和最后段删除比例

    rate_y=0.05;%曲线的最顶端和最底段删除比例

    [x_uni,index_x_uni]=unique(x);%找出有多少个不同的x坐标

    x_uni(1:floor(length(x_uni)*rate_x))=[];%除去前rate_x(如5%)的x坐标

    x_uni(floor(length(x_uni)*(1-rate_x)):end)=[];%除去后rate_x的x坐标

    index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];%除去前rate_x的x坐标

    index_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];%除去后rate_x的x坐标

    [mxu,a]=size(x_uni);

    [mx,b]=size(x);

    for ii=1:mxu

       if ii==mxu

           ytemp=y(index_x_uni(ii):mx);

       else

           ytemp=y(index_x_uni(ii):index_x_uni(ii+1));

       end

       %删除方差过大的异常点

       threshold1=mean(ytemp)-std(ytemp);

       threshold2=mean(ytemp)+std(ytemp);

       ytemp(find(ytemp

       ytemp(find(ytemp>threshold2))=[];

       %删除距顶端和底端较近的点

       thresholdy=(max_y-min_y)*rate_y;%y坐标向阈值

       ytemp(find(ytemp>max_y-thresholdy))=[];%删除y轴向距离顶端与底端距离小于rate_y的坐标

       ytemp(find(ytemp

       %剩下的y求均值

       y_uni(ii)=mean(ytemp);

    end

    %此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点

    x_uni(find(isnan(y_uni)))=[];

    y_uni(find(isnan(y_uni)))=[];

    %画图

    % figure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')

    axis([min_x,max_x,min_y,max_y])%根据输入设置坐标范围

    % 将最终提取到的x与y数据保存

    curve_val(1,:)=x_uni';

    curve_val(2,:)=y_uni;

    %% 对提取出的数据进行拟合(按实际情况进行修改)

    [p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);%多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)

    [y_fit,DELTA]=polyval(p,x_uni,s);%求拟合后多项式在x_uni对应的y_fit值

    hold on

    % figure,

    plot(x_uni,y_fit),title('拟合后的曲线')

    axis([min_x,max_x,min_y,max_y])%根据输入设置坐标范围

    % 输出数据到EXCEL

    delete('test.xlsx');

    data = [ x_uni y_fit ];

    [m n]=size(data);

    data_cell=mat2cell(data,ones(m,1),ones(n,1));%matrix转变成cell

    title={'温度差','电压mV'};%添加变量名称,x轴,y轴变量名

    result=[title;data_cell];%归纳变量名称和数据

    S = xlswrite( 'test1.xls',result,'Sheet1' ); %保存数据到excel

    保存文,运行(F5):

               28e45d379ceb767148015036ed8371f8.png            

    输入x最大最小值及y最大最小值:

               467852d73efc7e79b5b6381285b6bee2.png            

               ad92c15bbf8570904dbf20bb6e663ad9.png            

               66482669cae7108aa45137985fd09320.png            

               1986ef39940b8c622d9f64d58b20d83d.png            

    到此就成功用matlab从图片中精确提取出数据!!

    微信关注图中张十三的博客公众号,学习更多技术干货:

               8101eb0664c53dd6a03c8764e0c73cd1.png            

    展开全文
  • m文件基本结构:function[output]=functionname(input);m文件说明:函数定义行,H1备注说明函数用途,可作为索引。。。。 .运算是相同大小的矩阵,对应元素运算。关于/(右除)正常的,\(左除)不正常的。 [C,I]=max(....

    m文件基本结构:function[output]=functionname(input);
    m文件说明:函数定义行,H1备注说明函数用途,可作为索引。。。。

    .运算是相同大小的矩阵,对应元素运算。
    关于/(右除)正常的,\(左除)不正常的。

    [C,I]=max(...),返回最大值C,和位置索引I。

    库函数improd(A,B)=A.*B;

    关于 关系运算法:A==(或者>,<,>=,<=)B;ENTER;比较A B对应元素,若相等则为1,不等为0;

    逻辑流程控制:
    (1)if expression
    statements
    end
    (2)for index1=start1:perlength:end
    statements1
    end
    (3)while expression1(为真时,进入循环)
    statements
    end
    (4)break,遇到break跳出循环,执行循环外的第一句,仅跳出包含它的最内层。
    (5)switch switch_expression
    case case_expression
    statement(s)
    case{case_expression1,case_expression2,...}仅用于unit8,unit6,double类型图像
    statement(s)
    otherwise
    statement(s)
    end

    举例
    a=10;
    b=5;
    while a>0
    a=a-1;
    while b>0
    b=b-1;
    end
    end

    A==imread('C:\MATLAB7\work\DigitalImageProcessing(Gonzalez)\data\d.jpg');%读取图像
    for q=0:5:20
    filename=(sprintf('qq_%3d.jpg',q);
    imwrite(A,filename,'quality',q);
    end

    B=imread('C:\MATLAB7\work\DigitalImageProcessing(Gonzalez)\data\d.jpg');%读取图像
    A=im2bw(B);
    switch newclass
    case 'unit8'
    g=im2unit8(A);
    case 'uint16'
    g=im2unit16(A);
    case 'double'
    g=im2double(A);
    otherwise
    error('unknown or improper image class.')
    end

    关于代码优化:
    matlab为数组运算而设计的编程语言,故而尽量数组化运算。
    一维数组化:x=0:100;二维数组化[C,R]=meshgrid(x,y);

    计算运行时间:tic(开始计时) 程序 t1=toc;(结束计时)。t1就是运行时间。

    关于交互式I/O编程:
    (1)disp(.)用于显示.到屏幕上,而没有ans等 等号左边的内容,可显示任何东西文字字符串等等,文字要加''。
    (2)t=input(' 输入文字(对要输入对象的说明) ','s'),s为字符串或数组。若单纯输入数字,则使用n=str2num(t);
    若输入既有字符又有数字,用[a,b,c,d...]=strread(t,'%f%q%q',',')提取。
    >> t=input('输入名字,年龄,性别','s')
    输入名字,年龄,性别王,23,男
    t =
    王,23,男
    >> [a,b,c]=strread(t,'%c%f%c',',')

    (3)混合数组以及提取:c={'sun',[1 2 3;2 3 4],3};提取c{1}=sun......
    提取另一种方式,将c作为结构体,用c.char_string搜索字符串,用c.matrix搜索矩阵,用c.scalar搜索数字。

    转载于:https://www.cnblogs.com/endlesshunger/p/4438113.html

    展开全文
  • 匿名用户1级2014-05-02 回答在MATLAB中,选择结构可由两种语句来实现。(1) if语句if语句的最简单用法为:if 表达式;程序模块;endif语句的另一种用法为:if 表达式程序模块1else程序模块2end例1 使用if语句判断学生...

    匿名用户

    1级

    2014-05-02 回答

    在MATLAB中,选择结构可由两种语句来实现。

    (1) if语句

    if语句的最简单用法为:

    if 表达式;

    程序模块;

    end

    if语句的另一种用法为:

    if 表达式

    程序模块1

    else

    程序模块2

    end

    例1 使用if语句判断学生的成绩是否及格。

    程序:

    clear

    n=input(’输入n= ’)

    m=60;

    if n<m,

    r=’不及格’

    else

    r=’及格’

    end

    练习一:将例1写入M-文件编辑器,然后在command window 调用这个程序。

    当针对多个条件进行选择时,可以采用下面的格式:

    if 表达式1

    程序模块1

    elseif 表达式2

    程序模块2

    …… ……

    elseif 表达式n

    程序模块n

    else

    程序模块n+1

    end

    例2 将百分之的学生成绩转换为五分制输出。

    程序:

    clear

    n=input(’输入n= ’)

    if n>=90

    chji=’优秀’

    elseif n>=80

    chji=’良好’

    elseif n>=70

    chji=’中等’

    elseif n>=60

    chji=’及格’

    else

    chji=’不及格’

    end

    追问:

    有没有if....

    if....

    else.....

    else....

    这种形式呢,这种类似于c语言的嵌套形式

    追答:

    不会这么弱吧,都告诉你if的语法结构了,嵌套还不容易吗。

    追问:

    嗯,上面的能看懂啊,不过只有上面的两种么,这种呢

    689b4c438705dd3a260b31a5a818cc46.png

    我试了好像不可以啊

    a4ab634a88a7c4d8b346e0eaa33dd898.png

    追答:

    if 表达式;

    程序模块;

    end

    你需要加end,加了吗?

    追问:

    加了,上面是个示意,我想问一下,那样嵌套可以吗,谢谢啦

    追答:

    你得这样写

    if 表达式;

    if 表达式;

    .....

    else

    ....

    end

    end

    追问:

    哦,这样就好啦,谢谢

    追答:

    不客气

    展开全文
  • 如何导入matlab文件格式的数据到python工作环境 首先在导入sciopy.io包 import sciopy.io as sio 找到文件所在的路径 利用loadmat()方法导入数据 dataset = sio.loadmat("/kaggle/input/bearing_location_3...
  • MATLAB编写函数文件的实例

    千次阅读 2019-04-06 10:51:07
    MATLAB编写函数文件的实例 1.在M文件编辑器,编写“ssort”函数的代码 函数代码如下: function out=ssort(a) %ssort程序代码按照升序排列数据 %Define variables: %a input array to sort %ii index variable %...
  • matlab .m文件中定义多个函数

    千次阅读 2011-08-17 01:24:11
    matlab帮助系统的说明: Functions The main difference between a script and a function is thata function accepts input from and returns output
  • 1在MATLAB帮助文件中查找有关title的使用方法,并为y1对应的图形添加标题y1=sin(t);为y2对应图形添加标题y2=e-atcos(3t),其中a根据输入显示具体值(使用num2str函数)。请写出有关指令。t=0:pi/20:4*pi;y1=sin(t);a=...
  • [input_file, Fs] = audioread('movie.AVI');audiowrite('target_file.WAV', input_file, Fs);
  • VS2008C#调用Matlab生成的DLL文件

    万次阅读 热门讨论 2012-02-22 19:23:16
    1、创建一个简单的.m文件 打开Matlab 2009a,新建一个.m文件,输入如下代码: function result=twice(inputvar); result=2*inputvar; 将代码保存为twice.m文件...在Matlab 2009a的Command Window输入deploytool并
  • 将矩阵输出到txt文件中的方法,遍寻网络,始见真经!!! fid=fopen('C:Documentsand Settingscleantotal.ped','wt');%写入文件路径 matrix=input_mattrix%input_matrix为待输出矩阵 [m,n]=size(matrix); for....
  • [MATLAB] 读取ASII文件中的复数数据

    千次阅读 2019-06-20 10:57:36
    写了一段MATLAB代码,实现读入复数数据,并保存下来。复数数据没法直接读取,只能先读取string,再转换。 test.data 文件内容如下: -9+8i -18-2i 6-15i 7-12i -4+25i 8-20i % Read complex data file_name = 'test'...
  • matlab function报错:too many input arguments 使用matlab function的时候,我定义了一个函数ladder,调用的时候,出现报错: 网上的大牛说,产生这个错误的原因是因为函数定义名称和系统定义函数冲突。但是在我...
  • 学习Excel技术,关注微信公众号:excelperfectQ:如下图1所示,一个名为“InputFile.csv”文件,每行有6个数字,每个数字使用空格分隔开。图1现在,我要将以60至69开头的行放置到另一个名为“OutputFile.csv”的文件...
  • matlab脚本文件和函数文件的区别

    千次阅读 2019-05-28 20:34:36
    1.脚本文件(myScript) mynumber= input(‘Enter a number:’) %UNTITLED3 此处显示有关此函数的摘要 % 此处显示详细说明 switch mynumber case -1 disp=’负1′; case 0 disp=’0′; case 1 disp=’正1′; ...
  • 作为function的初学者表示,主程序和子函数要分别保存在两个m文件,并且在同一根目录下~~ %%% 主程序和子函数都出自于blog http://blog.sina.com.cn/s/blog_53f291190100b5ld.html %%% 由于上面的程序有两个问题,...
  • matlab 输出数据以16进制存入文件中

    万次阅读 2018-07-06 13:55:29
    function [] = dec2hex_file(input, row,N, file_path)%input : 输入数据(复数)%row :输入数据长度%N :实部虚部的位宽%file_path :输出数据的保存路径for i=1:row a(i) = real(input(i)); b(i) = imag...
  • matlab中的函数类型

    2020-04-05 16:41:09
    Matlab中的函数可以分为匿名函数、M文件主函数、嵌套函数、子函数、私有函数、和重载函数。 匿名函数通常只是由一句很简单的声明语句组成,使用匿名函数的优点是不需要维护一个M文件。 创建匿名函数的标准格式如下:...
  • 运行NNET工具箱时出错 解决办法: 1、找到subsasgn函数文件 line11 ...3、用上述运行结果的netname=后面的结果,也就是“net”代替找到subsasgn函数文件 line11语句inputname(1) 4、再次运行后没有出...
  • 我用matlab的load 加载vc编辑框中文件路径指向的txt文件,加载不成功。 首先我做的是用按钮将文件路径添加到...其实总的说来就是怎样将vc编辑框中的内容写到matlab中,或者说将编辑框变量与matlab中的变量关联起来。
  • Matlab中滤波器设计-滤波器设计.rar file:///C:/Users/lenovo/Desktop/input.fig利用FDAtool进行滤波器设计的技巧,包括低通、带通和高通滤波器的设计 验证程序 fs=200;%采样频率 t=/fs; s=sin sin sin;%混合sim...
  • 欢迎使用Markdown编辑器写博客 本Markdown编辑器使用StackEdit修改而来,用它写博客,将会带来全新的体验哦: ...导入导出Markdown文件 丰富的快捷键 快捷键 加粗 Ctrl + B 斜体 Ctrl + I...
  • matlab中,如果对函数m文件用中文命名,会出现上述错误,将文件名改为英文后,就没有错了。
  • 问题描述:在打开MATLAB中的.m文件时,出现“ Undefined function or method 'uiopen' for input arguments of type 'char'.”错误
  • MATLAB中要是用.mat文件的数据进行迭代计算,请问下面两种代码写法哪个更快? ``` load Data_1.mat for i =1:300 value = calculation(input1,input2,Data_1) ... end ``` calculation是一个我自定义...
  • 第3章 MATLAB程序设计 3.1 M文件 3.2 程序控制结构 3.3 函数文件 3.4 程序举例 3.5 程序调试;3.1 M文件 ;例3-1 分别建立命令文件... %清除工作空间的变量 f=input'Input Fahrenheit temperature; c=5(f-32)/9 然后在
  • 1 function的命名不能和matlab中的已有命令函数一致 2 每个function都会保存在一个m文件,并且m文件的名称和function的名称一致 3 function里面用到的函数都是固定单一的 4 function内部运行,如果一定要输出,那...
  • m文件的文件名应当是中文的,改成英文字母就可以了
  • matlab中rgb2hsi和hsi2rgb

    2021-04-02 16:32:49
    将此函数保存于某个.m文件中,即可直接调用 上述链接还有hsi2rgb function hsi = rgb2hsi(rgb) %RGB2HSI Converts an RGB image to HSI. % HSI = RGB2HSI(RGB) converts an RGB image to HSI. The input image % ...
  • MatlabGUI调用Simulink编译成可执行的exe文件的解决方法-GUIDE_fig.fig 我(们)曾经N次遇到过这样的提问: “为什么GUI里如果使用sim, simset等函数时,就不能编译成可执行的exe文件发布呢?” 我们也解释...

空空如也

空空如也

1 2 3 4 5 ... 7
收藏数 126
精华内容 50
关键字:

matlab中input文件

matlab 订阅