gai.m

这是MSK调制与解调的新方法

clear;clc;

Ts=0.125; fb=8;         % fb=1/ts; ts 码元宽度
fc=68;                 % 振荡频率
N=50;                 % 产生随机序列的点数
ntb=1000;              % 每个码元采样点数
tmin=1/ntb;           % 最小时间间隔
sizeofnoise=10;       % 加入高斯噪声大小(dB) 
k_noise=0.2;          % 加入工频干扰的幅度
fg=50;                % 加入工频干扰的频率


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

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

%将差分编码后的p序列进行串并转换分为两路bi和bq
l=length(a)/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:ntb*l
    msk(i)=mbi(i)-mbq(i);
end

subplot(3,1,1),plot(msk(1:3000))
grid on
subplot(3,1,2),plot(msk(3001:6000))
grid on
% 
 n1=0:4095;
 y1=fft(msk,4096);
 df=ntb*0.5/(Ts*4096);
 w1=df*n1;
subplot(3,1,3),plot(w1,abs(y1))
figure;
subplot(2,1,1),plot(mbi(1:3000))
grid on
subplot(2,1,2),plot(mbq(1:3000))
grid on
% 加入高斯白噪声
noise_msk=awgn(msk,10,'measured');
figure;subplot(2,1,1);plot(noise_msk(1:3000));grid on
       subplot(2,1,2);plot(noise_msk(3001:6000));grid on
% MSK信号的解调
t=tmin*Ts/0.5:tmin*Ts/0.5:length(a)*Ts;

% 加入50HZ干扰
noise=k_noise*sin(2*pi*fg*t);
noise_msk2=noise_msk+noise;
% 调制信号和载波x(t),y(t)相乘,得到I路,Q路信号
receive_i=noise_msk2.*cos(pi/(2*Ts).*t).*cos(2*pi*fc.*t);
receive_q=noise_msk2.*sin(pi/(2*Ts).*t).*sin(2*pi*fc.*t);
% 画出和本地载波相乘后的I路和Q路的波形图
figure(1);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/50);
%  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:N/2
    sout_i((i-1)*ntb+1)=receive_i((i-1)*ntb+1);
    for j=2:ntb
        sout_i((i-1)*ntb+j)=sout_i((i-1)*ntb+j-1)+receive_i((i-1)*ntb+j);
    end
end

for i=1:N/2
    sout_q((i-1)*ntb+1)=receive_q((i-1)*ntb+1);
    for j=2:ntb
        sout_q((i-1)*ntb+j)=sout_q((i-1)*ntb+j-1)+receive_q((i-1)*ntb+j);
    end
end
ii=1*Ts/0.5:Ts/0.5:N/2*Ts/0.5;
figure;
plot(t,sout_i)
title('sout_i')
figure;
plot(t,sout_q)
title('sout_q')
 % 画出滤波后I路和Q路的信号，并与输入的I路和Q路信号做比较
 ii=1:N/2;
figure;
subplot(2,2,1);
plot(t,receive_i_lowpass)
title('receive_i_lowpass waveform')
subplot(2,2,2),plot(t,receive_q_lowpass)
title('receive_q_lowpass waveform')
subplot(2,2,3),stem(ii,bi)
subplot(2,2,4),stem(ii,bq)       
% 抽样判决，I路信号
for i=1:N/2
    if sout_i((i-1)*ntb+ntb/2)>0
        out_i(i)=1;
    else
        out_i(i)=-1;
    end 
end
% 抽样判决，Q路信号
% 因为I路信号和Q路信号有一个码元周期的延时，导致Q路最后一个码元
% 在指定的长度内不能被恢复出来，程序中对其赋值-1处理。
for i=1:N/2
    if (i-1)*ntb+1.5*ntb>l*ntb
        out_q(i)=-1;
    else
    if sout_q((i-1)*ntb+1.5*ntb)>0
        out_q(i)=-1;
    else
        out_q(i)=1;
    end 
    end
end

figure;
subplot(2,2,1),plot(t,sout_i)
title('receive_i_lowpass waveform')
subplot(2,2,2),plot(t,sout_q)
title('receive_q_lowpass waveform')
subplot(2,2,3),stem(ii,bi);grid on
subplot(2,2,4),stem(ii,bq);grid on
       

% 完成I路信号和Q路信号的并串转换
 for i=0:N/2-1
    decis_data(i*2+1)=out_i(i+1);
    decis_data(i*2+2)=out_q(i+1);
 end
% 差分译码
 for i=2:N
     out_data(1)=a(1);
     if  decis_data(i)==1
         out_data(i)=-decis_data(i-1);
     else
         out_data(i)=decis_data(i-1);
     end
 end
out_data
% 计算误码率
count=0;
for i=1:N-1
    if a(i)==out_data(i)
        count=count+1;
    end
end
pe=(N-count-1)/(N-1);
pe