%雙極性不歸零碼、雙極性歸零碼和差分碼的時域及功率譜密度圖
close all
clear all
%采樣點數(shù)的設置
k=input('采樣點數(shù)N=2^k，請輸入k的值，默認為14:');
if isempty(k)
	k=14;
end
%每碼元采樣數(shù)的設置
L=input('請輸入每碼元的采樣數(shù)，要求必須為2的階乘，默認為64:');
if isempty(L)
    L=64;
end
N=2^k;
M=N/L;%M為碼元個數(shù)
dt=1/L;%時域采樣間隔
T=N*dt;%時域截斷區(qū)間
df=1.0/T;%頻域采樣間隔
Bs=N*df/2;%頻域截斷區(qū)間
t=linspace(-T/2,T/2,N);%產(chǎn)生時域采樣點
f=linspace(-Bs,Bs,N);%產(chǎn)生頻域采樣點
EP1=zeros(size(f));
EP2=zeros(size(f));
EP3=zeros(size(f));
for x=1:20%取20次樣值
    a=round(rand(1,M));%產(chǎn)生一個長度為M的隨機序列a，0和1等概出現(xiàn)
    nrz=zeros(L,M);%產(chǎn)生一個L行M列的nrz矩陣，初始化為全0矩陣
    rz=zeros(L,M);%產(chǎn)生一個L行M列的rz矩陣，初始化為全0矩陣
    cfm=zeros(L,M);%產(chǎn)生一個L行M列的cfm矩陣，初始化為全0矩陣
    for i=1:M 
        if a(i)==1
            nrz(:,i)=1;%使nrz矩陣第i列全部元素都為1
            rz(1:L/2,i)=1;%使rz矩陣第i列前L/2個元素為1
        else
            nrz(:,i)=-1;%使nrz矩陣第i列全部元素都為-1
            rz(1:L/2,i)=-1;%使rz矩陣第i列前L/2個元素為-1
        end
    end
    for i=2:M
        %產(chǎn)生差分碼，令cfm矩陣第i列全部元素都等于其第i-1列的第一個元素與矩陣a第i個元素相加后模2的結果
        cfm(:,i)=mod((a(i)+cfm(1,i-1)),2);
    end;
    %分別重排nrz、rz和cfm矩陣為1行N列的矩陣
    nrz=reshape(nrz,1,N);
    rz=reshape(rz,1,N);
    cfm=reshape(cfm,1,N);
    %做傅立葉變換并算出功率譜密度
    NRZ=t2f(nrz,dt);
    P1=NRZ.*conj(NRZ)/T;
    RZ=t2f(rz,dt);
    P2=RZ.*conj(RZ)/T;
    CFM=t2f(cfm,dt);
    P3=CFM.*conj(CFM)/T;
    %求功率譜密度的均值
    EP1=(EP1*(x-1)+P1)/x;
    EP2=(EP2*(x-1)+P2)/x;
    EP3=(EP3*(x-1)+P3)/x;
end
figure(1)%開啟一個編號為1的繪圖窗口
%設置窗口1左上角的位置在距屏幕左側0像素、下側400像素的地方，長為340像素，寬為300像素
set(1,'position',[0,400,340,300])
plot(t,nrz)%畫雙極性不歸零碼的時域圖
axis([-10,10,min(nrz)-0.1,max(nrz)+0.1])
title('雙極性不歸零碼')
xlabel('t(ms)')
ylabel('nrz(t)')
grid on
figure(2)%開啟一個編號為2的繪圖窗口
%設置窗口2左上角的位置在距屏幕左側340像素、下側400像素的地方，長為340像素，寬為300像素
set(2,'position',[340,400,340,300])
plot(t,rz)%畫雙極性歸零碼的時域圖
axis([-10,10,min(rz)-0.1,max(rz)+0.1])
title('雙極性歸零碼')
xlabel('t(ms)')
ylabel('rz(t)')
grid on
figure(3)%開啟一個編號為3的繪圖窗口
%設置窗口3左上角的位置在距屏幕左側680像素、下側400像素的地方，長為340像素，寬為300像素
set(3,'position',[680,400,340,300])
plot(t,cfm)%畫差分碼的時域圖
axis([-10,10,min(cfm),max(cfm)+0.1])
title('差分碼')
xlabel('t(ms)')
ylabel('cfm(t)')
grid on
figure(4)%開啟一個編號為4的繪圖窗口
P1B=30+10*log10(EP1+eps);%將功率譜密度的單位轉換成dB
%設置窗口4左上角的位置在距屏幕左側0像素、下側50像素的地方，長為340像素，寬為300像素
set(4,'position',[0,50,340,300])
plot(f,P1B)%畫雙極性不歸零碼的功率譜密度圖
axis([-5,5,-50,50])
title('雙極性不歸零碼功率譜密度圖')
xlabel('f(kHz)')
ylabel('P1(f)')
grid on
figure(5)%開啟一個編號為5的繪圖窗口
P2B=30+10*log10(EP2+eps);%將功率譜密度的單位轉換成dB
%設置窗口5左上角的位置在距屏幕左側340像素、下側50像素的地方，長為340像素，寬為300像素
set(5,'position',[340,50,340,300])
plot(f,P2B)%畫雙極性歸零碼的時域圖
axis([-5,5,-50,50])
title('雙極性歸零碼功率譜密度圖')
xlabel('f(kHz)')
ylabel('P2(f)')
grid on
figure(6)%開啟一個編號為6的繪圖窗口
P3B=30+10*log10(EP3+eps);%將功率譜密度的單位轉換成dB
%設置窗口6左上角的位置在距屏幕左側680像素、下側50像素的地方，長為340像素，寬為300像素
set(6,'position',[680,50,340,300])
plot(f,P3B)%畫差分碼的時域圖
axis([-5,5,-50,50])
title('差分碼功率譜密度圖')
xlabel('f(kHz)')
ylabel('P3(f)')
grid on