%數字通信QAM仿真
clear
%產生基帶符號流
M=16;       %進制 16QAM
N=40;       %信源比特流長度,為保證對齊,選擇是log2(M)的整倍數
bitstream=zeros(1,N);
symbolstream=zeros(1,N/log2(M));
for n1=1:N
    bitstream(n1)=binornd(1,0.5,1,1);
end
%格雷碼編碼
for n1=1:length(symbolstream)
    pack=[bitstream(4*(n1-1)+1) bitstream(4*(n1-1)+2) bitstream(4*(n1-1)+3) bitstream(4*(n1-1)+4)];
    if  pack==[0 0 0 0]
        symbolstream(n1)=1+j;
    elseif  pack==[0 0 1 0]
        symbolstream(n1)=2+j;
    elseif  pack==[0 0 1 1]
        symbolstream(n1)=2+2j;
    elseif  pack==[0 0 0 1]
        symbolstream(n1)=1+2j;
    elseif  pack==[0 1 0 1]
        symbolstream(n1)=-1+2j;
    elseif  pack==[0 1 0 0]
        symbolstream(n1)=-1+j;
    elseif  pack==[0 1 1 1]
        symbolstream(n1)=-2+2j;
    elseif  pack==[0 1 1 0]
        symbolstream(n1)=-2+j;
    elseif  pack==[1 1 1 0]
        symbolstream(n1)=-2-j;
    elseif  pack==[1 1 1 1]
        symbolstream(n1)=-2-2j;
    elseif  pack==[1 1 0 1]
        symbolstream(n1)=-1-2j;
    elseif  pack==[1 1 0 0]
        symbolstream(n1)=-1-j;
    elseif  pack==[1 0 0 0]
        symbolstream(n1)=1-j;
    elseif  pack==[1 0 1 0]
        symbolstream(n1)=2-j;
    elseif  pack==[1 0 0 1]
        symbolstream(n1)=1-2j;
    else symbolstream(n1)=2-2j;
    end
end
figure(1)
hold on
for n1=-2:2
    for n2=-2:2
        if ~((n1==0)|(n2==0)) 
            scatter(n1,n2);
        end
    end
end
axis([-3 3 -3 3]);
hold off
title('16QAM星座圖');
xlabel('同相支路')
ylabel('正交支路')
%產生載波
fc=10;           %載波頻率
fs=100;          %采樣頻率
T=1;             %每符號持續時間
N_samples=T*fs;  %每符號內的采樣點數
t=0:T/N_samples:(T-T/N_samples);
carrier=exp(j*2*pi*fc*t);
%產生信號脈沖g(t)
%gt=ones(1,length(carrier));
gt=1-cos(2*pi*(1/T)*t);

%生成調制信號S(t)
St_complex=zeros(1,length(carrier)*length(symbolstream));
for n1=1:length(symbolstream)
    St_complex((N_samples*(n1-1)+1):(N_samples*(n1-1)+N_samples))=(symbolstream(n1)*carrier).*gt;
end
figure(2)
St_real=real(St_complex);
plot(St_real)
title('QAM仿真波形圖 載波10Hz g(t)為升余弦脈沖');
xlabel('采樣點')
ylabel('幅度')
%求出功率譜,這種功率譜不是按照公式直接計算的,而是對所產生的一段時間內的信號的采樣點進行FFT所得到的
spectrum=(real(fft(St_real,10*length(St_real)))).^2;
S_spectrum=spectrum(1:length(spectrum)/2);
figure(3)
F=0:fs/(2*length(S_spectrum)):fs/2-fs/(2*length(S_spectrum));
plot(F,10*log10(S_spectrum))
axis([0 50 -70 60]);
title('QAM信號功率譜密度')
xlabel('頻率/Hz')
ylabel('功率')
figure(4)
index=400:1600;
plot(F(index),10*log10(S_spectrum(index)));
axis([4 16 0 60]);
title('QAM信號功率譜密度在載頻附近的放大圖')
xlabel('頻率/Hz')
ylabel('功率')