rz_dpsk_spectrum.asv

光通信中非归零型差分相位调制信号建模与仿真分析

%用于产生7级伪随机码序列
clear
pack
input7=[1 0 0 1 1 0 1];
%input7=[1 0 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 1 1 0 1];

p7=[0 0 0 1 0 0 1];
prbs7=[ones(1,120) input7];
for ii=1:120;
    a=input7;                              %用于保存输入
    b=input7.*p7;                          %产生下位的输入
    c=sum(b);
    d=mod(c,2);
    prbs7(121-ii)=d;
    input7=[d a(1) a(2) a(3) a(4) a(5) a(6)];
end
a=0;                                       %论证1 0出现概率是相等
b=0;
for ii=1:127;
    if prbs7(ii)==1
        a=a+1;
    else
        b=b+1;
    end
end
a=a/127;
b=b/127;
%脉冲初始情况的设定，
fs=2e10;                                    %输入电信号时钟，其输出光脉冲时钟为电信号两倍
ts=1/(2*fs);;                                 %脉冲时间
p0=1;                                       %设定输入光功率
f0=1.931e14;                                %光波频率
%用于产生各种波形,首先产生一个脉冲的离散点，然后得到全部序列的离散点
n=8;                                     %单个脉冲的取样点

%**************************************************************************
%profile为单个脉冲的包络函数
%t为产生全部脉冲的离散点，注意值域起始点的选定,占空比不同则起始点不同, 否则脉冲显示不完整        

t=linspace(0.5*ts,1.5*ts,n); 
profile=p0*(cos((pi/2)*sin(2*pi*fs*t)));                  %占空比为1/3
% t=linspace(0.0*ts,1.0*ts,n);                
% profile=p0*(sin((pi/2)*sin(2*pi*fs*t)))*exp(i*pi/2);      %占空比为2/3，即CSRZ-DPSK
% t=linspace(0.75*ts,1.75*ts,n);
% profile=p0*(cos((pi/4)*(sin(2*pi*2*fs*t)-1)))*exp(i*pi/4);    %占空比为1/2，注意此时调制频率为2*fs，不同于以上两种。

%**************************************************************************
%完成整个序列取样点的采集(加载数据信号到相位中)
m=length(prbs7);                           %input7为特定序列，换成prbs7则为伪随机序列

k=m*n;                                      %全部取样点的个数
sequence=zeros(1,k);                              %产生取样点矩阵                              
for ii=1:m;
   if sequence_c(ii)==1;                        %input7为特定序列，换成prbs7则为伪随机序列
       for jj=1:n;
      sequence(((ii-1)*n+jj))=profile(jj);             %数据为1则相位相同
       end
   else
       for jj=1:n;
      sequence(((ii-1)*n+jj))=exp(i*pi)*profile(jj);   %数据为0则相位相反
       end
   end
end
kk=length(sequence);
%产生载波信号
wt=linspace(0,(ts*m),kk);
sequence_c=sequence.*cos(2*pi*f0*wt);
%**************************************************************************
% %画出时域的图形
% %画出单个脉冲的形状
% figure(1);
% pulse=profile;                                     %pulse为不加载波的单个脉冲
% plot(t,abs(pulse).^2);
% figure(2);
% pulse_c=profile.*cos(2*pi*f0*t);                    %pulse_c为加载波的单个脉冲
% plot(t,abs(pulse_c).^2);
% %画出全部伪随机序列脉冲形状
% xt=1:length(sequence_c);
% xt=(ts/n)*xt;
% figure(3);
% plot(xt,abs(sequence).^2);                          %输出不加载波的脉冲序列
% figure(4);
% plot(xt,abs(sequence_c).^2);                        %输出加载波的脉冲序列

%**************************************************************************
%利用FFT进行时域到频域的变换
%注意：频域范围取决于采样点数n，曲线的精细程度取决于时域范围T。

dt=ts/n;                                              %dt为时域间隔
T=dt*length(sequence);                                %length(sequence)为脉冲在时域上的采样点数,也等于频域上的采样点数，即n*length(input)
df=1/T;                                               %T为时域范围，df为频域间隔
F=df*length(sequence);                                %F为频域范围，也等于1/dt
f_x=df*(-length(sequence)/2:length(sequence)/2-1);        

%f_x=linspace(-T/2*df,T/2*df,T)
%傅立叶变换，再fftshift移位
spectra_amp=fftshift(fft(sequence))/length(sequence); %傅立叶变换后要除以采样点数
spectra_pow=abs(spectra_amp).^2;
spectra_log=10*log10(abs(spectra_amp).^2+1e-6);       %对数坐标更接近实际显示
figure(5);
plot(f_x/1e9,spectra_pow);                            %纵坐标换成GHz,观察更方便