%2次平方倍频法,估计qpsk信号载频%
clear
clc
%固定参数设定:信息码长度,载频,采样率,mc次数%
dn=5; %信息码长度
fs=60*10^6; %载波频率
fn=600*10^6; %生成信号的采样率(很大,近似模拟信号)
mctime=100; %蒙特卡罗仿真次数
fcenter=fs; %接收机滤波器中心频率
Bi=20*10^6; %接收机滤波器带宽
fsr=50*10^6; %信号采样频率,带通采样
%*********************************%
[data]=qpsk1023(dn,fs,fn); %信源信号,data的长度=1023*dn*fn/10*10^6;
%*********************************%
db=-10;
s=awgn(data,db,'measured');
I=xcorr(s);
plot(I(170000:300000))
%设定二次滤波器参数%
fc2=10*10^6;
Bi2=10*10^6;
wn=[2*(fc2-(Bi2/2))/fsr 2*(fc2+(Bi2/2))/fsr]; %注意要乘2。因为1代表fn/2。
b=fir1(64,wn,'DC-0');
%*********************************%
clear wn fc2 Bi2
% %平方倍频法估计qpsk调制信号载频,功率谱估计法采用平均周期图法%
% %****************************平均周期图法功率谱估计的检测概率计算bpsk1023***************************%
% relativeerror=zeros(1,14);
% correcttime=zeros(1,14);
% for db=-20:-7
% ftotal=0;
% dbsign=db+21; %gate1的下标为db+29是信噪比和下标的统一
% for mc=1:mctime
% % %*******信号通过信道,接收机接收,带通采样*********%
% s=awgn(data,db,'measured');
% x=receiver(data,fcenter,Bi,fn,fsr); %x为接收待处理的信号
% % x=data;
% % %****************%
%
% x=x.^4; %4倍频信号
% % x=filter(b,1,x); %滤波
%
% fx=fft(x,2^15); %计算信号功率谱
% px=((abs(fx)).^2)/(length(x));
%
%
% plotglp(px,1,fsr)
% pfinal=px(1:round(length(px)/2)); %搜索频谱最大处,估计载频
% indexmax=find(pfinal==max(pfinal));
% festimat=indexmax*fsr/length(px);
% f=(250*10^6-festimat)/4;
% reerror=abs((f-fs)/(fs));
% relativeerror(dbsign)=relativeerror(dbsign)+reerror;
%
% if reerror<0.0001
% correcttime(dbsign)=correcttime(dbsign)+1;
% end
%
% clear s x fx px pfinal indexmax f reerror
% mc
% end
%
% relativeerror(dbsign)=relativeerror(dbsign)/mctime;
% clear mc dbsign
% %显示进度%
% for i=1:10
% db
% end
% clear i
% %********%
% end
% festimate1_re_10_q1023=relativeerror;
% festimate1_ct_10_q1023=correcttime;
% save festimate1_re_10_q1023 festimate1_re_10_q1023
% save festimate1_ct_10_q1023 festimate1_ct_10_q1023
% clear festimate1_re_10_q1023 relativeerror festimate1_ct_10_q1023 correcttime