基于Matlab的单频激光倍频的仿真，使用PRBS信号调制

%% 载波频率 % 为简化计算，我们产生了PRBS7的信号，其占12700个点即127位每位100个点，其采样频率为10e6 Hz。 % 而具体的实验中，我们使用的是1064nm的种子，其频率约为3*10e14 Hz。而实际信号的频率为800Mbps和1Gbps，即8*10e8和10e9 % 则实际情况下，一个PRBS信号频率为800/127=6.29M Hz，时长为1.59*e-7 s。 % 所以如果使用一个完整的PRBS7的周期作为计算，则光的频率对应为3*10e14 Hz /6.29M =4.77e7 角频率为2 \pi * 4.77e7 % =3e8。 %% 信号频率 fs=10^6;%采样频率 T=1/fs;%采样点数 L=12700;%信号长度 t=(0:L-1)*T;%时间序列 f=fs*(0:(L-1))/L; binwidth = fs/length(t); f_hz=[0:binwidth:fs-binwidth];%将点数转化为频率 %% 定义载波 w=3e8;%载波角频率 Z1=exp(1i*w*t);%载波 plot(t(1:500),Z1(1:500)) plot(f_hz,abs(fftshift(fft(Z1)))); %% %% 读取PRBS信号 load("Z.mat");%读取之前处理好的PRBS信号 plot(t,Z) % plot(f_hz,abs(fftshift(fft(Z)))); plot(f_hz,abs(fft(Z))); title('Amplitude Spectrum of PRBS7') xlabel('f (Hz)') ylabel('|P1(f)|') PRBS_F=fft(Z); plot(f_hz,PRBS_F) a=0.13;%用来调制PRBS信号的高频分量。 PRBS_F(a*length(PRBS_F):(1-a)*length(PRBS_F))=0; %将频率中的高频部分变为0 plot(f_hz,PRBS_F); %% 画出傅里叶谱的一半 P2 = PRBS_F/L; P1 = P2(1:L/2+1); P1(2:end-1) = 2*P1(2:end-1); f = fs*(0:(L/2))/L; plot(f,P1) title('Single-Sided Amplitude Spectrum of X(t)') xlabel('f (Hz)') ylabel('|P1(f)|') PRBS_lb=ifft(PRBS_F); plot(t,PRBS_lb); title('PRBS7') xlabel('t') ylabel('Signal strength') %% 调制 % y=exp(1i.*(2.*pi.*f_c.*t+beta.*pm)) beta=pi;%调制深度 Z2=exp(1i*w*t+1i*beta*Z)+randn(size(t))/10; plot(t,Z2); fft_z2=fft(Z2); plot(f_hz,abs(fftshift(fft_z2))/(length(Z2)/2)); title('PRBS7-FFT') xlabel('f(Hz)') ylabel('Signal strength') Z_RC=exp(1i*w*t+1i*beta*PRBS_lb)+randn(size(t))/10; % plot(t,Z_RC); fft_z2_RC=fft(Z_RC); plot(f_hz,abs(fftshift(fft_z2_RC))/(length(fft_z2_RC)/2)); title('PRBS7-FFT') xlabel('f(Hz)') ylabel('Signal strength') FFTC=abs(fftshift(fft_z2_RC))/(length(fft_z2_RC)/2); %% 将图像平移到中心，方便看图 C=max(FFTC,[],[1,2]);%获得极大值点 B=find(FFTC==C);%寻找极大值点对应的横坐标 FFTC1=FFTC(1:B*2-1); tll=(0:B*2-1-1)*T; f_hz_tll=[0:binwidth:tll-binwidth] plot(tll,FFTC1); plot(tll,20*log10(FFTC1)); title('PRBS7-FFT') xlabel('f') ylabel('Signal strength(dB)') %% 倍频 Z22=exp(2*(1i*w*t+1i*beta*Z))+randn(size(t))/10; plot(t,Z22); fft_z22=fft(Z22); plot(t,abs(fftshift(fft_z22))/(length(Z22))); Z12=exp(2*(1i*w*t))+randn(size(t))/10; plot(t,Z12); fft_z12=fft(Z12); plot(t,abs(fftshift(fft_z12))/(length(Z12)/2)); Z_RC2=exp(2*(1i*w*t+1i*beta*PRBS_lb))+randn(size(t))/10; plot(t,Z_RC2); fft_z2_RC2=fft(Z_RC2); plot(f_hz,abs(fftshift(fft_z2_RC2))/(length(fft_z2_RC2)/2)); title('SHG-PRBS7-FFT') xlabel('f(Hz)') ylabel('Signal strength') % FFTC2=abs(fftshift(fft_z2_RC2))/(length(fft_z2_RC2)/2); % C2=max(FFTC2,[],[1,2]);%获得极大值点 % B2=find(FFTC2==C2);%寻找极大值点对应的横坐标 % FFTC12=FFTC2(1:B2*2-1); % tll2=(0:B2*2-1-1)*T; % plot(tll2,FFTC12); % plot(tll2,20*log10(FFTC12)); %% 添加sin % sin函数 N1=64; si=0.5.*sin(2*pi*N1*t/(T*L))+0.5; si2=0.5.*sin(2*pi*N1*t/(T*L))+0.5; plot(t,si); % beta=pi/2; % Z1_sin=exp(1i*w*t+2*beta*1i*si); % plot(t,Z1_sin); % plot(f_hz,abs(fftshift(fft(Z1_sin)))); % Z_RC_sin=exp(1i*w*t+1i*beta*PRBS_lb+1i*2*beta*si);%+randn(size(t))/10; Z_RC_sin=exp(1i*w*t+1i*beta*si); plot(t,Z_RC_sin); fft_z2_sin=fft(Z_RC_sin); plot(f_hz,abs(fftshift(fft_z2_sin))/(length(Z_RC_sin)/2)); title('Sin-FFT') xlabel('f(Hz)') ylabel('Signal strength') %% 方波 N1=64; pm=0.5.*square(2*pi*N1*t/(T*L))+0.5; plot(t,pm); % Z_RC_pm=exp(1i*w*t+1i*beta*PRBS_lb+1i*beta*pm); Z_RC_pm=exp(1i*w*t+1i*beta*pm)+randn(size(t))/10; plot(t,Z_RC_pm); fft_z2_pm=fft(Z_RC_pm); plot(f_hz,abs(fftshift(fft_z2_pm))/(length(Z_RC_pm)/2)); title('Square-FFT') xlabel('f(Hz)') ylabel('Signal strength') % plot(pm);%产生了N1个方波 用来对比PRBS信号 %% 三角波 san= 0.5*sawtooth(2*pi*N1*t/(T*L))+0.5; plot(t,san); san_fft=fft(san); plot(t,abs(san_fft)); Z_san=exp(1i*w*t+1i*beta*san); plot(t,Z_san); Z_san_fft=fft(Z_san); plot(t,abs(Z_san_fft)/(length(Z_san)/2)); title('Sawtooth-FFT') xlabel('f(Hz)') ylabel('Signal strength')