numbits = 99;   %每次迭代的比特數
numchans = 2;   %信號的信道數
nsamp = 16;    %每個符號的采樣率
num=10;
numerrs = 0; % Number of bit errors seen so far
demod_ini_phase = zeros(1,numchans);  % 調制相位
mod_ini_phase = zeros(1,numchans);    % 解調相位
ini_state = complex(zeros(nsamp,numchans)); % 解調狀態

for i = 1 : num
    x = randint(numbits,numchans);   %二進制信號
    [y,phaseout] = mskmod(x,nsamp,[],mod_ini_phase);
    mod_ini_phase = phaseout;    %用于下次MSK調制使用
    [z, phaseout, stateout] =mskdemod(awgn(y,5,'measured'),nsamp,[],demod_ini_phase,ini_state);
    ini_state = stateout;  %用于下次解調使用
    demod_ini_phase = phaseout;   %用于下次解調使用
%統計誤比特率
numerrs = numerrs + length(find(x(:,1)~=z(:,1))) + length(find(x(:,2)~=z(:,2)));  
end
ber=numerrs/(numbits*numchans*num);