%% %% 原理讲解移步链接：贼详细的8PSK调制与解调详细过程 - 子木的文章 - 知乎https://zhuanlan.zhihu.com/p/47258287
%% 8PSK 基本参数
clc
clear all
close all
Fc=70e6;%Fc载波频率
%% 基带信号
data_num=300;
data = randint(1,data_num); %产生1行300列的矩阵，元数为0-1之间的随机数，包括0和1
figure(1)
subplot(211)
plot(data);title('原始基带信号')%以data分量为纵坐标，以元素序号为横坐标，用直线依次连接数据点，绘制曲线
ts=0:1:length(data)-1;
subplot(212)
stem(ts,data);title('原始基带信号');%“有时间向量的话，序列是冲激”
%% 8psk调制
fs_16k=16e3;%假设最开始采样频率是16e3
it=[];qt=[];
sit=[];sqt=[];
for i=1:3:length(data)
I(i) = data(i);
Q(i+1) = data(i+1);
U(i+2) = data(i+2);
if I(i)==0 && Q(i+1)==0 && U(i+2)==0
it=0.924*ones(1,1);%1行1列的全为1的矩阵
qt=0.383*ones(1,1);
elseif I(i)==0 && Q(i+1)==0 && U(i+2)==1
it=0.383*ones(1,1);
qt=0.924*ones(1,1);
elseif I(i)==0 && Q(i+1)==1 && U(i+2)==1
it=-0.383*ones(1,1);
qt=0.924*ones(1,1);
elseif I(i)==0 && Q(i+1)==1 && U(i+2)==0
it=-0.924*ones(1,1);
qt=0.383*ones(1,1);
elseif I(i)==1 && Q(i+1)==1 && U(i+2)==0
it=-0.924*ones(1,1);
qt=-0.383*ones(1,1);
elseif I(i)==1 && Q(i+1)==1 && U(i+2)==1
it=-0.383*ones(1,1);
qt=-0.924*ones(1,1);
elseif I(i)==1 && Q(i+1)==0 && U(i+2)==1
it=0.383*ones(1,1);
qt=-0.924*ones(1,1);
elseif I(i)==1 && Q(i+1)==0 && U(i+2)==0
it=0.924*ones(1,1);
qt=-0.383*ones(1,1);
end
sit=[sit it]; sqt=[sqt qt];
end %映射
figure(2)
subplot(221)
plot(sit);title('i路信号')
subplot(222)
plot(sqt);title('q路信号')
subplot(223)
plot(-fs_16k/2:fs_16k/length(sit):fs_16k/2-fs_16k/length(sit),fftshift(abs(fft(sit))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('I路频谱图')
subplot(224)
plot(-fs_16k/2:fs_16k/length(sqt):fs_16k/2-fs_16k/length(sqt),fftshift(abs(fft(sqt))));
xlabel('Frequency(Hz)');
title('Q路频谱图')
figure(3)
scatter(sit,sqt);%星座图
grid on
%% 16k-160k,插值：先插后滤，抽取：先抽后滤
fs_160k=160e3;
sit_160k=upsample(sit,10);
sqt_160k=upsample(sqt,10);
beta=0.4;%r=0时是矩形
span=length(sit_160k);
sps=10;
h = rcosdesign(beta,span,sps);
sit_160k_f=conv(sit_160k,h,'same');
sqt_160k_f=conv(sqt_160k,h,'same');
figure(4)
subplot(221)
plot(sit_160k_f);title('10倍插值i路')
subplot(222)
plot(sqt_160k_f);title('10倍插值q路')
subplot(223)
plot(-fs_160k/2:fs_160k/length(sit_160k_f):fs_160k/2-fs_160k/length(sit_160k_f),fftshift(abs(fft(sit_160k_f))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('I路频谱图(升余弦)')
subplot(224)
plot(-fs_160k/2:fs_160k/length(sqt_160k_f):fs_160k/2-fs_160k/length(sqt_160k_f),fftshift(abs(fft(sqt_160k_f))));
xlabel('Frequency(Hz)');
title('Q路频谱图(升余弦)')
%% 160k-1600k,插值：先插后滤，抽取：先抽后滤
fs_1600k=1600e3;
sit_1600k=upsample(sit_160k_f,10);
sqt_1600k=upsample(sqt_160k_f,10);
y_de7=fir1(127,1/10);
sit_1600k_f=conv(sit_1600k,y_de7,'same');
sqt_1600k_f=conv(sqt_1600k,y_de7,'same');
figure(5)
subplot(221)
plot(sit_1600k_f);title('100倍插值i路')
subplot(222)
plot(sqt_1600k_f);title('100倍插值q路')
subplot(223)
plot(-fs_1600k/2:fs_1600k/length(sit_1600k_f):fs_1600k/2-fs_1600k/length(sit_1600k_f),fftshift(abs(fft(sit_1600k_f))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('I路频谱图(低通)')
subplot(224)
plot(-fs_1600k/2:fs_1600k/length(sqt_1600k_f):fs_1600k/2-fs_1600k/length(sqt_1600k_f),fftshift(abs(fft(sqt_1600k_f))));
xlabel('Frequency(Hz)');
title('Q路频谱图(低通)');
%% 1600k-16M,插值：先插后滤，抽取：先抽后滤
fs_16M=16e6;
sit_16M=upsample(sit_1600k_f,10);
sqt_16M=upsample(sqt_1600k_f,10);
y_de7=fir1(127,1/10);
sit_16M_f=conv(sit_16M,y_de7,'same');
sqt_16M_f=conv(sqt_16M,y_de7,'same');
figure(6)
subplot(221)
plot(sit_16M_f);title('1000倍插值i路')
subplot(222)
plot(sqt_16M_f);title('1000倍插值q路')
subplot(223)
plot(-fs_16M/2:fs_16M/length(sit_16M_f):fs_16M/2-fs_16M/length(sit_16M_f),fftshift(abs(fft(sit_16M_f))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('I路频谱图(低通)')
subplot(224)
plot(-fs_16M/2:fs_16M/length(sqt_16M_f):fs_16M/2-fs_16M/length(sqt_16M_f),fftshift(abs(fft(sqt_16M_f))));
xlabel('Frequency(Hz)');
title('Q路频谱图(低通)');
%% 16M-25M
fs_25M=25e6;
%resample为信号降采样处理，理解如下：B=resample(x,90,250);
% 采样从250Hz降到90Hz，如果250在前,就是插值从90到250,可以看B的长度,250Hz采样4000个数据等于90hz采样1440个数据,这就是降采样。
sit_25M=resample(sit_16M_f,25,16);
sqt_25M=resample(sqt_16M_f,25,16);
figure(7)
subplot(221)
plot(sit_25M);title('25/16倍插值i路')
subplot(222)
plot(sqt_25M);title('25/16倍插值q路')
subplot(223)
plot(-fs_25M/2:fs_25M/length(sit_25M):fs_25M/2-fs_25M/length(sit_25M),fftshift(abs(fft(sit_25M))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('I路频谱图(低通)')
subplot(224)
plot(-fs_25M/2:fs_25M/length(sqt_25M):fs_25M/2-fs_25M/length(sqt_25M),fftshift(abs(fft(sqt_25M))));
xlabel('Frequency(Hz)');
title('Q路频谱图(低通)')
%% 25M-250M
fs_250M=250e6;
%resample为信号降采样处理，理解如下：B=resample(x,90,250);
% 采样从250Hz降到90Hz，如果250在前,就是插值从90到250,可以看B的长度,250Hz采样4000个数据等于90hz采样1440个数据,这就是降采样。
sit_250M=upsample(sit_25M,10);
sqt_250M=upsample(sqt_25M,10);
y_de7=fir1(127,1/10);
sit_250M_f=conv(sit_250M,y_de7,'same');
sqt_250M_f=conv(sqt_250M,y_de7,'same');
figure(8)
subplot(221)
plot(sit_250M_f);title('10倍插值i路')
subplot(222)
plot(sqt_250M_f);title('10倍插值q路')
subplot(223)
plot(-fs_250M/2:fs_250M/length(sit_250M_f):fs_250M/2-fs_250M/length(sit_250M_f),fftshift(abs(fft(sit_250M_f))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('I路频谱图(低通)')
subplot(224)
plot(-fs_250M/2:fs_250M/length(sqt_250M_f):fs_250M/2-fs_250M/length(sqt_250M_f),fftshift(abs(fft(sqt_250M_f))));
xlabel('Frequency(Hz)');
title('Q路频谱图(低通)')
%% 带通调制
fs_250M=250e6;
T=length(sqt_250M_f);%T=
t=0 : 1/fs_250M :T/fs_250M - 1/fs_250M;%fs = 100
% c=exp(j*2*pi*Fc*t); %载波信号,Fc=70e6
c1=cos(2*pi*Fc*t); %同相载波
c2=-sin(2*pi*Fc*t); %正交载波
psk8 = sit_250M_f.*c1 + sqt_250M_f.*c2;
figure(9),
% subplot(221);plot(t,c);title('载波信号')
subplot(221);plot(t,c1);title('同相载波')
subplot(222);plot(t,c2);title('正交载波')
subplot(223);plot(t,psk8); title('已调信号')
subplot(224);
plot(-fs_250M/2:fs_250M/length(psk8):fs_250M/2-fs_250M/length(psk8),fftshift(abs(fft(psk8))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('已调信号频谱图')
%% 信道
EsN0=4;%单位为db
snr1=10.^(EsN0/10);%单位为比例。Es/N0，信噪比转化为线性值
spnow=norm(psk8).^2/length(psk8);
sigma=sqrt(spnow/(2*snr1));
psk8_channel=psk8+sigma*randn(1,length(psk8));
% bpf_70M=load('bpf_70M.fcf');%两种滤波器调用方式
% BPF_Data=filter(bpf_70M,10,psk8_channel);
bpf70M_Coe=load('bpf70M.mat');%两种滤波器调用方式
BPF_Data=filter(bpf70M_Coe.bpf70M,1,psk8_channel);
figure(18)
subplot(211)
plot(BPF_Data);title('带通滤波之后的已调信号')
subplot(212)
plot(-fs_250M/2:fs_250M/length(BPF_Data):fs_250M/2-fs_250M/length(BPF_Data),fftshift(abs(fft(BPF_Data))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('带通滤波之后的已调信号(带通)')
%% 解调
fs_250M=250e6;
rit_250M_de=2*BPF_Data.*c1;
rqt_250M_de=2*BPF_Data.*c2;
y_de7=fir1(127,1/10);
rit_250M_f=conv(rit_250M_de,y_de7,'same');
rqt_250M_f=conv(rqt_250M_de,y_de7,'same');
figure(10)
subplot(221)
plot(rit_250M_f);title('I路解调信号')
subplot(222)
plot(rqt_250M_f);title('Q路解调信号')
subplot(223)
plot(-fs_250M/2:fs_250M/length(rit_250M_f):fs_250M/2-fs_250M/length(rit_250M_f),fftshift(abs(fft(rit_250M_f))));
xlabel('Frequency(Hz)');
ylabel('Amp');
title('I路解调信号频谱图(低通)')
subplot(224)
plot(-fs_250M/2:fs_250M/length(rqt_250M_f):fs_250M/2-fs_250M/length(rqt_250M_f),fftshift(abs(fft(rqt_250M_f))));
xlabel('
