实现北斗导航系统的boc信号的matlab调制

%function [signal_modu] = BOC_modulation(alpha,beta,N,t) % BOC_modulation(alpha,beta); % 因为扩频码位数原因，这里时间暂时只能设置为小于 1 ms， %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % BOC stands for Binary Offset Carrier % alpha and beta: alpha (respectively beta) is the ratio between the sub-carrier frequency % (resp. the code rate) and the reference frequency f0. Then : % f_sub = alpha*f0 and f_code = beta*f0 where f0 equals to 1.023 MHz. % N: N divide frequency of f0, the basic frequency, is regarded as data rate % fs: the sampling frequency ，the default is fs = 60 MHz； % t : duration time t of the signal, the default value is t = 1ms; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; close all; clc; alpha=1; beta=1; f0 = 1.023e6; t = 1e-3/beta; % 运算时间的确定，由于该程序用的是周期为1023的m序列，为了保持相关性，根据扩频码速率确定运算时间1e-3/beta N = 0.1; % sampling frequency % ------------------------ f_data = N*f0; % 数据消息速率 f_sub = alpha*f0; % 亚载波频率 f_code = beta*f0; % 扩频码速率 fs =200*f0 ; % 采样频率，16的话，图6横坐标不对64*alpha*f0；8*alpha*f0 n = 2*f_sub/f_code; delta_t = 0:1/fs:t-1/fs;%<1x16368 double> % ------------------------ % 数据信息加载模块 data_num = ceil(t*f_data); % 取大于浮点数的整数，以便后边根据载波长度进行截取 data = rand(1,data_num); % 随机数组1*103 for i = 1:data_num if data(i)<=0.5 data(i)=-1; else data(i) = 1; end end data_spread = repmat(data,fs/f_data,1); % 按照抽样频率进行扩展<160x103 double> data_temp = reshape(data_spread,1,[]);%<1x16480 double> data_avail = data_temp(1,1:length(delta_t)); % 截取与载波相同的位数进行运算<1x16368 double> % ------------------------ % 扩频码发生器模块（民用C/A码重复周期 1 ms，码见距离 1us，相当于300米） % 生产的周期为1023的M序列，为节省计算时间直接调用.mat文件，生成程序见.m文件 load codes codes; % （码速率为 1.023 MHz 试用，其他调制速率待解决 ） M_code = codes(1,1:f_code*t).*2-1; % M 扩频码<1x1023 double> M_spread = repmat(M_code,fs/f_code,1);%<16x1023 double> M_avail = reshape(M_spread,1,[]); % 经扩频符号扩展后的扩频码<1x16368 double> % ------------------------ % 方波发生器模块 %square_wave =square(2*pi*f_sub.*delta_t);%square(2*pi*f_sub.*delta_t) ;%sign(sin((2*pi*f_sub.*delta_t)));%<1x16368 double> BOC square_wave = sign(sin((2*pi*f_sub.*delta_t))); %square_wave = cos((2*pi*f_sub.*delta_t)); %square_wave = sign(exp((i*2*pi*f_sub.*delta_t)));%AltBOC %square_wave = sign(cos((2*pi*f_sub.*delta_t)))+i*sign(sin((2*pi*f_sub.*delta_t))); % ------------------------ % BOC信号 signal_modu = data_avail;%调制前 % ------------------------ figure(1); subplot(3,1,1);plot(delta_t,M_avail); v = axis;axis([0 5e-6 -1.5 1.5]); set(gca,'XTick',[0 5*1e-6 ]);% 对横坐标表示的方式进行定义 set(gca,'XTickLabel',{'0','5'}); xlabel('GOLD sequence /us','FontSize',10);%ylabel('A','Fontsize',12); subplot(3,1,2);plot(delta_t,square_wave); v = axis;axis([0 5e-6 -1.5 1.5]); set(gca,'XTick',[0 5*1e-6]); set(gca,'XTickLabel',{'0','5'}); xlabel('square wave /us','FontSize',10);%ylabel('A','Fontsize',12); %signal_modu = M_avail.*square_wave;%调制后data_avail.*M_avail.*square_wave; subplot(3,1,3);plot(delta_t,signal_modu); v = axis;axis([0 5e-6 -1.5 1.5]); xlabel('BOC signal /us','FontSize',10);%ylabel('A','Fontsize',12); set(gca,'XTick',[0 5*1e-6 10*1e-6 ]); set(gca,'XTickLabel',{'0','5'}); % figure(10); % subplot(3,1,1);plot(delta_t,fftshift(fft(abs(M_avail)))); % % v = axis;axis([0 20e-6 -1.5 1.5]); % % set(gca,'XTick',[0 5*1e-6 10*1e-6 15*1e-6 20*1e-6]);% 对横坐标表示的方式进行定义 % % set(gca,'XTickLabel',{'0','5','10','15','20'}); % xlabel('M sequence /us','FontSize',8);ylabel('A','Fontsize',12); % subplot(3,1,2);plot(delta_t,fftshift(fft(abs(square_wave)))); % % v = axis;axis([0 20e-6 -1.5 1.5]); % % set(gca,'XTick',[0 5*1e-6 10*1e-6 15*1e-6 20*1e-6]); % % set(gca,'XTickLabel',{'0','5','10','15','20'}); % xlabel('square wave /us','FontSize',8);ylabel('A','Fontsize',12); % subplot(3,1,3);plot(delta_t,fftshift(fft(abs(signal_modu)))); % % v = axis;axis([0 20e-6 -1.5 1.5]); % % xlabel('BOC signal /us','FontSize',8);ylabel('A','Fontsize',12); % % set(gca,'XTick',[0 5*1e-6 10*1e-6 15*1e-6 20*1e-6]); % % set(gca,'XTickLabel',{'0','5','10','15','20'}); % step1 = (-length(signal_modu)+1:length(signal_modu)-1)/length(signal_modu)*t; s_corr = xcorr(signal_modu); % figure(2); % plot(step1,s_corr/max(s_corr),'r'); %---------------------------------- % axis([-(alpha+5)*(1/(2*f_sub)) (alpha+5)*(1/(2*f_sub)) -1.2 1.2]); % ylabel('Normalized Correlation','FontSize',12); % xlabel('Delay（s）','FontSize',12); % title('归一化自相关函数','FontSize',12); % grid on; f=fs/2;%e740*f0 step2 = (-length(signal_modu)+1:length(signal_modu)-1)/length(signal_modu)*f; s_f = fftshift(fft(s_corr)); %BOC调制信号的自相关函数作傅里叶变换得到频谱密度 % figure(3); % plot(step2,10*log10(abs(s_f)));%/max(abs(s_f)) % xlabel('Frequency [Hz]','FontSize',12); % ylabel('PSD [dBW]','Fontsize',12); % title('BOC(10,5)调制前的实际信号的功率谱','FontSize',12); % %axis([-4e7 4e7 0 80]);%axis([-2e7 2e7 -70 10]); % grid on; figure(33); h1=spectrum.periodogram;%获得周期法对象的属性 psd(h1,signal_modu,'Fs',fs,'Centerdc',true);%：指示DC信号在twosided信号中间 title('PSD'); axis([-20 20 -120 -40]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % BOC 信号功率谱密度 % if rem(alpha/beta,2) == 0, % even = 'true'; % else % even = 'false'; % end delta_f =linspace(-80e6,80e6,9999) ;%delta_f = linspace(-15e6,15e6,10000); % difference with respect to L1 n = 2*f_sub/f_code; r = rem(n,2); % if r == 0, even = 'true'; end if r == 0, G = f_code* ( tan(pi*delta_f/(2*f_sub)).*sin(pi*delta_f/f_code)./(pi*delta_f+eps)).^2; else %G=8*f_code*(cos(pi*delta_f/f_code)).^2.*(1-cos(pi*delta_f/(f_code.*n)))./(pi^2*delta_f.^2.*(cos(pi*delta_f/(f_code.*n))).^2); G = f_code* ( tan(pi*delta_f/(2*f_sub)).*cos(pi*delta_f/f_code)./(pi*delta_f+eps)).^2; end % ------------------------ figure(4); %plot(delta_f,G); % Y 轴坐标为功率 h1 = plot(delta_f,10*log10(abs(G)/max(abs(G))), 'linewidth',1); % Y轴为对数坐标，使用的是对数刻度 xlabel('Frequency [MHz]','FontSize',12); ylabel('PSD [dBW]','Fontsize',12); legend('BOC(1,1)');%添加图例的标注 axis([-15e6 15e6 -30 10]) %axis([-15e6 15e6 -110 -60]); set(gca,'XTick',[-15*1e6 -10*1e6 -5*1e6 0 5*1e6 10*1e6 15*1e6]);% 对图中的横坐标标注按真实意义进行更改 set(gca,'XTickLabel',{'-15','-10','-5','0','5','10','15'}); % set(gca,'YTick',[ -30 -20 -10 -0 10]);% 对图中的横坐标标注按真实意义进行更改 % set(gca,'YTickLabel',{'-30','-20','-10','0','10'}); title('功率谱密度函数','FontSize',12); % % hold on; % % plot(delta_f.*pi/2, (sinc(pi*delta_f/(2*f_code))).^2./f_code); % G2=(sinc(pi*delta_f/(2*f_code))).^2./f_code; % h2 = plot(delta_f.*pi/2, 10*log10(abs(G2)),'linewidth', 1,'linestyle','-.');%.*pi/2 % xlabel('Frequency [MHz]','FontSize',12); % ylabel('PSD','Fontsize',12); % legend('GPS C/A'); % axis([-15e6 15e6 -100 -60]); % set(gca,'XTick',[-15*1e6 -10*1e6 -5*1e6 0 5*1e6 10*1e6 15*1e6]); % set(gca,'XTickLabel',{'-15','-10','-5','0','5','10','15'}); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 计算其自相关函数 % 1.通过功率谱密度与自相关函数的关系求 % (有点问题，相关函数正负峰数不对) % s_corr = fftshift(ifft((G)));%ifftshift(abs(ifft(G))); % step = linspace(-t,t,9999); % figure(2); % %plot(s_corr/max(s_corr),'r'); % plot(step,s_corr/max(s_corr),'r'); % axis([-(alpha+5)*(1/(2*f_sub)) (alpha+5)*(1/(2*f_sub)) -1.2 1.2]); fn=50e6;%*f0 signal_modu = data_avail.*M_avail.*square_wave.*cos(2*pi*fn*delta_t);%调制后 % a=0;b=sqrt(0.05);%均值为a，方差为b^2 % n=length(signal_modu); % noise=a+b*randn(1,n);%高斯白噪声 % *noise_f=fft(fftshift(fft(noise))); % figure; % plot(noise); % signal