gai_juanji.m

这是关于msk的调制解调。用卷积的方法进行调制 程序的最后还给出了误码率的计算

clear;clc;

Ts=0.25; fb=4;         % fb=1/ts; ts 码元宽度
fc=8;                 % 振荡频率
N=50;                 % 产生随机序列的点数
ntb=100;              % 每个码元采样点数
tmin=1/ntb;           % 最小时间间隔
sizeofnoise=5;       % 加入高斯噪声大小(dB) 

%产生伪随机序列,并将其转换为+1,-1的码元序列
temp=rand(1,N);
for s=1:N
    if temp(s)>=0.5
        a(s)=1;
    else
        a(s)=0;
    end
end

% 卷积码编码
trel=poly2trellis(3,[4 5 7]);
code=convenc(a,trel);

% -1和0的转换
len=length(code);
for i=1:len
    if code(i)==1
        a1(i)=1;
    else 
        a1(i)=-1;
    end
end

%差分编码
p(1)=a1(1);
for i=2:len
    if a1(i)==p(i-1)
        p(i)=1;
    else p(i)=-1;
    end
end

%将差分编码后的p序列进行串并转换分为两路bi和bq
l=length(a1)/2;
for m=1:l
    bq(m)=p(2*m);
    bi(m)=p(2*m-1);
end 
% MSK要求产生的Q路信号与I路信号相比延迟一个码元周期 ts
% 产生I路信号mbi  将mbi提前一个“码元持续时间Tb”

for j=1:l
            k=round((j-1+tmin)*ntb);
            n=round(j*ntb);
            Ti=(j-1.5+tmin)*Ts/0.5:tmin*Ts/0.5:(j-0.5)*Ts/0.5;
            mbi(k:n)=bi(j)*cos(pi/2*fb.*Ti).*cos(2*pi*fc.*Ti);
end

%产生Q路信号mbq
for j=1:l
            k=round((j-0.5+tmin)*ntb);
            n=round((j+0.5)*ntb);
            Tq=(j-1+tmin)*Ts/0.5:tmin*Ts/0.5:j*Ts/0.5;
            if n>l*ntb  %将最后一个Tb内的mbq信号循环回来
                mbq(1:ntb/2)=bq(j)*sin(pi/2*fb.*Tq(ntb/2+1:ntb)).*sin(2*pi*fc.*Tq(ntb/2+1:ntb));
                mbq(k:k+ntb/2-1)=bq(j)*sin(pi/2*fb.*Tq(1:ntb/2)).*sin(2*pi*fc.*Tq(1:ntb/2));
            else
            mbq(k:n)=bq(j)*sin(pi/2*fb.*Tq).*sin(2*pi*fc.*Tq);
            end
end

%两路信号相加得到msk信号
for i=1:100*l
    msk(i)=mbi(i)-mbq(i);
end

subplot(3,1,1),plot(msk(1:300))
grid on
subplot(3,1,2),plot(msk(301:600))
grid on
% 
 n1=0:4095;
 y1=fft(msk,4096);
 df=2*ntb/4096*Ts;
 w1=df*n1;
subplot(3,1,3),plot(w1,abs(y1))
figure;
subplot(2,1,1),plot(mbi(1:300))
grid on
subplot(2,1,2),plot(mbq(1:300))
grid on
% 加入高斯白噪声
noise_msk=awgn(msk,-5,'measured');
figure;subplot(2,1,1);plot(noise_msk(1:300));grid on
       subplot(2,1,2);plot(noise_msk(301:600));grid on
% MSK信号的解调
t=tmin*Ts/0.5:tmin*Ts/0.5:Ts*len;
% 调制信号和载波x(t),y(t)相乘,得到I路,Q路信号
receive_i=msk.*cos(pi/(2*Ts).*t).*cos(2*pi*fc.*t);
receive_q=msk.*sin(pi/(2*Ts).*t).*sin(2*pi*fc.*t);
% 画出和本地载波相乘后的I路和Q路的波形图
figure;subplot(2,1,1);plot(t,receive_i)
       title('receive_i waveform');grid on
       subplot(2,1,2);plot(t,receive_q)
       title('receive_q waveform');grid on
       
% 设计巴特沃兹低通滤波器，并画出其频率响应曲线
[b,aa]=butter(10,3/100);
freqz(b,aa);

% 得到的I路和Q路信号通过低通滤波器
% 滤除 2wc 高频成分，得到I路和Q路的原始信号
receive_i_lowpass=filter(b,aa,receive_i);
receive_q_lowpass=filter(b,aa,receive_q);

% % 采用积分的方法提取I路和Q路信息
% 
% for i=1:len
%     s1((i-1)*100+1)=receive_i((i-1)*100+1);
%     for j=2:100
%         s1((i-1)*100+j)=s1((i-1)*100+j-1)+receive_i((i-1)*100+j);
%     end
% end
% 
% for i=1:len
%     s2((i-1)*100+1)=receive_i((i-1)*100+1);
%     for j=2:100
%         s2((i-1)*100+j)=s2((i-1)*100+j-1)+receive_q((i-1)*100+j);
%     end
% end
% 画出滤波后I路和Q路的信号，并与输入的I路和Q路信号做比较
ii=Ts/0.5:Ts/0.5:l*Ts/0.5;
figure;subplot(2,2,1),plot(t,receive_i_lowpass)
       hold on
       stem(ii,bi)
       title('receive_i_lowpass waveform');grid on
       subplot(2,2,2),plot(t,receive_q_lowpass)
       hold on
       stem(ii,bq)
       title('receive_q_lowpass waveform');grid on
       subplot(2,2,3),stem(ii,bi);grid on
       subplot(2,2,4),stem(ii,bq);grid on
       
% 抽样判决，I路信号
for i=1:l
    if receive_i_lowpass((i-1)*ntb+ntb)>0
        out_i(i)=1;
    else
        out_i(i)=-1;
    end 
end
% 抽样判决，Q路信号
% 因为I路信号和Q路信号有一个码元周期的延时，导致Q路最后一个码元
% 在指定的长度内不能被恢复出来，程序中对其赋值-1处理。
for i=1:l
    if (i-1)*ntb+2ntb>l*ntb
        out_q(i)=-1;
    else
    if receive_q_lowpass((i-1)*ntb+2*ntb)>0
        out_q(i)=1;
    else
        out_q(i)=-1;
    end 
    end
end
% 完成I路信号和Q路信号的并串转换
 for i=0:l-1
    decis_data(i*2+1)=out_i(i+1);
    decis_data(i*2+2)=out_q(i+1);
 end
% 差分译码
 for i=2:len
     out_data(1)=a1(1);
     if  decis_data(i)==1
         out_data(i)=decis_data(i-1);
     else
         out_data(i)=-decis_data(i-1);
     end
 end

% -1和数字0之间的转换
for i=1:len
    if out_data(i)==-1
        out_data1(i)=0;
    else
        out_data1(i)=out_data(i);
    end
end

% 卷积码译码
decode=vitdec(out_data1,trel,8,'cont','hard');

decode
% 计算误码率
count=0;
for i=1:len-1
    if a1(i)==out_data(i)
        count=count+1;
    else
        i
    end
end
pe_msk=(len-1-count)/(len-1);
pe_msk

count1=0;
for i=1:N-8
    if a(i)==decode(i+8)
        count1=count1+1;
    end
end
pe_code=(N-8-count1)/(N-8);
pe_code