nrzrz.m

双极性不归零码、双极性归零码和差分码的时域及功率谱密度图的绘制

%双极性不归零码、双极性归零码和差分码的时域及功率谱密度图
close all
clear all
%采样点数的设置
k=input('采样点数N=2^k，请输入k的值，默认为14:');
if isempty(k)
	k=14;
end
%每码元采样数的设置
L=input('请输入每码元的采样数，要求必须为2的阶乘，默认为64:');
if isempty(L)
    L=64;
end
N=2^k;
M=N/L;%M为码元个数
dt=1/L;%时域采样间隔
T=N*dt;%时域截断区间
df=1.0/T;%频域采样间隔
Bs=N*df/2;%频域截断区间
t=linspace(-T/2,T/2,N);%产生时域采样点
f=linspace(-Bs,Bs,N);%产生频域采样点
EP1=zeros(size(f));
EP2=zeros(size(f));
EP3=zeros(size(f));
for x=1:20%取20次样值
    a=round(rand(1,M));%产生一个长度为M的随机序列a，0和1等概出现
    nrz=zeros(L,M);%产生一个L行M列的nrz矩阵，初始化为全0矩阵
    rz=zeros(L,M);%产生一个L行M列的rz矩阵，初始化为全0矩阵
    cfm=zeros(L,M);%产生一个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
        %产生差分码，令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