精华内容
下载资源
问答
  • MATLAB曲线拟合工具箱(cftool)介绍(完结)

    千次阅读 多人点赞 2021-07-23 18:47:24
    本文通过实例对MATLAB曲线拟合工具箱进行详细讲解,帮助大家更容易理解曲线拟合工具箱(cftool)。 目录1.实例介绍2. 进入系统辨识工具箱界面3. 加载数据4. 加载数据5. 选择拟合曲线的类型 1.实例介绍 已知 x = [0 ...

    本文通过实例对MATLAB曲线拟合工具箱进行详细讲解,帮助大家更容易理解曲线拟合工具箱(cftool)。

    1.实例介绍

    已知
    x = [0 0.2 0.50.8 0.9 1.3 1.4 1.9 2.1 2.2 2.5 2.6 2.9 3.0];
    y = [1.27792.1596 2.7311 2.5974 2.4068 1.6215 1.4178 0.9955 0.9666 0.8837 0.9639 1.00311.1233 1.1583];

    并且根据某种物理或数学关系确定y=f(x)的表达形式,并求出拟合结果对应的系数。

    2. 进入曲线拟合工具箱界面

    两种方法,第一种:
    打开app栏的曲线拟合工具箱(Curve Fitting),
    在这里插入图片描述
    第二种,直接在命令行窗口输入“cftool”:
    在这里插入图片描述
    进入界面后,弹出如下窗口:
    在这里插入图片描述

    3. 加载数据

    新建一个.m文件,并写入如下代码:

    clc;clear;
    x = [0 0.2 0.5 0.8 0.9 1.3 1.4 1.9 2.1 2.2 2.5 2.6 2.9 3.0];
    y = [1.2779 2.1596 2.7311 2.5974 2.4068 1.6215 1.4178 0.9955 0.9666 0.8837 0.9639 1.0031 1.1233 1.1583];
    

    作为要拟合曲线所需要的数据。

    4. 加载数据

    在上述窗口中选中相应数据和拟合选项:在这里插入图片描述

    5. 选择拟合曲线的类型

    通过下拉菜单选择拟合曲线的类型,
    在这里插入图片描述
    工具箱提供的拟合类型有:

    • Custom Equations:用户自定义的函数类型;
    • Exponential:指数逼近,有2种类型, a ∗ e x p ( b ∗ x ) a*exp(b*x) aexp(bx) a ∗ e x p ( b ∗ x ) + c ∗ e x p ( d ∗ x ) a*exp(b*x) + c*exp(d*x) aexp(bx)+cexp(dx)
    • Fourier:傅立叶逼近,有7种类型,基础型是 a 0 + a 1 ∗ c o s ( x ∗ w ) + b 1 ∗ s i n ( x ∗ w ) a0 + a1*cos(x*w) + b1*sin(x*w) a0+a1cos(xw)+b1sin(xw)
    • Gaussian:高斯逼近,有8种类型,基础型是 a 1 ∗ e x p ( − ( ( x − b 1 ) / c 1 ) 2 ) a1*exp(-((x-b1)/c1)^2) a1exp(((xb1)/c1)2)
    • Interpolant:插值逼近,有4种类型,Nearest neighbor、Linear、Cubic、Shape-preserving(PCHIP);
    • Linear Fitting:线性拟合;
    • Polynomial:多形式逼近;
    • Power:幂逼近,有2种类型, a ∗ x b a*x^b axb a ∗ x b + c a*x^b + c axb+c
    • Rational:有理数逼近;
    • Smoothing Spline:平滑逼近;
    • Sum of Sin Functions:正弦曲线逼近,有8种类型,基础型是 a 1 ∗ s i n ( b 1 ∗ x + c 1 ) a1*sin(b1*x + c1) a1sin(b1x+c1)
    • Weibull:只有一种, a ∗ b ∗ x ( b − 1 ) ∗ e x p ( − a ∗ x b ) a*b*x^(b-1)*exp(-a*x^b) abx(b1)exp(axb)

    6. 进行曲线拟合

    假设我们以 y = a ∗ s i n ( b ∗ x ) ∗ e x p ( c ∗ x ) + d y =a*sin(b*x)*exp(c*x)+d y=asin(bx)exp(cx)+d的表达形式进行拟合,则选择”Custom Equation”,在方框中输入相应的函数表达式,拟合过程及结果图像如下图所示:
    在这里插入图片描述
    可以发现,曲线和拟合程度较差。这是因为,对同一问题的拟合情况,每次可能都不一样,这取决对于参数a,b,c,d的StartPoint的选取。解决方法是在拟合过程中,将a,b,c和d也作为约束拟合的条件,例子中已知a,b,c,d的的取值范围(正负范围),可以进行如下操作,点击[Fit Options]按钮,弹出的窗口如下:
    在这里插入图片描述
    可以调整a,b,c,d参数的StartPoint,Lower,Upper三个选项来是拟合更加准确,比如说,将a,b,d的Lower选项设为0(a,b,d>0),将c的Upper选项设为0(c<0),设置如下图所示:
    在这里插入图片描述
    设置完毕之后,就会自动出现重新拟合之后的图像,如下图所示。
    在这里插入图片描述
    可以看到,拟合程度较之前有了很大的提高.因此,可以预见的是,在拟合过程中,设置好待拟合函数的参数的StartPoint,Lower和Upper三者的值可以使拟合更加准确。

    7. 拟合结果分析

    在左侧的Result中显示拟合模型参数以及拟合效果

    在这里插入图片描述
    拟合效果评测:

    • SSE:拟合误差平方和,接近0,表示与数据拟合的好,但是要小心过拟合;
    • R-Square:实测数据与推理数据之间的相关系数平方值,趋近于1较好;
    • RMSE:均方差;

    8.其他常用拟合方法

    当然,除了上面提到的拟合方法之外,还有两种常用的拟合方法:

    • Interpolant
      插值逼近,该方法的优势在于会连接所有点,而使其SSE为0,R-square为1,如下图所示:
      在这里插入图片描述
    • Smoothing Spline
      平滑逼近,该方法的会尽可能逼近所有点,使其SSE尽可能逼近0,R-square尽可能逼近1,如下图所示:
      在这里插入图片描述

    9.输出拟合参数

    如果希望只显示拟合图像,可以点击“文件”——>“Print to figure”
    在这里插入图片描述
    这样就可以只显示拟合图像了,如下图所示。
    在这里插入图片描述
    如果希望导出拟合后的曲线数据,可以点击“文件”——>“Generate Code”
    在这里插入图片描述
    生成代码后,默认函数名为createFit,可以自行修改,直接保存,就可以调用了。比如说,我要导出五次多项式Polynomial逼近结果,按照上述方式导出后,可以查看生成代码的信息:
    在这里插入图片描述
    其中,fitresult是函数的输出,是一个结构体,可以用fitresult.p1得到p1的系数,同理其它系数也可得。
    在这里插入图片描述
    如果想导出拟合后的曲线数据,只需要把横坐标传给fitresult就可以了:

    y = fitresult(x)';
    

    10.结论

    本文主要讨论了MATLAB曲线拟合工具箱(cftool)的拟合过程。通过工具箱模块可以非常方便地对曲线进行拟合,不需要太多的编程,曲线拟合方法多样,效果较好。

    ok,以上便是曲线拟合工具箱的全部内容了,如果对你有所帮助,记得点个赞哟~

    展开全文
  • 来喽!一、粒子群算法理论粒子群算法来源于鸟类集体活动的规律性,进而利用群体智能建立简化模型。它模拟的是鸟类的觅食行为...每个粒子有了初始位置与初始速度,然后迭代寻优。每一次迭代中,每个粒子通过跟踪两个...

    8f4917f943304f605dd2ed802c9649bf.png

    来喽!

    一、粒子群算法理论

    粒子群算法来源于鸟类集体活动的规律性,进而利用群体智能建立简化模型。它模拟的是鸟类的觅食行为,将求解问题的空间比作鸟类飞行的时间,每只鸟抽象成没有体积和质量的粒子,来表征一个问题的可行解。

    1.1 粒子群算法建模

    粒子群算法首先在给定的解空间中随机初始化粒子群,待优化问题的变量数决定了解空间的维数。每个粒子有了初始位置与初始速度,然后迭代寻优。每一次迭代中,每个粒子通过跟踪两个极值来更新自己的解空间中的位置和速度,一个是单个粒子本身在迭代中找到的最优粒子(个体极值),一个是所有粒子在迭代过程中的最优解粒子(全局极值)。

    1.2 粒子群算法特点

    (1)基于群体智能理论的优化算法,高效的并行算法。

    (2)粒子群算法随机初始化种群,使用适应值来评价个体的优劣程度,进行一系列的随机搜索。但粒子群算法根自己的速度来决定搜索,不需要进行交叉变异等操作,避免了复杂的操作。

    (3)每个粒子在算法结束时仍保持个体极值,所以除了得到问题的最优解以外,还能得到若干次优解,给出更多方案。

    (4)粒子群算法特有的记忆使其可以动态的跟踪当前搜索情况并调整搜索方向。

    二、粒子群算法种类

    2.1 基本粒子群算法

    假设在一个D维的目标搜索空间中,有N个粒子组成群落,其中第i个粒子表示为一个D维向量:

    第i个粒子飞行速度也是一个D维向量:

    第i个粒子迄今为止搜索到最优位置称为个体极值:

    整个粒子群迄今为止搜索到的最优位置为全局极值:

    在找到这两个最优值后,粒子群更新自己的速度和位置:

    其中,

    是学习因子(加速常数);
    为[0,1]范围内的均匀随机数;
    是粒子速度;
    增加了粒子飞行的随机性。更新速度的式子由三部分组成,第一部分为‘惯性’,反应了粒子群运动的‘习惯’;第二部分是‘认知’,反映了对自身历史经验的记忆;第三部分是‘社会’,反映了粒子间协同合作知识共享的群体历史经验。

    2.2 标准粒子算法

    引入了两个新的概念:

    探索:指粒子在一定程度上离开原先的轨迹,向新的方向搜索,体现了开拓未知区域的能力,探索能力是全局搜索能力。

    开发:指粒子在一定程度上继续在原先的轨迹上进行进一步搜索,是局部搜索能力。

    1998年,Shi Yuhui等人提出了带有惯性权重的改进粒子群算法:

    第一个式子有三个部分,第一部分保证算法的全局收敛性,第二、三部分保证局部收敛能力。w较大,全局收敛能力强局部收敛能力减弱;w较小,全局收敛能力减弱,局部收敛能力较强。当w=1时,与基本粒子算法相同,实验结果表明:w在0.8~1.2之间有更快的收敛速度,w>1.2时,容易陷入局部极值。

    在搜索时还可以对w进行动态调整,采用较多的是Shi提出的线性递减权值策略:

    表示最大进化代数;
    表示惯性最大权重;
    表示惯性最小权重。

    2.3 压缩粒子群算法

    Clerc等人利用约束因子来控制最终的收敛,可以有效搜索不同区域:

    为压缩因子:

    实验结果表明:与使用惯性权重的粒子群算法比较,使用具有压缩因子的粒子群算法有更快的收敛速度。

    2.4 离散粒子群算法

    Kennefy和Eberhart提出了离散二进制版的粒子群算法。将离散的问题映射到了连续的粒子运动空间,计算上扔保留速度-位置更新的运算规则。粒子在空间状态的取值只限于0,1两个值,而速度的每一维

    代表每一位
    取值为1的可能性,因此
    的公式保持不变,但是
    只能在[0,1]之间取值:

    686d14c77d38f36471bcf6a936e3f811.png

    r是从

    分布产生的随机数。

    三、粒子群算法流程

    4813d2d11f4e7b486c2bbdabd61d3e3a.png

    四、MATLAB实例及仿真

    参数说明:

    粒子种群规模N:一般为100~200。

    惯性权重w:控制算法的开发和探索能力,一般取值为[0.8,1.2]。

    加速常数

    :调节向
    方向飞行的最大步长,分别决定了粒子个体经验和群体经验对粒子运动轨迹的影响。
    粒子缺乏认知能力,
    个体之间没有信息交互,所以一般设置
    ,通常取
    。这样个体经验和群体经验有了同样的影响力。

    粒子最大速度

    :粒子速度在每一维都有一个最大速度限制值,用来对粒子速度进行钳制,该值太大,粒子也许会飞过优秀区域,太小粒子们可能会陷入局部最优,无法移动足够远跳出局部最优。研究发现 设定
    和调整惯性权重的作用是等效的,所以
    一般用于初始化进行设定,二不再对最大速度进行细致的选择和调节。

    邻域结构设定:全局版的粒子群算法将整个群体作为粒子的邻域,具有搜索快的优点但是容易陷入局部最优。局部版本粒子群算法将位置相近的个体作为粒子的邻域,收敛速度慢,不易陷入局部最优。实际中可先采用全局粒子群算法寻找最优解方向,再采用局部粒子群算法细致搜索。

    边界条件处理:边界吸收的方法。

    例1:计算函数

    的最小值,其中个体x的维数n=10,这是一个简单的平方和函数,只有一个极小点x=(0,0,...,0),理论上最小值分f(0,0,...,0)=0。

    解:

    clear all;
    close all;
    clc;
    N=100;                         %群体粒子个数
    D=10;                          %粒子维数
    T=200;                         %最大迭代次数
    c1=1.5;                        %学习因子1
    c2=1.5;                        %学习因子2
    w=0.8;                         %惯性权重
    Xmax=20;                       %位置最大值
    Xmin=-20;                      %位置最小值
    Vmax=10;                       %速度最大值
    Vmin=-10;                      %速度最小值
    %初始化个体
    x=rand(N,D)*(Xmax-Xmin)+Xmin;
    v=rand(N,D)*(Vmax-Vmin)+Vmin;
    %初始化个体最优位置和最优值
    p=x;
    pbest=ones(N,1);
    for i=1:N
        pbest(i)=func1(x(i,:));
    end
    %初始化全局最优位置和最优值
    g=ones(1,D);
    gbest=inf;
    for i=1:N
        if (pbest(i)<gbest)
            g=p(i,:);
            gbest=pbest(i);
        end
    end
    gb=ones(1,T);
    %按照公式依次迭代直到满足精度或者迭代次数
    for i=1:T
        for j=1:N
            if (func1(x(j,:))<pbest(j))
                p(j,:)=x(j,:);
                pbest(j)=func1(x(j,:));
            end
            if (pbest(j)<gbest)
                g=p(j,:);
                gbest=pbest(j);
            end
            v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))+c2*rand*(g-x(j,:));
            x(j,:)=x(j,:)+v(j,:);
            %边界条件处理
            for ii=1:D
                if (v(j,ii)<Vmin)||(v(j,ii)>Vmax)
                    v(j,ii)=rand*(Vmax-Vmin)+Vmin;
                end
                if (x(j,ii)<Xmin)|(x(j,ii)>Xmax)
                    x(j,ii)=rand*(Xmax-Xmin)+Xmin;
                end
            end
        end
        %记录全局最优值
        gb(i)=gbest;
    end
    g;                         %最优个体
    gb(end);                   %最优值
    figure
    plot(gb)
    xlabel('迭代次数')
    ylabel('适应度值')
    title('适应度进化曲线')
    %适应度函数
    function result=func1(x)
    summ=sum(x.^2);
    result=summ;
    end

    eabed64558542dcba992c8a20e7bec2f.png

    例2:求函数

    的最小值,其中x的取值范围是[-4,4],y的取值范围是[-4,4]。这是一个有多个局部极值的函数。函数图形如图:
    clear all;
    close all;
    clc;
    x=-4:0.02:4;
    y=-4:0.02:4;
    N=size(x,2);
    for i=1:N
        for j=1:N
            z(i,j)=3*cos(x(i)*y(j))+x(i)+y(j)*y(j);
        end
    end
    mesh(x,y,z)
    xlabel('x')
    ylabel('y')

    f5a7e45103cc84328cfbdb2e5d224088.png

    解:

    clear all;
    close all;
    clc;
    N=100;                         %群体粒子个数
    D=2;                           %粒子维数
    T=200;                         %最大迭代次数
    c1=1.5;                        %学习因子1
    c2=1.5;                        %学习因子2
    Wmax=0.8;                      %惯性权重最大值
    Wmin=0.4;                      %惯性权重最小值
    Xmax=4;                        %位置最大值
    Xmin=-4;                       %位置最小值
    Vmax=1;                        %速度最大值
    Vmin=-1;                       %速度最小值
    %初始化个体
    x=rand(N,D)*(Xmax-Xmin)+Xmin;
    v=rand(N,D)*(Vmax-Vmin)+Vmin;
    %初始化个体最优位置和最优值
    p=x;
    pbest=ones(1,D);
    for i=1:N
        pbest(i)=func2(x(i,:));
    end
    %初始化全局最优位置和最优值
    g=ones(1,D);
    gbest=inf;
    for i=1:N
        if (pbest(i)<gbest)
            g=p(i,:);
            gbest=pbest(i);
        end
    end
    gb=ones(1,T);
    %按照公式依次迭代直到满足精度或者迭代次数
    for i=1:T
        for j=1:N
            %更新个体最优位置和最优值
            if (func2(x(j,:))<pbest(j))
                p(j,:)=x(j,:);
                pbest(j)=func2(x(j,:));
            end
            %更新全局最优位置和最优值
            if (pbest(j)<gbest)
                g=p(j,:);
                gbest=pbest(j);
            end
            %计算动态惯性权重值
            w=Wmax-(Wmax-Wmin)*i/T;
            %更新位置和速度
            v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))+c2*rand*(g-x(j,:));
            x(j,:)=x(j,:)+v(j,:);
            %边界条件处理
            for ii=1:D
                if (v(j,ii)<Vmin)||(v(j,ii)>Vmax)
                   v(j,ii)=rand*(Vmax-Vmin)+Vmin;
                end
                if (x(j,ii)<Xmin)||(x(j,ii)>Xmax)
                    x(j,ii)=rand*(Xmax-Xmin)+Xmin;
                end
            end
        end
        gb(i)=gbest;
    end
    g;%最优个体
    gb(end);%最优值
    figure
    plot(gb)
    xlabel('迭代次数')
    ylabel('适应度值')
    title('适应度进化曲线')
    %适应度函数
    function value=func2(x)
    value=3*cos(x(1)*x(2))+x(1)+x(2)^2;
    end
            

    32150e88a0a8e765431b0aaa37c06bf6.png

    例3:用离散粒子群算法求函数

    的最小值,其中x的取值范围是[0,9],这是一个有局部多个极值的函数,图形如图所示:
    clear all;
    close all;
    clc;
    x=0:0.01:9;
    y=x+6*sin(4*x)+9*cos(5*x);
    plot(x,y)
    xlabel('x')
    ylabel('f(x)')

    e7b71574f3ece27f7be56e4a6eced536.png

    解:

    clear all;
    close all;
    clc;
    N=100;                         %群体粒子个数
    D=20;                           %粒子维数
    T=200;                         %最大迭代次数
    c1=1.5;                        %学习因子1
    c2=1.5;                        %学习因子2
    Wmax=0.8;                      %惯性权重最大值
    Wmin=0.4;                      %惯性权重最小值
    Xs=9;                        %位置最大值
    Xx=0;                       %位置最小值
    Vmax=10;                        %速度最大值
    Vmin=-10;                       %速度最小值
    %初始化个体
    x=randi([0,1],N,D);
    v=rand(N,D)*(Vmax-Vmin)+Vmin;
    %初始化个体最优位置和最优值
    p=x;
    pbest=ones(N,1);
    for i=1:N
        pbest(i)=func3(x(i,:),Xs,Xx);
    end
    %初始化全局最优位置和最优值
    g=ones(1,D);
    gbest=inf;
    for i=1:N
        if (pbest(i)<gbest)
            g=p(i,:);
            gbest=pbest(i);
        end
    end
    gb=ones(1,T);
    %按照公式依次迭代直到满足精度或者迭代次数
    for i=1:T
        for j=1:N
            %更新个体最优位置和最优值
            if (func3(x(j,:),Xs,Xx)<pbest(j))
                p(j,:)=x(j,:);
                pbest(j)=func3(x(j,:),Xs,Xx);
            end
               %更新全局最优位置和最优值
            if (pbest(j)<gbest)
                g=p(j,:);
                gbest=pbest(j);
            end
              %计算动态惯性权重值
              w=Wmax-(Wmax-Wmin)*i/T;
              %更新位置和速度
            v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))+c2*rand*(g-x(j,:));
           %边界条件处理
            for ii=1:D
                if (v(j,ii)<Vmin)||(v(j,ii)>Vmax)
                   v(j,ii)=rand*(Vmax-Vmin)+Vmin;
                end
            end
            vx(j,:)=1./(1+exp(-v(j,:)));
            for jj=1:D
               if vx(j,jj)>rand
                   x(j,jj)=1;
               else
                   x(j,jj)=0;
               end
            end
        end
        gb(i)=gbest;
    end
    g;%最优个体
    m=0;
    for j=1:D
        m=g(j)*2^(j-1)+m;
    end
    f1=Xx+m*(Xs-Xx)/(2^D-1);%最优值
    figure
    plot(gb)
    xlabel('迭代次数')
    ylabel('适应度值')
    title('适应度进化曲线')
    %适应度函数
    function result=func3(x,Xs,Xx)
    m=0;
    D=length(x);
    for j=1:D
        m=x(j)*2^(j-1)+m;
    end
    f=Xx+m*(Xs-Xx)/(2^D-1);%译码成十进制数
    fit=f+6*sin(4*f)+9*cos(5*f);
    result=fit;
    end
      

    3e4179b64fd905aa4ce474a9f5ba9d37.png

    例4:0-1背包问题,有N件物品和容积为V的包,第i件物品的容积是c(i),价值是w(i),求将这些物品放入包中,使物体的总容积不超过背包容积,且总价值和最大。假设物品数量为10,背包容量为300.每件物品的体积为[95,75,23,73,50,22,6,57,89,98];价值为[89,59,19,43,100,72,44,16,7,64];

    解:

    clear all;
    close all;
    clc;
    N=100;                          %群体粒子个数
    D=10;                           %粒子维数
    T=200;                          %最大迭代次数
    c1=1.5;                         %学习因子1
    c2=1.5;                         %学习因子2
    Wmax=0.8;                       %惯性权重最大值
    Wmin=0.4;                       %惯性权重最小值
    Vmax=10;                        %速度最大值
    Vmin=-10;                       %速度最小值
    V=300;                          %背包容量
    C=[95,75,23,73,50,22,6,57,89,98];            %物品体积
    W=[89,59,19,43,100,72,44,16,7,64];           %物品价值
    afa=2;                                       %惩罚函数系数
    %初始化个体
    x=randi([0,1],N,D);
    v=rand(N,D)*(Vmax-Vmin)+Vmin;
    %初始化个体最优位置和最优值
    p=x;
    pbest=ones(N,1);
    for i=1:N
        pbest(i)=func4(x(i,:),C,W,V,afa);
    end
    %初始化全局最优位置和最优值
    g=ones(1,D);
    gbest=eps;
    for i=1:N
        if (pbest(i)>gbest)
            g=p(i,:);
            gbest=pbest(i);
        end
    end
    gb=ones(1,T);
    %按照公式依次迭代直到满足精度或者迭代次数
    for i=1:T
        for j=1:N
             %更新个体最优位置和最优值
            if (func4(x(j,:),C,W,V,afa)>pbest(j))
                p(j,:)=x(j,:);
                pbest(j)=func4(x(j,:),C,W,V,afa);
            end
            %更新全局最优位置和最优值
            if (pbest(j)>gbest)
                g=p(j,:);
                gbest=pbest(j);
            end
             %计算动态惯性权重值
              w=Wmax-(Wmax-Wmin)*i/T;
              %更新位置和速度
              v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))+c2*rand*(g-x(j,:));
         %边界条件处理
            for ii=1:D
                 if (v(j,ii)<Vmin)||(v(j,ii)>Vmax)
                 v(j,ii)=rand*(Vmax-Vmin)+Vmin;
                 end
            end
            vx(j,:)=1./(1+exp(-v(j,:)));
            for jj=1:D
                 if vx(j,jj)>rand
                 x(j,jj)=1;
                 else
                 x(j,jj)=0;
                 end
            end
        end
        gb(i)=gbest;
    end
    g;%最优个体
    figure
    plot(gb)
    xlabel('迭代次数')
    ylabel('适应度值')
    title('适应度变化曲线')
    %适应度函数
    function result=func4(f,C,W,V,afa)
    fit=sum(f.*W);
    TotalSize=sum(f.*C);
    if TotalSize<=V
        fit=fit;
    else
        fit=fit-afa*(TotalSize-V);
    end
    result=fit;
    end
            

    f81bbfbb038657ea56255f7c16d77998.png

    希望我明天继续爱学习!!!

    展开全文
  • OFDM完整仿真过程及解释(MATLAB

    万次阅读 多人点赞 2019-04-19 17:03:45
    因为是复制过来,如果出现图片显示... OFDM完整仿真过程及解释(MATLAB) - 子木的文章 - 知乎 https://zhuanlan.zhihu.com/p/57967971 目录: 一、说明 二、ofdm总体概述 ...六、OFDM的MATLAB仿真程序 一、说...

    因为是复制过来,如果出现图片显示不完整以及需要源程序请点击下面链接查看原文:

    OFDM完整仿真过程及解释(MATLAB) - 子木的文章 - 知乎

    点击这里访问原文

    后面的更新没有同步,点上面链接可以看更新部分。

     

    目录:

    一、说明

    二、ofdm总体概述

    三、基本原理

    四、过程中涉及的技术

    五、OFDM基本参数的选择

    六、OFDM的MATLAB仿真程序

     

    一、说明

    0.能找到这篇文章,说明对ofdm已经有一点了解,所以其原理就不再赘述,这篇代码的目的只是希望能对ofdm整个过程有一个理解;

    1.看书上ofdm介绍挺简单的,自己来仿真才发现很多知识点都不知道;

    2.花了很长时间才理清整个ofdm过程,网上的程序都是一段一段的,不能直接理解整个过程。所以想着自己来做一个完整过程的仿真,加深理解;

    3.基带信号能完成整个过程,但是想加进频带传输这一部分,就完整了;

    4.信道部分想用瑞利信道的,程序写出来了,但是误差和信道估计这一块还不是很明白,所以就先用的高斯信道;

    5.不足之处欢迎指正。。。。

    二、概述

    OFDM是一种特殊的多载波传输方案,它可以被看作是一种调制技术,也可以被当作一种复用技术。

    简单来说:OFDM是一种多载波的传输方法,它将频带划分为多个子信道并行传输数据,将高速数据流分成多个并行的低速数据流,然后调制到每个信道的子载波上进行传输。由于它将非平坦衰落无线信道转化成多个正交平坦衰落的子信道,从而可消除信道波形间的干扰,达到对抗多径衰落的目的。

    正交频分复用(OFDM)是对多载波调制(MCM)的一种改进,在。它的特点是:各子载波相互正交,所以扩频调制后的频谱可以相互重叠,不但减少了子载波间的相互干扰,还大大提高了频谱利用率。

    选择OFDM的一个很大的原因是该系统能够很好的对抗频率选择性衰落和窄带干扰。在单载波系统中,一次衰落或者干扰会导致整个链路失效,但是在多载波系统中,某一时刻只会有少部分的子信道受到深衰落的影响。

    三、基本原理

    3.1 OFDM系统收发机的典型(根据实际需要可添/删部分)框图如下:

    OFDM收发机框图

    其中,上半部分对应于发射机链路,下半部分对应于接收机链路。

    发送端将被传输的数字信号转换成子载波幅度和相位的映射,并进行离散傅里叶变换(IDFT),将数据的频谱表达式变到时域上。IFFT和IDFT变换的作用相同,只是有更高的计算效率,所以适用于所有的应用系统。接收端进行与发送端相反的操作,用FFT变换分解,子载波的幅度和相位最终转换回数字信号。

    这里理解为传输的频域信号是因为IFFT是从频域到时域,实际上这里IFFT充当的是一个实现子载波正交的作用,具体可以推导其DFT公式。知乎里公式编辑太麻烦了。

     

    3.2 OFDM调制与解调

    一个OFDM符号之内包括多个经过调制的子载波的合成信号,其中每个子载波都可以收到psk(相移键控)和qam(正交幅度调制)的调制。

    OFDM发射机将信息比特流映射成一个psk或qam符号序列,之后将串行的符号序列转换为并行符号流。每N个经过串并转换的符号被不同的子载波调制。

    OFDM符号是N个并行符号的复合信号,若单个串行符号的传输时间(周期)是Ts,则一个OFDM符号的持续时间(周期)Tsym=N*Ts。

    频域调制信号X[k]的频率为:fk=k/Tsym,子载波数量为N,则k=0,1,2.....N-1。(由DFT原理推导)

    四、过程中涉及的技术

    为什么要用?怎么用?

    4.1 保护间隔

    多径信道会对OFDM符号造成ISI影响,破坏了子载波间的正交性。故需要采取一些方法来消除多径信道带来的符号间干扰(ISI)影响,即插入保护间隔。

    保护间隔有两种插入方法:一种是补零(zp),即在保护间隔中填充0;另一种是插入循环前缀(cp)或循环后缀(cs)实现OFDM的循环扩展(为了某种连续性)。

    zp是在保护间隔内不插入任何信号,但是在这种情况下,由于多径传播的影响,会产生载波间干扰(ICI),即不同的子载波间会产生干扰。

    一般采用cp。cp是将OFDM后部的采样复制到前面,长度为Tcp,故每个符号的长度为Tsym=Tsub+Tcp,Tsub为数据部分子载波数。Tcp大于或等于多径时延,符号间的ISI影响将被限制在保护间隔中,因此不会影响下一个OFDM的FFT变换。

    4.2交织

    交织的作用是将突发错误转换为随机错误,有利于前向纠错码的译码,提高了整个通信系统的可靠性。交织由两个变换过程组成。第一次变换保证了相邻的编码比特被映射到不相邻的子载波上。第二次变换保证了相邻的编码比特被分别映射到星座图的重要和非重要比特上,避免出现长时间的低比特位映射。

    交织块的长度Ncbps,对qpsk、16qam、64qam分别为2、4、6,s=Ncbps/2,d=16。

    4.3信道编码

    由于移动通信存在干扰和衰落,在信号传输过程中将出现差错,故对数字信号必须采用纠、检错技术,即纠、检错编码技术,以增强数据在信道中传输时抵御各种干扰的能力,提高系统的可靠性。

    这里的信道编码一般采用卷积编码,Viterbi译码。

    卷积编码是现代数字通信系统中常见的一种前向纠错码,区别于常规的线性分组码,卷积编码的码字输出不仅与当前时刻的信息符号输入有关,还与之前输入的信息符号有关。

    4.4 扩频

    “扩频通信技术是一种信息传输方式,其信号所占有的频带宽度远大于所传信息必需的最小带宽;频带的扩展是通过一个独立的码序列来完成,用编码及调制的方法来实现的,与所传信息数据无关;在接收端则用同样的码进行相关同步接收、解扩及恢复所传信息数据”

    根据香农定理,带宽和信噪比可用互换,扩频扩展了带宽,则对信噪比的要求可降低。

    4.5 导频

    导频不携带信息,导频是双方已知的数据,是用来做信道估计的。

    在接收机中,虽然利用接收到的段训练序列、长训练序列可以进行信道均衡、频率偏差校正,但符号还会存在一定的剩余偏差,且偏差会随着时间的累积而累积,会造成所有子载波产生一定的相位偏移。因此,还需要不断地对参考相位进行跟踪。要能实现这个功能,需要在52个非0子载波中插入导频符号。

    4.6 RF(射频)调制

    OFDM调制器的输出产生了一个基带信号,将此基带信号与所需传输的频率进行混频操作,利用模拟技术或数字上变频可完成。由于数字调制技术提高了处理I、Q信道之间的匹配性和数字IQ调制器相位的准确性,将会更加精确。

    五、OFDM基本参数的选择

    5.1 各种OFDM参数的选择就是需要在多项要求冲突中进行折衷考虑。通常来说,首先要确认3个参数:带宽、比特率、及保护间隔。

    5.1.1 按照惯例,保护间隔的时间长度应该为应用移动环境信道的时延扩展均方根值的2~4倍。

    5.1.2 确定保护间隔之后,则OFDM符号周期长度就确定了。为了最大限度的减少由于插入保护比特所带来的信噪比的损失,OFDM符号周期长度远远大于保护间隔长度。但是符号周期又不能任意大,否则就需要更多的子载波,带宽不变,子载波间隔就变小,系统的实现复杂度就提高了,而且还加大了系统的峰值平均功率比,同时系统对频率偏差更加敏感。因此,一般选择符号周期长度是保护间隔的5倍,这样,由于插入保护比特所造成的信噪比损耗只有1dB左右。

    5.1.3 确定保护间隔和符号周期长度之后,子载波的数量可由-3dB带宽除以子载波间隔(即去掉保护间隔之后的符号周期的倒数)得到。或者可由所要求比特速率除以每个子信道的比特速率来确定子载波的数量。每个信道中所传输的比特速率可由调制类型、编码速率、和符号速率来确定。

    5.2 有用符号持续时间T

    T对子载波之间间隔、译码的等待周期都有影响,为了保持数据的吞吐量,子载波数目和FFT的长度要有相对较大的数量,这就导致符号持续时间变长。总之,符号周期长度的选择以保证信道的稳定为前提。

    5.3 子载波数

    N=1/T

    其数值与FFT处理过的复数点数相对应,需适应数据速率和保护间隔的要求。

    5.4 调制模式

    OFDM系统的调制模式基于功率和频谱利用率来选择,可采用qam、psk。

    为了使所有的点有相同的平均功率,二进制序列映射后的复数要归一化。(BPSK\QPSK\16QAM\64QAM分别对应乘以1、1/根号2、1/根号10、1/根号42),解调的时候再变回去。

    5.5 以具体实例说明;

    要求:(1)比特率为25Mbit/s(2)可容忍的时延扩展为200ns(3)带宽小于18MHz。

    1)由200ns时延扩展得保护间隔为800ns;

    2)由保护间隔800ns得符号周期长度6*800ns=4.8us;

    3)子载波的间隔选取4.8-0.8=4us的倒数,即250KHz;

    4)由所要求的比特速率与OFDM符号速率的比值,每个符号需要传送的比特:25Mbit/s)/(1/4.8us)=120bit。

    5)为了完成上面120bit/符号,有两种选择:利用16QAM和码率为1/2的编码方法,这样每个子载波携带2bit的有用信息,因此需要60个子载波;另一种是利用QPSK和码率为3/4的编码方法,每个子载波携带1.5bit信息。因此需要80个子载波,然而80个子载波意外着带宽:80*250KHz=20MHz,大于所给带宽要求,故取第一种,即60个子载波。可利用64点IFFT来实现,剩余4个子载波补0.

    六、OFDM的MATLAB仿真主程序

    clc;
    clear;
    
    %————————————————————————————————————————————————————————%
    %q1:ifft点数难道不是应该等于子载波数吗?子载波数与ifft点数的关系?
    %a:ifft点数等于子载波数
    %q2:对矩阵进行fft?
    %a:y可以是一向量或矩阵,若y为向量,则Y是y的FFT,并且与y具有相同的长度。若y为一矩阵,则Y是对矩阵的每一列向量进行FFT。
    %q3:怎么对ofdm信号上变频
    %————————————————————————————————————————————————————————%
    
    %% 参数设置
    
    N_sc=52;      %系统子载波数(不包括直流载波)、number of subcarrier
    N_fft=64;            % FFT 长度
    N_cp=16;             % 循环前缀长度、Cyclic prefix
    N_symbo=N_fft+N_cp;        % 1个完整OFDM符号长度
    N_c=53;             % 包含直流载波的总的子载波数、number of carriers
    M=4;               %4PSK调制
    SNR=0:1:25;         %仿真信噪比
    N_frm=10;            % 每种信噪比下的仿真帧数、frame
    Nd=6;               % 每帧包含的OFDM符号数
    P_f_inter=6;      %导频间隔
    data_station=[];    %导频位置
    L=7;                %卷积码约束长度
    tblen=6*L;          %Viterbi译码器回溯深度
    stage = 3;          % m序列的阶数
    ptap1 = [1 3];      % m序列的寄存器连接方式
    regi1 = [1 1 1];    % m序列的寄存器初始值
    
    
    %% 基带数据数据产生
    P_data=randi([0 1],1,N_sc*Nd*N_frm);
    
    
    %% 信道编码(卷积码、或交织器)
    %卷积码:前向纠错非线性码
    %交织:使突发错误最大限度的分散化
    trellis = poly2trellis(7,[133 171]);       %(2,1,7)卷积编码
    code_data=convenc(P_data,trellis);
    
    
    %% qpsk调制
    data_temp1= reshape(code_data,log2(M),[])';             %以每组2比特进行分组,M=4
    data_temp2= bi2de(data_temp1);                             %二进制转化为十进制
    modu_data=pskmod(data_temp2,M,pi/M);              % 4PSK调制
    % figure(1);
    scatterplot(modu_data),grid;                  %星座图(也可以取实部用plot函数)
    
    %% 扩频
    %————————————————————————————————————————————————————————%
    %扩频通信信号所占有的频带宽度远大于所传信息必需的最小带宽
    %根据香农定理,扩频通信就是用宽带传输技术来换取信噪比上的好处,这就是扩频通信的基本思想和理论依据。
    %扩频就是将一系列正交的码字与基带调制信号内积
    %扩频后数字频率变成了原来的m倍。码片数量 = 2(符号数)* m(扩频系数)
    %————————————————————————————————————————————————————————%
    
    code = mseq(stage,ptap1,regi1,N_sc);     % 扩频码的生成
    code = code * 2 - 1;         %将1、0变换为1、-1
    modu_data=reshape(modu_data,N_sc,length(modu_data)/N_sc);
    spread_data = spread(modu_data,code);        % 扩频
    spread_data=reshape(spread_data,[],1);
    
    %% 插入导频
    P_f=3+3*1i;                       %Pilot frequency
    P_f_station=[1:P_f_inter:N_fft];%导频位置(导频位置很重要,why?)
    pilot_num=length(P_f_station);%导频数量
    
    for img=1:N_fft                        %数据位置
        if mod(img,P_f_inter)~=1          %mod(a,b)就是求的是a除以b的余数
            data_station=[data_station,img];
        end
    end
    data_row=length(data_station);
    data_col=ceil(length(spread_data)/data_row);
    
    pilot_seq=ones(pilot_num,data_col)*P_f;%将导频放入矩阵
    data=zeros(N_fft,data_col);%预设整个矩阵
    data(P_f_station(1:end),:)=pilot_seq;%对pilot_seq按行取
    
    if data_row*data_col>length(spread_data)
        data2=[spread_data;zeros(data_row*data_col-length(spread_data),1)];%将数据矩阵补齐,补0是虚载频~
    end;
    
    %% 串并转换
    data_seq=reshape(data2,data_row,data_col);
    data(data_station(1:end),:)=data_seq;%将导频与数据合并
    
    %% IFFT
    ifft_data=ifft(data); 
    
    %% 插入保护间隔、循环前缀
    Tx_cd=[ifft_data(N_fft-N_cp+1:end,:);ifft_data];%把ifft的末尾N_cp个数补充到最前面
    
    %% 并串转换
    Tx_data=reshape(Tx_cd,[],1);%由于传输需要
    
    %% 信道(通过多经瑞利信道、或信号经过AWGN信道)
     Ber=zeros(1,length(SNR));
     Ber2=zeros(1,length(SNR));
    for jj=1:length(SNR)
        rx_channel=awgn(Tx_data,SNR(jj),'measured');%添加高斯白噪声
        
    %% 串并转换
        Rx_data1=reshape(rx_channel,N_fft+N_cp,[]);
        
    %% 去掉保护间隔、循环前缀
        Rx_data2=Rx_data1(N_cp+1:end,:);
    
    %% FFT
        fft_data=fft(Rx_data2);
        
    %% 信道估计与插值(均衡)
        data3=fft_data(1:N_fft,:); 
        Rx_pilot=data3(P_f_station(1:end),:); %接收到的导频
        h=Rx_pilot./pilot_seq; 
        H=interp1( P_f_station(1:end)',h,data_station(1:end)','linear','extrap');%分段线性插值:插值点处函数值由连接其最邻近的两侧点的线性函数预测。对超出已知点集的插值点用指定插值方法计算函数值
    
    %% 信道校正
        data_aftereq=data3(data_station(1:end),:)./H;
    %% 并串转换
        data_aftereq=reshape(data_aftereq,[],1);
        data_aftereq=data_aftereq(1:length(spread_data));
        data_aftereq=reshape(data_aftereq,N_sc,length(data_aftereq)/N_sc);
        
    %% 解扩
        demspread_data = despread(data_aftereq,code);       % 数据解扩
        
    %% QPSK解调
        demodulation_data=pskdemod(demspread_data,M,pi/M);    
        De_data1 = reshape(demodulation_data,[],1);
        De_data2 = de2bi(De_data1);
        De_Bit = reshape(De_data2',1,[]);
    
    %% (解交织)
    %% 信道译码(维特比译码)
        trellis = poly2trellis(7,[133 171]);
        rx_c_de = vitdec(De_Bit,trellis,tblen,'trunc','hard');   %硬判决
    
    %% 计算误码率
        [err,Ber2(jj)] = biterr(De_Bit(1:length(code_data)),code_data);%译码前的误码率
        [err, Ber(jj)] = biterr(rx_c_de(1:length(P_data)),P_data);%译码后的误码率
    
    end
     figure(2);
     semilogy(SNR,Ber2,'b-s');
     hold on;
     semilogy(SNR,Ber,'r-o');
     hold on;
     legend('4PSK调制、卷积码译码前(有扩频)','4PSK调制、卷积码译码后(有扩频)');
     hold on;
     xlabel('SNR');
     ylabel('BER');
     title('AWGN信道下误比特率曲线');
    
     figure(3)
     subplot(2,1,1);
     x=0:1:30;
     stem(x,P_data(1:31));
     ylabel('amplitude');
     title('发送数据(以前30个数据为例)');
     legend('4PSK调制、卷积译码、有扩频');
    
     subplot(2,1,2);
     x=0:1:30;
     stem(x,rx_c_de(1:31));
     ylabel('amplitude');
     title('接收数据(以前30个数据为例)');
     legend('4PSK调制、卷积译码、有扩频');
    

    七、能看到这里,如果有丢丢帮助的话,emmmm点个赞~呗

    原文:

    整个过程

    本来对每一步都有讲解注释的,但是程序编辑多了感觉不美观,就删掉了。比如扩频,其原理、作用、如何实现~

    三、代码及说明

    1.尽量把每一句程序都注释,能达到初学者拿到程序就能懂的程度;

    2.下面这段程序是上变频之前的,包含了画图,对ofdm信号有一个直观的感受(与上面图片中的流程可能冲突,这里仅仅是为了画图解释,所以这也是最开始学容易绕晕的地方)

    clear;
    %% 参数设置
    sub_carriers=2048;%子载波数
    T = 1 / sub_carriers;
    time = [0:T:1-T];% Nifft份,每份相隔T
    
    Lp=4984;
    P_Tx=(rand(1,Lp)>0.5);%(bits)%产生1个长为Lp的数据包:
    conv_out=convolutional_en(P_Tx);%(卷积编码):
    interleave_table = interleav_matrix(ones(1,2*(Lp+8)));
    interleav_out = interleaving(conv_out ,interleave_table);%(交织器)
    
    x=qpsk(interleav_out);%(4QAM 调制)
    L=length(x);%信号长度
    
    s=48;
    symbol_used_len=L/s;%把输入分为S个符号,每个符号长为symbol_used_len
    %循环前缀的长度
    cp=256;
    %每一个OFDM符号的抽样值应补‘0’个数zeros_pad
    zeros_pad=sub_carriers-symbol_used_len;
    %每一个OFDM符号一侧应该补‘0’个数zeros_pad_side
    zeros_pad_side=zeros_pad/2;
    
    %对输入信号进行分割,分割为s个符号,再对每个符号进行FFT运算,实现OFDM解调,并保证能量不变
    time_domain_x_link=[];
    for I=0:(s-1)
        %对输入进行分割 
        x_temp=x(I*symbol_used_len+1:I*symbol_used_len+symbol_used_len);
        %对每个分割的部分进行补零操作,使其长为sub_carriers
        x_temp_pad=[zeros(1,zeros_pad_side),x_temp,zeros(1,zeros_pad_side)];
        %对每个符号进行IFFT运算
        time_domain_x_temp=ifft(x_temp_pad)*sqrt(sub_carriers);
        %对每个符号添加循环前缀
        time_domain_x_cp_temp=[time_domain_x_temp(sub_carriers-cp+1:sub_carriers),time_domain_x_temp];
        %将符号连接成为串行数据流
        time_domain_x_link=[time_domain_x_link,time_domain_x_cp_temp];
    
    end
    sum_xI = real(time_domain_x_link);
    sum_xQ = imag(time_domain_x_link);
    
    figure;
    num=1000;%画出前num个点  
    xaxis   = zeros(length(time(1:num)));
    plot(time(1:num), sum_xI(1:num), 'b:', time(1:num), sum_xQ(1:num), 'g:', time(1:num), abs(sum_xI(1:num)+j*sum_xQ(1:num)), 'k-', time(1:num), xaxis, 'r-');
    ylabel('y'),xlabel('t'),
    title(['前', num2str(num),'个点经ifft的QAM符号实部之和虚部之和以及实部与虚部的绝对值波形']),
    legend('实部之和','虚部之和', '绝对值');

    3.与上面图片流程相符的代码

    代码前面的问题也是我在这个过程中遇到的,困扰了好久,可以带着问题看看。欢迎讨论。

    clc;
    clear;
    
    %————————————————————————————————————————————————————————%
    %q1:fft点数难道不是应该等于子载波数吗?子载波数与ifft点数的关系?
    %q2:对矩阵进行fft?
    %q3:怎么对ofdm信号上变频
    %q4:基带速率是多少?怎么实现?
    %q5传输频带是多少?怎么实现?
    %q6子载波间隔是多少?怎么实现?
    %q7符号周期是多少?怎么实现?
    %————————————————————————————————————————————————————————%
    
    %% 参数设置
    
    N_sc=52;      %系统子载波数(不包括直流载波)、number of subcarrier
    N_fft=64;            % FFT 长度
    N_cp=16;             % 循环前缀长度、Cyclic prefix
    N_symbo=N_fft+N_cp;        % 1个完整OFDM符号长度
    N_c=53;             % 包含直流载波的总的子载波数、number of carriers
    M=4;               %4PSK调制
    SNR=0:1:25;         %仿真信噪比
    N_frm=10;            % 每种信噪比下的仿真帧数、frame
    Nd=6;               % 每帧包含的OFDM符号数
    P_f_inter=6;      %导频间隔
    data_station=[];    %导频位置
    L=7;                %卷积码约束长度
    tblen=6*L;          %Viterbi译码器回溯深度
    stage = 3;          % m序列的阶数
    ptap1 = [1 3];      % m序列的寄存器连接方式
    regi1 = [1 1 1];    % m序列的寄存器初始值
    
    
    %% 基带数据数据产生
    P_data=randi([0 1],1,N_sc*Nd*N_frm);
    
    
    %% 信道编码(卷积码、或交织器)
    %卷积码:前向纠错非线性码
    %交织:使突发错误最大限度的分散化
    trellis = poly2trellis(7,[133 171]);       %(2,1,7)卷积编码
    code_data=convenc(P_data,trellis);
    
    
    %% qpsk调制
    data_temp1= reshape(code_data,log2(M),[])';             %以每组2比特进行分组,M=4
    data_temp2= bi2de(data_temp1);                             %二进制转化为十进制
    modu_data=pskmod(data_temp2,M,pi/M);              % 4PSK调制
    % figure(1);
    scatterplot(modu_data),grid;                  %星座图(也可以取实部用plot函数)
    
    %% 扩频
    %————————————————————————————————————————————————————————%
    %扩频通信信号所占有的频带宽度远大于所传信息必需的最小带宽
    %根据香农定理,扩频通信就是用宽带传输技术来换取信噪比上的好处,这就是扩频通信的基本思想和理论依据。
    %扩频就是将一系列正交的码字与基带调制信号内积
    %扩频后数字频率变成了原来的m倍。码片数量 = 2(符号数)* m(扩频系数)
    %————————————————————————————————————————————————————————%
    
    code = mseq(stage,ptap1,regi1,N_sc);     % 扩频码的生成
    code = code * 2 - 1;         %将1、0变换为1、-1
    modu_data=reshape(modu_data,N_sc,length(modu_data)/N_sc);
    spread_data = spread(modu_data,code);        % 扩频
    spread_data=reshape(spread_data,[],1);
    
    %% 插入导频
    P_f=3+3*1i;                       %Pilot frequency
    P_f_station=[1:P_f_inter:N_fft];%导频位置(导频位置很重要,why?)
    pilot_num=length(P_f_station);%导频数量
    
    for img=1:N_fft                        %数据位置
        if mod(img,P_f_inter)~=1          %mod(a,b)就是求的是a除以b的余数
            data_station=[data_station,img];
        end
    end
    data_row=length(data_station);
    data_col=ceil(length(spread_data)/data_row);
    
    pilot_seq=ones(pilot_num,data_col)*P_f;%将导频放入矩阵
    data=zeros(N_fft,data_col);%预设整个矩阵
    data(P_f_station(1:end),:)=pilot_seq;%对pilot_seq按行取
    
    if data_row*data_col>length(spread_data)
        data2=[spread_data;zeros(data_row*data_col-length(spread_data),1)];%将数据矩阵补齐,补0是虚载频~
    end;
    
    %% 串并转换
    data_seq=reshape(data2,data_row,data_col);
    data(data_station(1:end),:)=data_seq;%将导频与数据合并
    
    %% IFFT
    ifft_data=ifft(data); 
    
    %% 插入保护间隔、循环前缀
    Tx_cd=[ifft_data(N_fft-N_cp+1:end,:);ifft_data];%把ifft的末尾N_cp个数补充到最前面
    
    %% 并串转换
    Tx_data=reshape(Tx_cd,[],1);%由于传输需要
    
    %% 信道(通过多经瑞利信道、或信号经过AWGN信道)
     Ber=zeros(1,length(SNR));
     Ber2=zeros(1,length(SNR));
    for jj=1:length(SNR)
        rx_channel=awgn(Tx_data,SNR(jj),'measured');%添加高斯白噪声
        
    %% 串并转换
        Rx_data1=reshape(rx_channel,N_fft+N_cp,[]);
        
    %% 去掉保护间隔、循环前缀
        Rx_data2=Rx_data1(N_cp+1:end,:);
    
    %% FFT
        fft_data=fft(Rx_data2);
        
    %% 信道估计与插值(均衡)
        data3=fft_data(1:N_fft,:); 
        Rx_pilot=data3(P_f_station(1:end),:); %接收到的导频
        h=Rx_pilot./pilot_seq; 
        H=interp1( P_f_station(1:end)',h,data_station(1:end)','linear','extrap');%分段线性插值:插值点处函数值由连接其最邻近的两侧点的线性函数预测。对超出已知点集的插值点用指定插值方法计算函数值
    
    %% 信道校正
        data_aftereq=data3(data_station(1:end),:)./H;
    %% 并串转换
        data_aftereq=reshape(data_aftereq,[],1);
        data_aftereq=data_aftereq(1:length(spread_data));
        data_aftereq=reshape(data_aftereq,N_sc,length(data_aftereq)/N_sc);
        
    %% 解扩
        demspread_data = despread(data_aftereq,code);       % 数据解扩
        
    %% QPSK解调
        demodulation_data=pskdemod(demspread_data,M,pi/M);    
        De_data1 = reshape(demodulation_data,[],1);
        De_data2 = de2bi(De_data1);
        De_Bit = reshape(De_data2',1,[]);
    
    %% (解交织)
    %% 信道译码(维特比译码)
        trellis = poly2trellis(7,[133 171]);
        rx_c_de = vitdec(De_Bit,trellis,tblen,'trunc','hard');   %硬判决
    
    %% 计算误码率
        [err,Ber2(jj)] = biterr(De_Bit(1:length(code_data)),code_data);%译码前的误码率
        [err, Ber(jj)] = biterr(rx_c_de(1:length(P_data)),P_data);%译码后的误码率
    
    end
     figure(2);
     semilogy(SNR,Ber2,'b-s');
     hold on;
     semilogy(SNR,Ber,'r-o');
     hold on;
     legend('4PSK调制、卷积码译码前(有扩频)','4PSK调制、卷积码译码后(有扩频)');
     hold on;
     xlabel('SNR');
     ylabel('BER');
     title('AWGN信道下误比特率曲线');
    
     figure(3)
     subplot(2,1,1);
     x=0:1:30;
     stem(x,P_data(1:31));
     ylabel('amplitude');
     title('发送数据(以前30个数据为例)');
     legend('4PSK调制、卷积译码、有扩频');
    
     subplot(2,1,2);
     x=0:1:30;
     stem(x,rx_c_de(1:31));
     ylabel('amplitude');
     title('接收数据(以前30个数据为例)');
     legend('4PSK调制、卷积译码、有扩频');

    4.上面就是整个基带传输过程,其实上变频和下变频也很简单,将信号分为IQ路,分别乘cos和-sin即可,关于这一步,可以参考我另一篇8PSK调制的文章,在那篇文章里有相似的原理。链接:https://blog.csdn.net/qq_41687938/article/details/89514982和贼详细的8PSK调制与解调详细过程 - 子木的文章 - 知乎 https://zhuanlan.zhihu.com/p/47258287

    5.本来打算解释解释原理的,但是想着网上资料很多,就不献丑了,想打王者了~

    展开全文
  • '适应度进化曲线' ) ; 完整代码下载: https://download.csdn.net/download/g425680992/10502951 粒子群算法的驱动因素 粒子群算法是一种随机搜索算法 。粒子的下一个位置受到自身历史经验和全局历史...

    从鸟群觅食行为到粒子群算法

    这里写图片描述

    鸟群寻找食物的过程中,鸟与鸟之间存在着信息的交换,每只鸟搜索目前离食物最近的鸟的周围区域是找到食物的最简单有效的办法。

    粒子群算法(以下简称PSO)就是模拟鸟群觅食行为的一种彷生算法 。 解=粒子=鸟 (鸟的位置象征着离食物的距离,粒子的位置也象征着离最优解的距离,是评价解质量的唯一标准), 找食物=找最优解,一个西瓜=一个粒子找到的历史最优解,一块肉=整个粒子群找到历史最优解 ,

    就像鸟的飞行路线会受到自己曾经寻找到的最优食物和鸟群曾经找到过的最优食物的双重影响一样,算法中,每一次迭代,粒子通过两个"极值"(全局历史最优解gBest和个体历史最优解pBest)来更新自己的速度,该速度又是更新粒子位置的关键,而粒子的位置象征着离最优解的距离,也是评价该粒子(解)的唯一标准 。

    粒子群算法的核心

    该算法的核心是如何根据pBest与gBest来更新粒子的速度和位置,标准粒子群给出了如下的更新公式:

    $ V_{t+1} =w \cdot V_t +c_1r_1\cdot(pBest-X_t) +c_2r_2\cdot(gBest-X_t) $

    X t + 1 = X t + V t + 1 X_{t+1} = X_t+V_{t+1} Xt+1=Xt+Vt+1

    $其中 , t:代数 , X是位置,V是速度,w是惯性权重,c是学习因子,r是随机数 $

    这里写图片描述

    如上图所示,假设这是一个在2维平面内寻找最优解的待求解问题,某一时间的某一粒子 X t X_t Xt处在原点位置 。则该粒子更新后的速度如上图所示 。 更新公式可以分为三个部门:

    • Part.1 : "惯性"或"动量"部分,反映粒子有维持自己先前速度的趋势
    • Part.2 : "认知"部门 , 反映粒子有向自身历史最优位置逼近的趋势
    • Part.3 : "社会"部门 , 反映粒子有向去群体历史最优位置逼近的趋势

    例 : 求解函数最小值

    ​ 求$f(x)=\sum_{i=1}{n}x_i2,(-20 \leq x\leq 20,n=10) $ 的最小值 ?

    % author zhaoyuqiang
    clear all ;
    close all ;
    clc ;
    N = 100 ; % 种群规模
    D = 10 ; % 粒子维度
    T = 100 ; % 迭代次数
    Xmax = 20 ;
    Xmin = -20 ;
    C1 = 1.5 ; % 学习因子1
    C2 = 1.5 ; % 学习因子2
    W = 0.8 ; % 惯性权重
    Vmax = 10 ; % 最大飞行速度
    Vmin = -10 ; % 最小飞行速度
    popx = rand(N,D)*(Xmax-Xmin)+Xmin ; % 初始化粒子群的位置(粒子位置是一个D维向量)
    popv = rand(N,D)*(Vmax-Vmin)+Vmin ; % 初始化粒子群的速度(粒子速度是一个D维度向量) 
    % 初始化每个历史最优粒子
    pBest = popx ; 
    pBestValue = func_fitness(pBest) ; 
    %初始化全局历史最优粒子
    [gBestValue,index] = max(func_fitness(popx)) ;
    gBest = popx(index,:) ;
    for t=1:T
        for i=1:N
            % 更新个体的位置和速度
            popv(i,:) = W*popv(i,:)+C1*rand*(pBest(i,:)-popx(i,:))+C2*rand*(gBest-popx(i,:)) ;
            popx(i,:) = popx(i,:)+popv(i,:) ;
            % 边界处理,超过定义域范围就取该范围极值
            index = find(popv(i,:)>Vmax | popv(i,:)<Vmin);
            popv(i,index) = rand*(Vmax-Vmin)+Vmin ; %#ok<*FNDSB>
            index = find(popx(i,:)>Xmax | popx(i,:)<Xmin);
            popx(i,index) = rand*(Xmax-Xmin)+Xmin ;
            % 更新粒子历史最优
            if func_fitness(popx(i,:))>pBestValue(i)    
               pBest(i,:) = popx(i,:) ;
               pBestValue(i) = func_fitness(popx(i,:));
            end
           if pBestValue(i) > gBestValue
                gBest = pBest(i,:) ;
                gBestValue = pBestValue(i) ;
           end
        end
        % 每代最优解对应的目标函数值
        tBest(t) = func_objValue(gBest); %#ok<*SAGROW>
    end
    figure
    plot(tBest);
    xlabel('迭代次数') ;
    ylabel('适应度值') ;
    title('适应度进化曲线') ;
    

    完整代码下载:https://download.csdn.net/download/g425680992/10502951

    这里写图片描述

    粒子群算法的驱动因素

    粒子群算法是一种随机搜索算法 。粒子的下一个位置受到自身历史经验和全局历史经验的双重影响,全局历史经验时刻左右着粒子的更新,群体中一旦出现新的全局最优,则后面的粒子立马应用这个新的全局最优来更新自己,大大提高了效率,相比与一般的算法(如遗传算法的交叉),这个更新过程具有了潜在的指导,而并非盲目的随机 。

    自身历史经验和全局历史经验的比例尤其重要,这能左右粒子的下一个位置的大体方向,所以,粒子群算法的改进也多种多样,尤其是针对参数和混合其他算法的改进 。

    总体来说,粒子群算法是一种较大概率收敛于全局最优解的,适合在动态、多目标优化环境中寻优的一种高效率的群体智能算法。

    展开全文
  • 初始值K=31.5 % The Bode plot of the liquid level control system with % gain and phase margin. The gain K must be input at the % command level before executing this m-file. Automatic % labeling of the...
  • 在直接数字域设计中,我们常常需要用到PID算法,而PID算法投入单片机使用时,往往需要硬件的支持,在调试时非常麻烦。本文通过Matlab仿真的手段实现PID,方便了开发者对系统的设计和实时调试。
  • matlab人脸识别论文

    万次阅读 多人点赞 2019-10-11 17:41:51
    在开始训练前,所有的权值和阈值都应该用一些不同的小随机数进行初始化。BP算法主要包括两个阶段: 2.2.1向前传播阶段 ①从样本集中取一个样本(Xp,Yp),将Xp输入网络,其中Xp为输入向量,Yp为期望输出向量。 ②...
  • 实验环境:matlab2016a数据集:数据集大小3*3000,表示3000个样本,每个样本包含2个特征,第三行表示样本所属的分类。对于此次实验编写的BP神经网络和RBF神经网络,均将原始数据集分为训练集和测试集两部分,训练集...
  • MATLAB】关于ode45的一部分用法

    万次阅读 多人点赞 2018-06-24 16:22:36
    当难以求得微分方程的解析解时,可以求其数值解(解析解就是给出解的具体函数形式,从解的表达式中就可以算出任何对应;数值解就是用数值方法求出近似解,给出一系列对应的自变量和解)。  Matlab中求微分方程...
  • 具有初始值initial()状态空间模型 5 根表roots() 23 根表solve() 24 根轨迹rlocus() 27 相位裕度,增益裕margin() 52 PID控制,余量 27 PID控制,步进参考 27 PID控制,灵敏度 28岁 Lesson05_4.slx PID控制,...
  • ROC曲线MATLAB实现以及AUC

    千次阅读 2017-11-26 10:00:50
    ROC曲线(Receiver Operating Characteristic Curve)是利用Classification模型真正率(True Positive Rate)和假正率(False Positive Rate)作为坐标轴,图形化表示分类方法的准确率的高低。 ROC(Receiver ...
  • Matlab一元非线性回归分析

    万次阅读 2018-12-27 21:11:11
    (2)根据散点图的走势,确定回归方程的具体形式,特别是参数个数的设定和设定初始值; (3)调用NonLinearModel的fit方法进行模型拟合; (4)模型改进,去除异常值的操作; (5)进行残差分析,验证模型。 ...
  • MATLAB实现中频正交采样(数字下变频)

    万次阅读 多人点赞 2019-01-09 22:08:12
    % 调用MATLAB butter函数快速设计滤波器 [BH,BW]=freqz(Bb,Ba); % 绘制频率响应曲线 Bf=filter(Bb,Ba,y); % 进行低通滤波 By=fft(Bf,n); % 对滤波后信号做n点FFT变换 figure(3); plot(t,y,'blue',t,Bf,'red')...
  • 遗传算法的matlab实现

    万次阅读 多人点赞 2017-06-23 18:14:07
    遗传算法是一种全局最优化算法,是运用了进化论...仿真平台使用的是matlab,主要使用的是谢菲尔德大学的matlab遗传算法工具箱。 具体代码如下: clc clear all close all lbx=-1;ubx=1; %函数自变量x范围【-1,1】 lby=-
  • 神经网络中权重的初始值如何设定?

    千次阅读 2019-07-16 11:25:34
    在神经网络的学习中,权重的初始值特别重要。实际上,设定什么样的权重初始值,经常关系到神经网络的学习能否成功。本文将介绍在Sigmod、tanh和Relu激活函数下权重初始值的推荐值。
  • “init”是对参数 [abcde] 的初始猜测。 如果为空,它们将根据输入数据自动确定,“w”是权重向量(默认为ones(size(x)))。 输出:“f”是拟合函数的。 “X”是估计参数。 “err”是标准化错误。 “它”是迭代...
  • [来源维基百科] % % 输入: % 点数:点数列表 2xN % epsilon:距离维度,指定之间的相似度% 原始曲线和近似(较小的 epsilon, % 曲线更相似) % 输出: % 结果:近似曲线 2xM (M<=N) 的点列表%
  • 差分进化算法之Matlab实现

    千次阅读 多人点赞 2019-03-25 18:04:24
    适应度函数变化曲线为: matlab代码为: % clear all; % close all; % size=50;%群体个数 Codel=2;%所求的变量个数 MinX(1)=-5;%未知量范围 MinX(2)=-5; MaxX(1)=5; MaxX(2)=5; G=200;%迭代次数 F=1.2;%...
  • 解决此类问题有以下几个步骤 1首先作出散点图确定函数的类别 2根据已知数据确定待定参数的初始值利用 Matlab软件计算最佳参数 3根据可决系数比较拟合效果计算可决系数的公式为;一数据预处理;1.输入和查看数据集;5个...
  • Matlab曲线拟合

    千次阅读 2014-05-06 23:49:05
    Excel拟合曲线方程    Excel可以通过画散点图,添加趋势线,对数据进行简单的对数,线性,多项式,指数,幂函数曲线拟合,可能不适用于它们的组合,也就是交复杂的线性;另一种方法是使用加载项中的规划求解,...
  • Matlab对指定参数的曲线进行非线性拟合

    万次阅读 多人点赞 2019-01-02 20:12:28
    Matlab拟合曲线的方式 Matlab拟合曲线的方式有很多种,有三次样条插值、线性插值、多项式拟合等等。多项式拟合由于函数由f(x)=anxn+an−1xn−1+...+a1x+a0f(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_1x+a_0f(x)=an​xn+an−1...
  • 要求二:用RMSE实现两种算法的性能比较, 得到两种算法的RMSE曲线对比图,横坐标为噪声方差,纵坐标为RMSE。 1.2、toa定位代码 TOA: % %% the simulation of TOA localization algorithm clear all; clc; %...
  • matlab曲线的颜色代码FisherLuDAlessandroWilson_AnalysisCode Fisher等人的分析代码。 2019年手稿。 图1: analyticsOpenLoopTuningScript.m 视觉调整曲线的初始处理 plotOrderedReceptiveFields.m 绘制E-PG视觉...
  • 双方博弈复制动态方程(合并了初始值和参数优化后的曲线,两者相比较,进行参数敏感性分析)与matlab数值仿真——matlab2016a版本 1.主函数 %% 一、初始值+将Δa_A、Δa_B分别增加至10、11 %①初始值 a_A1=14;a_B1=...
  • MATLAB app designer设计人机交互界面——二阶线性动态电路可视化分析的研究 这是我第一次尝试写博客,我试着给出电路课上要求的电路实验编程。但是电路的类型有点儿多,所以我只以二阶动态电路RCL,进行全响应...
  • 下面将讲解Koch曲线的形成原理及matlab程序实现(对空间上的任意两个点,画出它们的Koch曲线) 形成过程 对于二维空间上的任意一条线段,将该线段的两个端点称为首点和尾点,我们需要利用这首点和尾点来得到新的五个...
  • 数组说明中初始值表可以省略,这时其初始值为相应数组元素的值(如果其数组元素还没有值则初值为缺失值)。 数组说明中的数组元素名列表可以省略,这时其元素也有对应的变量名,变量名为数组名后附加序号,比如: ...
  • Matlab代码实践——BP神经网络

    千次阅读 2019-07-16 12:03:11
    net.trainParam.mu 的初始值(缺省为0.001)  trainlm  net.trainParam.mu_dec 的减小率(缺省为0.1)  trainlm  net.trainParam.mu_inc 的增长率(缺省为10)  trainlm  net.trainParam.mu_max ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 6,021
精华内容 2,408
关键字:

matlab曲线初始值

matlab 订阅