clear all
clc

%%%%%%%%%%%%%%%%%%%%%%%   不進行衰害補償的系統參數設置   %%%%%%%%%%%%%%%%%%

%各種損傷說明
cng=input('信道噪聲增益選擇：   0（無損失）   0.6（輕微損失）   2（硬損傷） : ');
cdi=input('多徑干擾選擇：   0（無干擾）   1（輕微干擾）   2（硬干擾） :');
fo=input('發射機載頻頻偏：   0   0.01   ：');
po=input('發射機載頻相位補償參數：   0   0.7   0.9   ：');
toper=input('時間補償參數（采樣周期的百分比）：   0   20   30   ：');
so=input('采樣周期偏移量（發射機符號周期低于接收機的單位數）：   0   1   ：');



%%%%%%%%%%%%%%%%%%%%%%%%%%     發射機     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%將字符串轉變為4-PAM序列
str='01234 I wish I were a leader in the future!';
m=letters2pam(str);     %調用letters2pam子函數
N=length(m);

%進行過采樣
M=100;     %過采樣因子
mup=zeros(1,N*M);
mup(1:M:end)=m;     

%用Hamming窗進行脈沖成形
p=hamming(M);     
x=filter(p,1,mup);     %卷積

% figure(1),plotspec(x,1/M)     %作出脈沖成形后的頻譜圖

%調制
t=1/M:1/M:length(x)/M;
fc=20;     %載波頻率
c=cos(2*pi*fc*t);     %載波
r=c.*x;     %調制后的信號

% figure(2),plotspec(r,1/M)     %作出調制后的頻譜圖





%%%%%%%%%%%%%%%%%%%%%%%%%%   引入各種損害參數    %%%%%%%%%%%%%%%%%%%%%%%%%

if cdi<0.5
    mc=[1 0 0];     %無干擾
elseif cdi<1.5
    mc=[1 zeros(1,M) 0.28 zeros(1,2.3*M) 0.11];     %輕微干擾
else mc=[1 zeros(1,M) 0.28 zeros(1,1.8*M) 0.44];     %硬干擾
end
mc=mc/(sqrt(mc*mc'));
dv=filter(mc,1,r);     %經過信道后對發射機信號進行濾波
nv=dv+cng*(randn(size(dv)));     %加入高斯信道噪聲
to=floor(0.01*toper*M);     %采樣周期偏移
rnv=nv(1+to:end);     %指定信號延遲
rt=(1+to)/M:1/M:length(nv)/M;     %采用延遲信息修正時間參數
rM=M+so;     %接收機采樣偏移





%%%%%%%%%%%%%%%%%%%%%%%%%%%%     信道     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%產生一個時變衰落信道
lr=length(nv);
fp=[ones(1,floor(0.2*lr)),0.5*ones(1,floor(lr-0.2*lr))];     %平坦衰落信道
r1=nv.*fp;

%使用自動增益控制對信號進行修正
ds=pow(r1);     %期望的信號平均功率
g=zeros(1,lr);     %初始化增益
g(1)=1;
nr=zeros(1,lr);
mu=0.0003;     %步長
for i=1:lr-1;
    nr(i)=g(i)*r1(i);
    g(i+1)=g(i)-mu*(nr(i)^2-ds);
end
r2=nr;

% figure(3),plotspec(r2,1/M)     %作出經過AGC后的頻譜圖




%%%%%%%%%%%%%%%%%%%%%%%%%%%     接收機     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%采用AM進行方式解調
c2=cos(2*pi*fc*t);     %同步載波進行混頻
x2=r.*c2;     

% figure(4),plotspec(x2,1/M)     %作出解調后的頻譜圖

fl=50;     %LPF中沖激響應的項數
fbe=[0 0.5 0.6 1];
damps=[1 1 0 0];     
b=remez(fl,fbe,damps);     %產生LPF沖激響應
x3=2*filter(b,1,x2);     %進行低通濾波

% figure(5),plotspec(x3,1/M)     %作出經過LPF后的頻譜圖

% l=0.5*fl+0.5*M+0.5*M;     %初始時延l對應于LPF、相關濾波器和信號周期長度的一半
% ul=floor((length(x3)-l-1)/(4*rM));     %作出眼圖
% figure(6),plot(reshape(x3(125:ul*4*rM+124),4*rM,ul))

%與脈沖信號進行相關操作以加強脈沖序列的波峰
y=filter(fliplr(p)/(pow(p)*M),1,x3);     %對接收信號脈沖進行卷積、歸一化

% figure(7),plot(reshape(y(125:ul*4*rM+124),4*rM,ul))      %作出眼圖

%對采樣信號及增量設置延遲
z=y(0.5*fl+M:M:end);     %下采樣

figure(8),plot([1:length(z)],z,'.')     %軟判決后的星座圖

%判決裝置和匹配特性
mprime=quantalph(z,[-3,-1,1,3])';      %量化
cluster_variance=(mprime-z)*(mprime-z)'/length(mprime);     
lmp=length(mprime);
percentage_symbol_error=100*sum(abs(sign(mprime-m(1:lmp))))/lmp     %誤碼率

%重建信號
reconstructed_message=pam2letters(mprime);