clear;

N=21;             %累加器的位數；

K=fix((2^N)*0.222);%頻率控制字，即累加的步長；

M=10;             %截取累加器的高10位； 

D=8;              %8bit DAC

g=gcd(K,2^N);

pe=2^N/(g);

n=1:pe;

pp=pe+1;          %頻譜分析點數

add_y=mod(n*K,2^N);  %累加器的輸出表達式；

error=mod(n*K,2^(N-M));

rom_y=cos(2*pi*(add_y-error)/(2^N));     %相位到幅度映射的輸出；

pat = - 1 + 1/ (2^ (D - 1) ):1/ (2^ (D - 1) ):1 - 1/ (2^ (D -1) ); %設定量化區間

codebook = - 1 + 1/ (2^D):1/ (2^(D - 1) ) :1 - 1/ (2^D); %設定量化碼本值

[Pyy,wy]=periodogram(rom_y,[],'onesided',pp);

figure;%figure1 只有相位截斷，沒有量化 figure1

psdplot(Pyy/max(Pyy),wy);%      歸一化顯示出圖；

rom_x=cos(2*pi*add_y/(2^N));

[index,quants2] = quantiz(rom_x,pat,codebook);%量化

[Pzz,wz]=periodogram(quants2,[],'onesided',pp);

figure;%沒有相位截取和ＤＡＣ的非線性等，ｆｉｇｕｒｅ2,只有量化

psdplot(Pzz/max(Pzz),wz);%      歸一化顯示出圖；

[index,quants1] = quantiz(rom_y,pat,codebook);%量化

[Pxx,wx]=periodogram(quants1,[],'onesided',pp);    %頻譜分析；

figure;       %ｆｉｇｕｒｅ3 相位截斷加量化

psdplot(Pxx/max(Pxx),wx);%     歸一化顯示出圖；

grid on