BPSK调制解调 科斯塔斯环载波同步 加纳算法码元同步 Matlab .m文件

%%%%%%%%%%%%%%%%%%%%%%%%%% %程序功能：BPSK信号解调 2018.03.07 % %程序流程：码元同步 载波同步 判决输出 % %要 求：采样率为码元速率的整数倍 % %%%%%%%%%%%%%%%%%%%%%%%%%% clc; close all; clear all; Fs = 32e6; %信号采样率 Fb = 500000; %码元速率 Fc = 2e6; %实际载波频率 ts = 1/Fs; %时间分辨率 wfc = Fc+10000; %初始频率 Rate = Fs/Fb; %每个码元样点个数 num = 5e4; %样点个数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%生成BPSK信号%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%带通滤波器设计%%%% Fstop1 = Fc-Fb-0.2e6; % First Stopband Frequency Fpass1 = Fc-Fb; % First Passband Frequency Fpass2 = Fc+Fb; % Second Passband Frequency Fstop2 = Fc+Fb+0.2e6; % Second Stopband Frequency Dstop1 = 0.0001; % First Stopband Attenuation Dpass = 0.057501127785; % Passband Ripple Dstop2 = 0.0001; % Second Stopband Attenuation dens = 20; % Density Factor % Calculate the order from the parameters using FIRPMORD. [N, Fo, Ao, W] = firpmord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 1 ... 0], [Dstop1 Dpass Dstop2]); % Calculate the coefficients using the FIRPM function. bfir1 = firpm(N, Fo, Ao, W, {dens}); LData = ceil(num/Rate); %码元数量 Symbs = zeros(1,LData); for i=1:1:LData Symbs(i) = randi(2,1)-1; %随机产生符号 end constellation_map=[0 pi]; %星座图 Pskmodu = constellation_map(Symbs+1); angl = zeros(num,1); for i=1:1:num angl(i) = Pskmodu(floor((i-1)/Rate)+1); end SNR = -30; BPSK_Sig = zeros(1,num); for k=1:1:num BPSK_Sig(k) = 10000*(cos(2*pi*Fc*k*ts+angl(k))+sqrt(10^(SNR/10))*randn(1,1)); %%产生信号并加噪声 end Data = conv(BPSK_Sig,bfir1); if(num>length(Data)) num = length(Data); %判断是否超过数据长度 若超过则num等于数据长度 end %低通滤波器设计 Fpass = 0.5e6; % Passband Frequency Fstop = 2e6; % Stopband Frequency Dpass = 0.057501127785; % Passband Ripple Dstop = 0.0001; % Stopband Attenuation dens = 20; % Density Factor % Calculate the order from the parameters using FIRPMORD. [N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]); % Calculate the coefficients using the FIRPM function. bfir2 = firpm(N, Fo, Ao, W, {dens}); %滤波器系数 FirCoeNum = length(bfir2); %滤波器长度 phase_Save = zeros(floor(num/Rate),1); %保存的角度 CodeS_Save = zeros(floor(num/Rate),1); %输出码元 Freq_Out = zeros(num,1); %载频输出 PhaseDert_Out = zeros(num,1); %相位差输出 Data_DoFreq = zeros(1,num); %下变频后数据 Data_LoPass = zeros(1,num); %低通滤波后数据 phase = 0; %相位 ss = 0; dert_w = 0; dert_f = 0; temp = 0; XX = 0; YY = 0; %%%考虑到滤波器因素跳过前面一些点 for k=1:1:num Data_DoFreq(k) = Data(k)*exp(-1j*(phase)); %下变频 并将数据存储在Data_DoFreq中 if(k>=FirCoeNum) %因为滤波需要先去的数据 所以等累积的点数大于FircoeNum时再开始滤波 Data_LoPass(k) = Data_DoFreq(k+1-FirCoeNum:k)*bfir2'; %低通滤波 并存储数据 if(mod(k,Rate)==0) if(k>64*5) a = real(Data_LoPass(k-round(XX))); b = real(Data_LoPass(k-64-round(XX))); c = real(Data_LoPass(k-32-round(XX))); if(a~=0&&b~=0) iout = c*(a/abs(a)-b/abs(b)); else iout = 0; end XX =XX+0.0002*iout + 0.0001*temp; temp = iout; if(XX<0.5) XX = XX+64; end if(XX>64.5) XX = XX-64; end end bx = round(XX); phase_now = atan(imag(Data_LoPass(k-bx))/real(Data_LoPass(k-bx))); %求相位 ss = ss+1; %自加项 phase_Save(ss) = phase_now; %保存相位差 CodeS_Save(ss) = Data_LoPass(k-bx); %%科斯塔斯环频率跟踪 if(ss>5) %PID 调节 dert_f = phase_Save(ss) - phase_Save(ss-1); if(abs(dert_f)>2) %%若差过大则可能是噪声 dert_f = 0; end dert_w = 400 * phase_now+10*sum(phase_Save(ss-4:ss-1)); wfc = wfc + dert_w +2048*dert_f; end end end phase = 2*pi*ts*(wfc)+phase; %相位 Freq_Out(k) = wfc; end %%判决可能发生 0 1 颠倒 for i=1:1:length(CodeS_Save) if(CodeS_Save(i)>0) CodeS_Save(i) = 1; else CodeS_Save(i) = 0; end end subplot(2,2,1); plot(Freq_Out,'-'); title('载波频率'); subplot(2,2,2); plot(phase_Save,'-'); title('相位偏差'); subplot(2,2,3); plot(real(Data_LoPass),'.-'); title('码元抽取'); hold on; t(1:1:num) = 0; plot(t,'r'); hold on; for i=1:1:num/64-1 ipo = [i*64+64-floor(XX),i*64+64-floor(XX)+1j*real(Data_LoPass(i*64+64-floor(XX)))]; plot(ipo,'r'); hold on; end N = 200; subplot(4,2,6); plot(Symbs(N:N+50),'.-'); title('原始符号'); subplot(4,2,8); plot(CodeS_Save(N+4:N+50+4),'.-'); title('解调符号（可能反转）');