%图1:原始双极性不归零码
%图2:原始双极性不归零码的功率谱
%图3:预编码后的双极性不归零码
%图4:It(串并转换后的上支路，方波)
%图5:Qt(串并转换后的下支路，方波)
%图6:Itt(It与余弦相乘)
%图7:Qtt(Qt与正弦相乘)
%图8:MSK时域信号
%图9:MSK功率谱
%图10:判决后的双极性不归零码
clear all
clc
global dt df t f N
close all
pi=3.1415926
fc=5 %单位:MHz
N=2^13 % 采样点数
L=64 % 每码元的采样点数
M=N/L % 码元数
Rb=2 % 发送码元的信息速率:2Mbps
Tb=1/Rb % 码元宽度:0.5us
dt=Tb/L;
df=1/(dt*N);
T=N*dt % 截短时间
B=N*df/2 % 系统带宽
t=[-T/2+dt/2:dt:T/2]; % 时域横坐标
f=[-B+df/2:df:B]; % 频域横坐标
EP=zeros(size(f));
for ii=1:10;
for j=1:10;
%原始双极性不归零码
b=sign(randn(1,M));
for i=1:L,s(i+[0:M-1]*L)=b;end
figure(1)
plot(t,s);
grid on % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('原始双极性不归零码')
xlabel('t (us)')
ylabel('x(t) (V)')
%原始双极性不归零码的功率谱
P=t2f(s);
P=P.*conj(P)/T;
EP=(EP*(j-1)+P)/j; %RZ功率谱的累计平均
end
Ps=10*log10(EP+eps);% 换算为db值
figure(2)
plot(f,Ps)
grid on % 打开网格
axis([-5*Rb,5*Rb,-50,10]) % 限定范围
title('原始双极性不归零码的功率谱')
xlabel('f (MHz)')
ylabel('H(f) (dB/MHz)')
%预编码
a(1)=b(1);
for i=M:-1:2,a(i)=b(i)*b(i-1);end
for i=1:L,sa(i+[0:M-1]*L)=a;end
figure(3)
plot(t,sa);
grid on % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('预编码后的双极性不归零码')
xlabel('t (us)')
ylabel('sa(t) (V)')
send=sa;
%串并转换
%It
It=zeros(size(t));
for k=0:2*L:N-1;
kk=1:2:2*L;
kkk=1:L;
It(k+kk)=send(k+kkk+L);
It(k+kk+1)=send(k+kkk+L);
end
for k=N:-1:L+1,It(k)=It(k-L);end
figure(4)
plot(t,It);
grid on;
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('It')
xlabel('t (us)')
ylabel('It (V)')
%Qt
Qt=zeros(size(t));
for k=0:2*L:N-1;
kk=1:2:2*L;
kkk=1:L;
Qt(k+kk)=send(k+kkk);
Qt(k+kk+1)=send(k+kkk);
end
figure(5)
plot(t,Qt);
grid on;
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('Qt')
xlabel('t (us)')
ylabel('Qt (V)')
%Itt
Itt=It.*cos(pi*t/2/Tb);
figure(6)
plot(t,Itt);
grid on % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('Itt')
xlabel('t (us)')
ylabel('Itt (V)')
%Qtt
Qtt=Qt.*sin(pi*t/2/Tb);
figure(7)
plot(t,Qtt);
grid on % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('Qtt')
xlabel('t (us)')
ylabel('Qtt (V)')
%MSK时域波形
msk=Itt.*cos(2*pi*fc*t)-Qtt.*sin(2*pi*fc*t);
figure(8)
plot(t,msk);
grid on % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('MSK时域信号')
xlabel('t (us)')
ylabel('msk (V)')
%MSK功率谱
PP=t2f(msk);
Pa=PP.*conj(PP)/T;
EP=(EP*(ii-1)+Pa)/ii; %RZ功率谱的累计平均
end
Pmsk=10*log10(EP+eps);% 换算为db值
figure(9)
plot(f,Pmsk)
grid on % 打开网格
axis([-2*fc,2*fc,-70,0]) % 限定范围
title('MSK功率谱')
xlabel('f (MHz)')
ylabel('H(f) (dB/MHz)')
%接收端
r=msk;
%接收端上支路
RI=r.*cos(2*pi*fc*t);
RI=real(f2t(t2f(RI)));
%接收端下支路
RQ=-r.*sin(2*pi*fc*t);
RQ=real(f2t(t2f(RQ)));
%取样
RIt=RI(2*L:2*L:N); % 上支路取样,每2Tb取样一个值，在整数处取样,-31-32,共M/2=64个
RQt=RQ(L:2*L:N); % 下支路取样,每2Tb取样一个值，在整数+0.5处取样,-31.5-31.5,共M/2=64个
%码型串并转换
Rt=zeros(1,M);
Rt(2:2:M)=RIt(1:M/2);
Rt(1:2:M-1)=RQt(1:M/2);
%判决
Rt=sign(Rt);
clear j; % 清除原复数变量的定义
d(1)=j;
for i=2:M,d(i)=d(i-1)*j;end
e=Rt.*d;
for i=1:2:M,e(i)=imag(e(i));end
f=b-e;
%判决后的双极性不归零码
for i=1:L,sy(i+[0:M-1]*L)=e;end
figure(10)
plot(t,sy);
grid on % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2]) % 限定范围
title('判决后的双极性不归零码')
xlabel('t (us)')
ylabel('sy(t) (V)')
