%file :8pskmodem.m
%产生二进制数字序列a，长度为no_seq
clear all
clc
no_seq=15000;
Tb=1;
fs=8; %采样频率
ts=1/fs;
N=Tb*fs; %一个比特内的样点数
a=randint(1,no_seq); %产生no_seq长度的二进制序列
No_sample=length(a)*N; %样点总数
%内插及显示原码序列
a2=kron(a,[ones(1,N)]);
t=(0:No_sample-1)/8*ts;
figure(1);
subplot(311);
plot(t*1000,a2);
title('原码序列显示');
xlabel('time/ms');
ylabel('scope/V');
axis([0 1000 0 1.5]);
grid;
%产生8psk信号
%产生:111 pi/8 110 3pi/8 010 5pi/8 011 7pi/8
%001 9pi/8 000 11pi/8 100 13pi/8 101 15pi/8
mapping=[11*pi/8 9*pi/8 5*pi/8 7*pi/8 13*pi/8 15*pi/8 3*pi/8 1*pi/8 ];
% mapping=[0 pi/4 3*pi/4 pi/2 7*pi/4 3*pi/2 pi 5*pi/4];
rpsk=[];
for k=1:3:no_seq
index=0;
for jj=k:k+2
index=2*index+a(jj);
end
index=index+1;
theta=mapping(index);
rpsk=[rpsk exp(j*theta)];%调制后的8PSK信号
end
%插值补零
crpsk=upsample(rpsk,N);
%成型滤波滤波器参数及产生
R=0.5; %
N_T=4; %
rate=8; %
T=1;
hh=rcosfir(R,N_T,rate,T,'sqrt');
%滤波器归一化，使Es=1
hh=hh./(max(hh))*(1-R+4*R/pi);
subplot(312);
plot(t(1:65)*1000,hh);
title('hh成型滤波器波形');
xlabel('t/ms');
grid;
axis([0 1000 -1.5 1.5]);
%成型滤波
sig = conv(crpsk,hh);
% sig = filter(hh,1,upsample(sig,N));
%画出开始部分若干个码元内成型滤波后的RPSK波形
sig2=sig((length(hh)+1)/2:end-(length(hh)-1)/2);
subplot(313);
plot(t(1:40000)*1000,sig2);
title('sig 成型滤波波形');
xlabel('t/ms');
grid;
axis([0 1000 -1.5 1.5]);
%功率谱
n=N;
df=fs/(no_seq/3);
ff = [0:df:df*((no_seq/3)-1)]-fs/2;
% sig2=sig((length(hh)+1)/2:end-(length(hh)-1)/2);
nsig = reshape(sig2,N*no_seq/3/n,n)';
for mm=1:n
SIG(mm,:) = (abs(fft(nsig(mm,:))/fs)).^2;
end
NSIG = sum(SIG,1);
NSIG = NSIG/max(NSIG);
figure(3);
semilogy(ff,fftshift(NSIG),'r-')
axis([-4 4 1e-6 1]);grid on;
ylabel('功率谱/V');
xlabel('频率/Hz');
title('RPSK功率谱');
grid;
%作出接收信号的眼图和星座图
h1=eyediagram(sig,4*fs,fs);
% h4=scatterplot(sig2);grid on
% QPSK经过有噪信道，解调，误码率统计
snr_in_db1=0:12;
Eb=1/3;
for jj = 1:length(snr_in_db1)
snr=10^(snr_in_db1(jj)/10);
Nsigma = sqrt(fs*Eb/snr/2);
resig = sig + Nsigma*(randn(size(sig))+j*randn(size(sig)));
% resig = sig +0;
%匹配滤波解调
resig2 = conv(hh,resig);
%抽样
desig = resig2(length(hh):N:end-(length(hh)));
% 误符号计算误码率
de_sig = reshape(desig,length(desig),1);
[mm nn] = max(real((de_sig*exp(-j.*mapping)).'));
nn = nn - 1;
%解映射
deb_sig = zeros(1,length(a));
deb_sig(1:3:end) = floor(nn./4);
deb_sig(2:3:end) = floor((nn - floor(nn./4).*4)./2);
deb_sig(3:3:end) = rem(nn,2);
error_rate(jj)=sum(xor(deb_sig,a))/length(a);
end
snr11=10.^(snr_in_db1./10);
h2=scatterplot(desig);
grid on;
p = 1/3*erfc(sqrt(3.*snr11).*sin(pi/8)); % 理论曲线
figure(5);
semilogy(snr_in_db1,error_rate,'r*','LineWidth',2);hold on;
semilogy(snr_in_db1,p,'b','LineWidth',2);
grid on;
