MATLAB连续系统离散仿真

img

要编写MATLAB程序，有t又有s，要怎么写，还有差分方程怎么写进程序里？

不知道你这个问题是否已经解决, 如果还没有解决的话:
  • 这有个类似的问题, 你可以参考下: https://ask.csdn.net/questions/376952
  • 这篇博客也不错, 你可以看下matlab索引超过矩阵错误，这个程序哪里有问题呀~
  • 除此之外, 这篇博客: Matlab实现 维特比译码中的 话不多说，直接展示代码 部分也许能够解决你的问题, 你可以仔细阅读以下内容或者直接跳转源博客中阅读:
    1. 首先是利用卷积码进行编码
    2. 然后加入awgn噪声
    3. 最后利用维特比译码
    4. 绘制误比特率图像

    同时我也上传所有的代码文件

    Matlab 实现维特比译码

    在使用时，首先运行encode，对于随机产生的N=1000个比特流进行卷积编码
    然后直接运行show errorbit，它会调用维特比译码的函数viterbi_hard()，同时绘制出在加入awgn白噪声后，不同信噪比下的误比特率。

    卷积码进行编码
    生成向量为 g1: [1,1,1], g2: [1,0,1]

    encode

    %%  初始化参数
N = 1000; %序列长度

code_in = randi(2,1,N)-1;
%code_in = [1,1,0,1,1,1,0,0,1];
N =  length(code_in);
g1 = [1,1,1];
g2 =  [1,0,1];

%% 计算 卷积码
%x_g1

x0 = g1(1) * code_in;
x1 = g1(2) * code_in;
x2 = g1(3)* code_in;
x = xor( x0(2:N), x1(1:N-1) );
x = [ x0(1), x, x1(N) ];
x3 = xor(x(3:N+1),x2(1:N-1));
x_g1 = [x(1:2) ,x3, x2(N) ];
%x_g2
x0 = g2(1) * code_in;
x1 = g2(2) * code_in;
x2 = g2(3)* code_in;
x = xor( x0(2:N), x1(1:N-1) );
x = [ x0(1), x, x1(N) ];
x3 = xor(x(3:N+1),x2(1:N-1));
x_g2 = [x(1:2) ,x3, x2(N) ];

%% 合并两路

 x = zeros(1,size(x_g1, 2)+size(x_g2, 2));
 x(1:2:end) = x_g1;
 x(2:2:end) = x_g2;
x = x(1:length(x) - 4);

    维特比译码
    将状态转移图存在graph中
    decode

    function m = viterbi_hard(x)   
%% 
% a 00
% b 10
% c 01
% d 11
%x = [ 1 1 0 1 0 1 0 0 0 1];
sizex = size(x);
s = sizex(2)/2;
x = [0,0,0,0, x];
% to record the value 
val_a = 0;
val_b = 0;
val_c = 0;
val_d = 0;

% graph       aa0    ab1    bc0   bd1    ca0     cb1   dc0   dd1      
gra =          [ 0,0;    1,1;    1,0;    0,1;    1,1;    0,0;    0,1;    1,0  ];

% to record route 
ma = zeros(1,s+2);
mb = zeros(1,s+2);
mc = zeros(1,s+2);
md = zeros(1,s+2);
%% stage 1 and 2

val_a = val_a + dis(gra(1,:), x(1:2));
ma(1)=0;
val_b = val_b + dis(gra(2,:), x(1:2));
mb(1)=1;

mc = mb;
md = mb;
val_c = val_b + dis(gra(3,:), x(3:4));
mc(2)=0;
val_d = val_b + dis(gra(4,:), x(3:4));
md(2)=1;

val_a = val_a + dis(gra(1,:), x(3:4));
ma(2)=0;
val_b = val_b + dis(gra(2,:), x(3:4));
mb(2)=1;

for i = 3:s+2
%     val_a_t =val_a;
%     val_b_t =val_b;
%     val_c_t =val_c;
%     val_d_t =val_d;
%     tempa = ma;
%     tempb = mb;
%     tempc = mc;
%     tempd = md;
    % for val_a
        if val_a + dis(gra(1,:), x(2*i-1:2*i)) >= val_c + dis(gra(5,:),x(2*i-1:2*i))
            tempa = mc; 
            val_a_t = val_c + dis(gra(5,:),x(2*i-1:2*i));
            tempa(i)=0;
        else
            val_a_t = val_a + dis(gra(1,:),x(2*i-1:2*i));
            tempa = ma;
            tempa(i)=0;
        end
      %for val_b
         if val_a + dis(gra(2,:), x(2*i-1:2*i)) >= val_c + dis(gra(6,:),x(2*i-1:2*i))
            tempb = mc; 
            val_b_t = val_c + dis(gra(6,:),x(2*i-1:2*i));
            tempb(i)=1;
        else
            val_b_t = val_a + dis(gra(2,:),x(2*i-1:2*i));
            tempb = ma;
            tempb(i)=1;         
         end
        
         %for val_c 
            if val_b + dis(gra(3,:), x(2*i-1:2*i)) >= val_d + dis(gra(7,:),x(2*i-1:2*i))
            tempc = md; 
            val_c_t = val_d + dis(gra(7,:),x(2*i-1:2*i));
            tempc(i)=0;
            else
            val_c_t = val_b + dis(gra(3,:),x(2*i-1:2*i));
            tempc = mb;
            tempc(i)=0;
            end
            
      %for val_c
            if val_b + dis(gra(4,:), x(2*i-1:2*i)) >= val_d + dis(gra(8,:),x(2*i-1:2*i))
            tempd = md; 
            val_d_t = val_d + dis(gra(8,:),x(2*i-1:2*i));
            tempd(i)=1;
        else
            val_d_t = val_b + dis(gra(4,:),x(2*i-1:2*i));
            tempd = mb;
            tempd(i)=1;
            end
    val_a =val_a_t;        
    val_b =val_b_t;
    val_c =val_c_t;
    val_d =val_d_t;
    ma = tempa;
    mb = tempb;
    mc = tempc;
    md = tempd;        
end

if val_a <= val_b
    m = ma;
    t = val_a;
else
    m = mb;
   t = val_a;
end

if val_c <= t
    m = mc;
    t =val_c;
end

if val_d <= t
    m = md;
    t = val_d;
end
%% 去掉最开始的00
m = m(3:s+2);

    译码中间用到dis 函数，测量两段序列的海明距离

    dis

    %% function to compute distances between x,y in the form of [0,0]  [1,1]
function d = dis(x, y)   
  d = sum(xor(x,y));

    加入噪声，绘制误比特率随snr变化的曲线
    show error bit

    
errbit  = zeros(1,21);
for j = -5:15 

y = awgn(x, j, 'measured');
for i = 1:2*N
    if y(i) >=0.5
        y(i)=1;
    else
        y(i)=0;
    end
end


% 译码
m=viterbi_hard(y); 
errbit(j+6) = sum(m~=code_in)/N ;

end

logerr = 10*log10(errbit);
plot( -5:15, logerr );

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^