%16QAM中頻調制
%升余弦窗成形，滾降系數0.35，符號率1MSybol/s,中頻頻率21MHz
%給出各級濾波器的系數，時域/頻域響應以及信號經過各濾波器的時域/頻域圖

%*************************產生QAM基帶信號*************************%
%產生偽隨機序列PN

N=500; %二進制數據長度，長度不足會造成星座點缺失
x=randint(1,N,2); 

%數據分組,串并變換
x1=x(1:2);
x2=x(3:4); %完成第一組轉換
for i=1:(N/4-1)  %完成所有點的轉換  
    x1=[x1 (x(i*4+1:i*4+2))];
    x2=[x2 (x(i*4+3:i*4+4))];
end

%二-十進制轉換（00-0,01-1,10-2,11-3)
xi=x1(1)*2+x1(2);
xq=x2(1)*2+x2(2); %完成第一組轉換
n=length(x1);
for i=1:n/2-1    %完成所有點的轉換  
    xi=[xi (x1(i*2+1)*2+x1(i*2+2))];
    xq=[xq (x2(i*2+1)*2+x2(i*2+2))];
end

%信號映射
for i=1:n/2
    switch(xi(i))
        case 0
            xi(i)=-3;
        case 1
            xi(i)=-1;
        case 2
            xi(i)=1;
        case 3
            xi(i)=3;
    end
    
    switch(xq(i))
        case 0
            xq(i)=-3;
        case 1
            xq(i)=-1;
        case 2
            xq(i)=1;
        case 3 
            xq(i)=3;
    end
end

%繪制基帶信號圖
figure(1);
plot(xi,xq,'*'),title('16AQM基帶信號星座圖');
axis([-4,4,-4,4]);


%***************************脈沖成形*************************************%

%升余弦窗，符號率1MSym/s，設采樣頻率Fs=4MHz
Fs=4000000; %采樣頻率設為符號率4倍
t=1/Fs;     %周期
R=0.35;     %滾降系數0.35
N_T=7;      %窗長，設置為標量7
RATE=4;     %采樣率
p=rcosfir(R,N_T,RATE,t,'sqrt'); %Design a raised cosine FIR filter，設計升余弦窗

figure(2);  %繪制升余弦窗時域圖
stem(p),title('升余弦窗時域圖');
figure(3);
freqz(p),title('升余弦窗頻域圖');
 
%I路和Q路信號成形
xi=upsample(xi,4); %信號內插
xq=upsample(xq,4);
yi=conv(p,xi);  %I路信號與升余弦窗卷積
yq=conv(p,xq);  %Q路信號與升余弦窗卷積
figure(4);      %繪制成形信號圖
subplot(2,1,1);
plot(yi),title('I路成形信號');
subplot(2,1,2);
plot(yq),title('Q路成形信號');
figure(5);
plot(db(freqz(yi)/max(abs(freqz(yi))))),title('I路成形信號頻譜'),xlabel('頻率/Hz'),ylabel('幅度/db');
figure(6);
plot(db(freqz(yq)/max(abs(freqz(yq))))),title('Q路成形信號頻譜'),xlabel('頻率/Hz'),ylabel('幅度/db');

%********************************HBF********************************%
%設計半帶濾波器，階數為15。
B=firhalfband(14,0.01,'dev'); %15階半帶濾波器時域響應
figure(7);
stem(B);
title('半帶濾波器時域圖');
figure(8);
plot(db(abs(freqz(B))));
title('半帶濾波器頻域圖');


%I、Q兩路信號分別與HBF卷積
%yi=upsample(yi,4);
%yq=upsample(yq,4);
yih=conv(B,yi);
yqh=conv(B,yq);
figure(9);   %繪制經過HBF后的信號
subplot(2,1,1);
plot(yih),title('HBF濾波后I路信號'); %時域信號
subplot(2,1,2);
plot(yqh),title('HBF濾波后Q路信號');
figure(10);
%freqz(yih);
plot(db(freqz(yih)/max(abs(freqz(yih))))),title('HBF濾波后I路信號頻譜'),xlabel('頻率/Hz'),ylabel('幅度/db') %頻域信號
figure(11);
plot(db(freqz(yqh)/max(abs(freqz(yqh))))),title('HBF濾波后Q路信號頻譜'),xlabel('頻率/Hz'),ylabel('幅度/db');

%*****************************CIC******************************%
%設計CIC
Rc=3; %interpolation factor,內插系數
Mc=1; %defferential delay,延時
N=5;  %number of sections 單元數量
Hm=mfilt.cicinterp(Rc,Mc,N); 
%figure(12);
%stem(Hm);
%I、Q兩路信號分別經過CIC濾波
y_fi=filter(Hm,yih);
y_fq=filter(Hm,yqh);
yih=double(yih); yic=double(y_fi); yic=yic/max(abs(yic));%數據類型轉換，并歸一化
yqh=double(yqh); yqc=double(y_fq); yqc=yqc/max(abs(yqc));
figure(13);
subplot(2,1,1);
plot(yic),title('I路信號經過CIC后的時域響應'),axis([1,900,-1.5,1.5]);
subplot(2,1,2);
plot(yqc),title('Q路信號經過CIC后的時域響應'),axis([1,900,-1.5,1.5]);
figure(14);
%freqz(yic);
plot(db(freqz(yic)/max(abs(freqz(yic))))),title('I路信號經過CIC后的頻域響應'),xlabel('頻率/Hz'),ylabel('幅度/db');
figure(15);
plot(db(freqz(yqc)/max(abs(freqz(yqc))))),title('Q路信號經過CIC后的頻域響應'),xlabel('頻率/Hz'),ylabel('幅度/db');


%****************************中頻信號*****************************%
%中頻21MHz
fc=21000000;
f0=2000000*32;
l=1:length(yic);   %信號序列長度
%F=cos(2*pi*l*fc/f0);  %載波信號
%figure(16);
%plot(abs(freqz(F)));
%figure(17);
%plot(F);
Yi=yic.*cos(2*pi*l*fc/f0);
Yq=yqc.*sin(2*pi*l*fc/f0);
figure(18);
subplot(2,1,1);
plot(Yi/max(abs(Yi))),title('I路中頻信號時域圖');
subplot(2,1,2);
plot(Yq/max(abs(Yq))),title('q路中頻信號時域圖');
figure(19)
plot(db(fft(Yi)/max(abs(fft(Yi))))),title('I路中頻信號頻域圖');
figure(20);
plot(db(fft(Yq)/max(abs(fft(Yq))))),title('Q路中頻信號頻域圖');
%I、Q兩路信號合成
Y=Yi+Yq;
figure(21);
plot(Y),title('中頻輸出信號');
figure(23);
plot(abs(fft(Y))),title('中頻信號頻域圖');