function sigf=sigfre()
clear
pack
%**************************************************************************
%input7為7位，特定碼序列
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];
%**************************************************************************
%prbs為127位，7級偽隨機碼序列
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;%1在128位偽隨機碼中占的百分比
b=b/127;%0在128位偽隨機碼中占的百分比
%**************************************************************************
%脈沖初始情況的設定，
fs=2e10;                                    %輸入電信號時鐘，其輸出光脈沖時鐘為電信號兩倍
ts=1/(2*fs);                                %脈沖時間ts=25ps
p0=1;                                       %設定輸入光功率
f0=1.930897e14;                             %光波頻率
n=8;                                     %單個脈沖的取樣點
%**************************************************************************
%pulse為單個脈沖的包絡函數
%t為產生單個脈沖的離散點         
fwhm=ts/1.6;                                %單個脈沖的半高全寬的取值比較關鍵！
t=linspace(-0.0*ts,1.0*ts,n);               %單個脈沖的離散點的取樣點 
% pulse=p0*(sin((pi/2)*exp(-0.6931*(2*(t-0.5*ts)/fwhm).^4)));                 
pulse=p0*ones(1,n); 
%**************************************************************************
%完成整個序列取樣點的采集(加載數據信號到相位中)
%input7為特定序列，換成prbs7則為偽隨機序列，
%注意：如換為input7，則以下有兩處需要替換，并且上面prbs7的產生要注釋掉。
m=length(prbs7);                                %m為序列中脈沖個數
k=m*n;                                          %k為全部取樣點的個數，即序列中脈沖個數*單個脈沖的取樣點
sequence=zeros(1,k);                            %產生取樣點矩陣，sequence為不加載波的序列                              
for ii=1:m;
   if prbs7(ii)==1;                             %input7為特定序列，換成prbs7則為偽隨機序列
       for jj=1:n;
      sequence(((ii-1)*n+jj))=pulse(jj);      %數據為1則相位相同
       end
   else
       for jj=1:n;
      sequence(((ii-1)*n+jj))=-pulse(jj);     %數據為0則相位相反
       end
   end
end
%**************************************************************************
%加上載波信號
wt=linspace(0,(ts*m),k);                        %整個時域里取總的k個采樣點
sequence_c=sequence.*exp(i*2*pi*f0*wt);         %其中wt是一個小采樣點的時間長度,sequence_c為加載波的序列，exp(i*2*pi*f0*wt)為光載波
pulse_c=pulse.*exp(i*2*pi*f0*t);                    %pulse_c為加載波的單個脈沖
%**************************************************************************
%畫出時域的圖形
%畫出單個脈沖的形狀
% figure(1);
% plot(t/1e-12,abs(pulse).^2);                        %pulse為不加載波的單個脈沖
% ylabel('Normalized Power (au)');xlabel('Time (ps)');
% title('NRZ-DPSK pulse');
% figure(2);
% plot(t/1e-12,abs(pulse_c).^2);
% ylabel('Normalized Power (au)');xlabel('Time (ps)');
% title('NRZ-DPSK pulse');
% %畫出全部偽隨機序列脈沖形狀
% figure(3);
% plot(wt/1e-12,abs(sequence).^2);                          %輸出不加載波的脈沖序列
% set(gca,'Xtick',[0:ts/1e-12:m*ts/1e-12]);
% ylabel('Normalized Power (au)');xlabel('Time (ps)');
% title('NRZ-DPSK pulse');
% figure(4);
% plot(wt/1e-12,abs(sequence_c).^2);                        %輸出加載波的脈沖序列
% set(gca,'Xtick',[0:ts/1e-12:m*ts/1e-12]);
% ylabel('Normalized Power (au)');xlabel('Time (ps)');
% title('NRZ-DPSK pulse');

%**************************************************************************
%利用FFT進行時域到頻域的變換
%注意：頻域范圍取決于采樣點數n，n越小，頻域橫坐標越小
%注意：曲線的精細程度取決于時域范圍T，T=dt*n*m=ts*m，其中ts為定值，則脈沖個數m越大，越精細
dt=ts/n;                                              %dt為一個小采樣點的時域間隔
T=dt*k;                                               %T是總的時域長度,length(sequence)為脈沖在時域上的采樣點數,也等于頻域上的采樣點數，即n*length(input)
df=1/T;                                               %T為時域范圍，df為頻域間隔
F=df*k;                                               %F為總的頻域范圍，也等于1/dt
f_x=df*(-k/2:k/2-1);                                  %頻率在頻段上的取樣點      
%傅立葉變換，再fftshift移位 
sigf=fftshift(fft(sequence_c));
num1=length(y);
spectra_amp=fftshift(fft(sequence_c))/k;              %傅立葉變換后要除以采樣點數
spectra_pow=abs(spectra_amp).^2;
spectra_log=10*log10(abs(spectra_amp).^2+1e-10);      %對數坐標更接近實際顯示
%**************************************************************************
%畫出頻譜圖
figure(3);
plot(f_x/1e9,spectra_log);                            %縱坐標換成GHz,觀察更方便
axis([-120,120,-100,0]);
set(gca,'Xtick',[-120,-80,-40,0,40,80,120]);
ylabel('Power (dB)');xlabel('Frequence (GHz)');
title('NRZ-DPSK Spectrum');
%找到對應于載波頻率的數組位置
x=max(spectra_pow);
leg=length(spectra_pow);
spe=spectra_pow;
b=1;
for ii=1:leg
    if spe(ii)>=x
        x=spe(ii);
        b=ii;
    else 
        b=b;
    end
end
b;
x;
wl0=(3.0e+8)/(1.55e-6);
length(f_x);
f_x;
wl0=(3.0e+8)/f0;