%matlab實現16QAM調制解調
%調制qammod(a,16)； 解調用qamdemod  畫星座圖modmap('qam',4) 
clear;      %設置參數
M = 16;
fd=200;
fc=4*fd;
fs=4*fc;
rolloff = 0.5;
nsamp = 4;         
num=1000;          %符號個數
% ******* -----產生QAM基帶信號 -----*******************************
x = randint(1,num,M);
y = qammod(x,M);
subplot(2,2,1);
plot(real(y),imag(y),'*'),title('基帶星座圖');
i=real(y);
q=imag(y);
%*************************************************************
% ****************----- 成形濾波 -----*****************************

filtorder = 33;                      % 濾波器長
delay = (filtorder-1)/(nsamp*2);     % 群延時
rrcfilter = rcosine(fd,nsamp*fd,'fir/sqrt',rolloff,delay);  
yi=conv(rrcfilter,upsample(i,nsamp));   
yq=conv(rrcfilter,upsample(q,nsamp));

%************************************************************
% *****************----- 加載波 -----************************
yii=interp(yi,fs/fc);%內插
yqq=interp(yq,fs/fc);
nn=1:length(yii);
yic=yii.*cos(2*pi*nn*fc/fs);% 加載波
yqc=yqq.*sin(2*pi*nn*fc/fs);% 加載波
iq = yic+yqc;

%******************-----進入信道-----******************************
iq = awgn(iq,8,'measured');      %信道中 載波噪聲  
%*********************************************************
%********************----- 去載波 ------************************

ii=iq.*cos(2*pi*nn*fc/fs);
qq=iq.*sin(2*pi*nn*fc/fs);

b=fir1(32,0.2);        
yif=filter(b,1,ii); 
yqf=filter(b,1,qq); 

yiw=2.*downsample(yif,fs/fc,0);
yqw=2.*downsample(yqf,fs/fc,0);
subplot(2,2,3);
plot(1:500,yiw(1+4:500+4),'r',1:500,yi(1:500),'g'),title('紅載波后，綠載波前');
%***********************************************************
% ******************----- 匹配 -----***********************
icc=conv(rrcfilter,yiw);
qcc=conv(rrcfilter,yqw);

icr=downsample(icc,nsamp,0);
qcr=downsample(qcc,nsamp,0);

ic=icr(2+delay*2:length(icr)-delay*2+1);      %去掉冗余
qc=qcr(2+delay*2:length(icr)-delay*2+1);

%********************************************************
%*********************畫出重構星座圖
subplot(2,2,2);
plot(ic,qc,'*'),title('重構星座圖');

%************************ 判 決******************************
iq2=ic+j.*qc;
z = qamdemod(iq2,M);
err=abs(z-x);       %與原信號比較
subplot(2,2,4);
plot(err);
numerr=(length(find(err>0)))/num;       %畫出錯誤符號
numerr
%*********************** 結束 ********************************