clc;
clear ;
close all;
Fc=2; % 载波频率
Ts=1; % 采样周期
Fs=1/Ts; % 采样频率
bl=Fs; % 低通滤波器带宽
n=32; % 抽样点数
N=500; % 码元数
dt=Ts/Fc/n; % 采样间隔
T=N*Ts; % 信号长度
t=0:dt:T-dt; % 时间向量
Lt=length(t); % 时间向量长度
EbNo = 3:0.5:12; % 设定EbNo范围
EsN0 = 10.^(EbNo/10);
RA = zeros(1,length(EbNo)); %#ok<PREALL> % 初始化误码率向量
%%%%****产生二进制信源****%%%%
a=randn(1,N); % 产生N个随机数
d=sign(a); % 将大于0的变为1小于0的变为-1
d1=sigexpand(d,Fc*n); % 将序列d周期变为Ts(d中每个数之间插Fc*n-1个0)
wt=ones(1,Fc*n); % 产生宽度为Ts的矩形窗
D=conv(d1,wt); % 卷积产生基带信号
figure(1);
subplot(3,1,1);
plot(t,d1(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('t(s)');
ylabel('A');
title('基带信号');
grid;
subplot(3,1,2);
plot(t,D(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('t(s)');
ylabel('A');
title('基带信号时域图');
grid;
[f,Df]=T2F(t,D(1:Lt));% 进行傅里叶变换
subplot(3,1,3);
plot(f,10*log10(abs(Df).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('基带信号频谱图');
grid;
figure(6);
subplot(2,2,1);
plot(t,D(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('基带信号时域图');
grid;
subplot(2,2,2);
plot(f,10*log10(abs(Df).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频密度(dB/Hz)');
title('基带信号频谱图');
grid;
% 串并转换
i=zeros(1,N); % 建立一个1行N列的矩阵
q=zeros(1,N);
for j=1:N
if rem(j,2)==1 % j除以2余数是否为1
i((j+1)/2)=d(j);
else
q(j/2)=d(j);
end
end
%%%=======I路信号=======%%%
di=sigexpand(i,2*Fc*n);
wt1=ones(1,2*Fc*n);
I=conv(di,wt1);
figure(2);
subplot(2,3,1);
plot(t,I(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('I路基带信号时域图');
grid;
[f1,IF]=T2F(t,I(1:Lt));
figure(2);
subplot(2,3,4);
plot(f1,10*log10(abs(IF).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('I路基带信号频谱图');
grid;
%%%=======Q路信号========%%%
dq=sigexpand(q,2*Fc*n);
wt1=ones(1,2*Fc*n);
Q=conv(dq,wt1);
Qdly=[-ones(1,n*Fc),Q(1:end-n*Fc)]; %进行延时，在前面添-1
figure(2);
subplot(2,3,2);
plot(t,Qdly(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('Q路基带信号时域图');
grid;
[f2,QF]=T2F(t,Qdly(1:Lt));
figure(2);
subplot(2,3,5);
plot(f2,10*log10(abs(QF).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('Q路基带信号频谱图');
grid;
figure(5);
subplot(2,4,1);
plot(t,I(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('I路基带信号时域图');
grid;
subplot(2,4,5);
plot(t,Qdly(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('Q路基带信号时域图');
grid;
% 载波
y1=cos(2*pi*Fc*t);
y2=sin(2*pi*Fc*t);
figure(2);
subplot(2,3,3);
plot(t,y1);
axis([0,7,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('载波信号时域图');
grid;
[f3,y1f]=T2F(t,y1);
figure(2);
subplot(2,3,6);
plot(f3,10*log10(abs(y1f).^2/T));
% p=2/T *10*log10(abs(h1tf)为求功率谱的公式
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('载波信号频谱图');
grid;
% 生成OQPSK信号
QI=I(1:Lt).* y1; %上下支路分别调制
QQ=Qdly(1:Lt).* y2;
figure(3);
subplot(2,2,1);
plot(t,QI);
axis([0,7,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('I路频带信号时域图');
grid;
[f4,QIf]=T2F(t,QI);
figure(3);
subplot(2,2,3);
plot(f4,10*log10(abs(QIf).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('I路频带信号功率谱图');
grid;
figure(3);
subplot(2,2,2);
plot(t,QQ);
axis([0,7,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('Q路频带信号时域图');
grid;
[f5,QQf]=T2F(t,QQ);
figure(3);
subplot(2,2,4);
plot(f5,10*log10(abs(QQf).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('Q路基带信号频谱图');
grid;
Oqpsk=QI+QQ; %两路信号相加得到调制后的信号
figure(4);
subplot(2,2,1);
plot(t,Oqpsk);
axis([0,7,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('已调信号时域图');
grid;
[f6,Oqpskf]=T2F(t,Oqpsk);
subplot(2,2,3);
plot(f6,10*log10(abs(Oqpskf).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('已调信号频谱图');
grid;
%%%====信道加入高斯白噪声进行接收解调====%%%
% 产生高斯白噪声
p1=-22;
Noise = wgn(1,Lt,p1); % 生成1*Lt个高斯白噪声，功率为-22dBW（分贝瓦）
% 接收信号Y
Y = Oqpsk + Noise;
figure(4);
subplot(2,2,2);
plot(t,Y);
axis([0,7,-2,2]);
xlabel('时间(S)');
ylabel('幅度');
title('接收信号时域图 ');
grid;
[f7,Yf]=T2F(t,Y);
subplot(2,2,4);
plot(f7,10*log10(abs(Yf).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('功率谱密度(dB/Hz)');
title('接收信号频谱图 Pn=-22dB');
grid;
% 相干解调
% 通过乘法器1
RI = Y .* y1; % 按元素相乘
figure(5);
subplot(2,4,2);
plot(t,RI);
axis([0,14,-2,2]);
xlabel('时间(S)');
ylabel('幅度');
title('通过乘法器I路信号时域图');
grid;
%通过低通滤波器
[f8,RIF]=T2F(t,RI);
[t1,RIL]=lpf(f8,RIF,bl);
subplot(2,4,3);
plot(t1,RIL)
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('通过低通滤波器I路信号时域图');
grid;
%抽样判决
cyI=RIL(Fc*n:2*Fc*n:end);
pjI=sign(cyI); %判决
di1=sigexpand(pjI,2*Fc*n);
I1=conv(di1,wt1); %得到上支路信号
I1dly=[-ones(1,n*Fc),I1(1:end-n*Fc)]; %上支路信号延时
subplot(2,4,4);
plot(t,I1dly(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('I路抽样判决后信号时域图');
grid;
% 通过乘法器2
RQ = Y .* y2;
subplot(2,4,6);
plot(t,RQ);
axis([0,14,-2,2]);
xlabel('时间(S)');
ylabel('幅度');
title('Q路通过乘法器信号时域图');
grid;
%通过低通滤波器
[f9,RQF]=T2F(t,RQ);
[t2,RQl]=lpf(f9,RQF,bl);
subplot(2,4,7);
plot(t2,RQl)
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('Q路信号通过低通滤波器时域图');
grid;
%抽样判决
cyQ=RQl(Fc*n:2*Fc*n:end);
pjQ=sign(cyQ); %判决
dq1=sigexpand(pjQ,2*Fc*n);
Q1 = conv(dq1,wt1);
subplot(2,4,8);
plot(t,Q1(1:Lt))
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('Q路抽样判决后信号时域图');
grid;
pjIdly=I1dly(Fc*n:2*Fc*n:end);
% 并串转换
d1=zeros(1,N);
for s=1:N/2
d1(2*s-1)=pjIdly(s);
d1(2*s)=pjQ(s);
end
%输出解调信号
rd1=sigexpand(d1,Fc*n);
R=conv(rd1,wt);
QR=[R(2*n*Fc+1:N*n*Fc),ones(1,2*n*Fc)]; % 去除延时
figure(6);
subplot(2,2,3);
plot(t,QR(1:Lt));
axis([0,14,-1.5,1.5]);
xlabel('时间(S)');
ylabel('幅度');
title('解调信号时域图 Pn=-22dB');
grid;
[f10,QRf]=T2F(t,R(1:Lt));
subplot(2,2,4);
plot(f10,10*log10(abs(QRf).^2/T));
axis([-8,8,-50,15]);
xlabel('频率(Hz)');
ylabel('频谱密度(dB/Hz)');
title('解调信号频谱图 Pn=-22dB');
grid;
%眼图
eyediagram(Y,32,2,8);%eyediagram (x, n, period, offset)增加了偏置因子。
%eyediagram函数认为信号的第( offset+1)个采样值之后每n个值为一周期，
%且该周期为period的整数倍。Off set必须是非负整数，其范围是{0, n-1 }
x=I(1:Fc*N*n)+1i*Qdly(1:Fc*N*n);
%星座图
scatterplot(x);
axis([-1.5
