本次采用的的是频率同步方法为：costas环
MATLAB源代码：
%%%%%%%%%%%%%%%%%%%%%%%%%%
%程序功能：BPSK信号解调 2018.03.07 %
%程序流程：码元同步 载波同步 判决输出 %
%要 求：采样率为码元速率的整数倍 %
%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
close all;
clear all;
Fs = 40e6; %信号采样率
Fb = 100e3; %码元速率
Fc = 0.8e6; %实际载波频率
ts = 1/Fs; %时间分辨率
wfc = Fc+1e3; %初始频率
Rate = Fs/Fb; %每个码元样点个数
num = 10e4; %样点个数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%生成BPSK信号%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%带通滤波器设计%%%%
% Fstop1 = Fc-Fb-200e3; % First Stopband Frequency
% Fpass1 = Fc-Fb; % First Passband Frequency
% Fpass2 = Fc+Fb; % Second Passband Frequency
% Fstop2 = Fc+Fb+200e3; % 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 = -200;
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
%%%% 带通滤波器 fs 40mhz fstop1 0.2mhz fpass1 0.8mhz fpass2 1.2mhz fstop2
%%%% 1.8mhz 阶数 ：169
load('bandpass.mat')
Data = conv(BPSK_Sig,bandpass);
if(num>length(Data))
num = length(Data); %判断是否超过数据长度 若超过则num等于数据长度
end
%低通滤波器设计
% Fpass = 0.2e6; % Passband Frequency
% Fstop = 1e6; % 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}); %滤波器系数
% 低通滤波器 fs:40mhz fpass:0.1mhz fstop:1 wpass: 1 wstop:80
load('lowpass.mat')
FirCoeNum = length(lowpass); %滤波器长度
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)*lowpass'; %低通滤波 并存储数据
if(mod(k,Rate)==0)
if(k>Rate*5)
a = real(Data_LoPass(k-round(XX)));
b = real(Data_LoPass(k-64-round(XX)));
c = real(Data_LoPass(k-32-round(XX)));
a1(k) = a;
b1(k) = b;
c1(k) = c;
if(a~=0&&b~=0)
iout = c*(a/abs(a)-b/abs(b));
iout1(k) = c*(a/abs(a)-b/abs(b));
else
iout = 0;
end
XX =XX+0.0002*iout + 0.0001*temp;
xx2(k) = XX;
temp = iout;
if(XX<0.5)
XX = XX+Rate;
end
if(XX>Rate+0.5)
XX = XX-Rate;
end
end
bx = round(XX);
xx1(k) = XX ;
phase_now = atan(imag(Data_LoPass(k-bx))/real(Data_LoPass(k-bx))); %求相位
phase_now1(k) = phase_now;
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; %相位修正
if phase >= 2*pi
phase = phase - 2*pi;
end
phase1(k) = phase;
Freq_Out(k) = wfc;
end
% figure()
% plot(real(Data_LoPass))
figure()
plot(Freq_Out)
% figure()
% plot(phase1)
% figure()
% plot(abs(fft(cos(phase1))))
%%判决可能发生 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
% figure()
% stem(CodeS_Save)
% hold on
% stem(Symbs,'r*')
% su
bp
lot(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('解调符号（可能反转）');
参考的别人的，然后就是基于源代码的simulink实现：
希望大家一起讨论，因为有些地方我还不是很理解，请多多交流。

DD

这个写错了，今天开会 被指正了，这个不是科斯塔斯环哈！