%********************Schmidl&Cox algorithm*******************
close all;
clear all;
clc;
% global Y1 Y %将变量定义为全局变量
%参数定义
N=32; %训练序列1长度
Ng=14; %循环前缀的长度 (保护间隔的长度)
Ns=Ng+N; %包括循环前缀的符号长度
Nuw=64; %UW 序列长度
Ndata=960; %数据负载长度
%************利用查表法生成复随机序列**********************
QAMTable=[7+7i,-7+7i,-7-7i,7-7i];
buf=QAMTable(randi([0 3],N/2,1)+1);
%********奇数子载波的位置插入零、偶数载波为复随机序列*******
x=zeros(N,1);
x(1:2:end)=buf(1:1:end);
%**************利用IFFT变换生成Schmidl训练符号***************
sch = ifft(x); %[A;A]的形式 “训练符号1”
sch = sch / max(sch);
%*****************添加一个后缀符号*************
src = QAMTable(randi([0 3],N,1)+1).';
%randint(N,1,4)+1从[1,4]中随机抽取数,组成N*1的矩阵
sym = ifft(src); %后缀序列
sig =[sch sym]; %N_width*2矩阵 zeros(N*width,1)
%**********************添加循环前缀*************************
tx = [sig(N - Ng +1:N,:);sig]; % 将sig后Ng行元素添加到前头
t = reshape(tx,1,size(tx,1)*size(tx,2));
%**********************添加UW序列*************************
uw=exp(1i*pi*(1:Nuw).^2./Nuw);%UW序列
%**********************数据负载及QPSK映射*************************
data = randi([0 3],Ndata,1); %数据负载
hMod = comm.QPSKModulator('PhaseOffset',pi/4);
data = step(hMod,data);
data = reshape(data,1,Ndata);
%**********************组帧*************************
frame=[t uw data uw]; %发送序列 1*数据数 15648
frame=reshape(frame,1,size(frame,1)*size(frame,2));
%**********************上采样*************************
L = 8; %采样倍数
frame_up =upsample(frame,L); % 1*(数据数*L)
%**********************成型滤波*************************
beta=0.3; %滚降系数,0~1间,决定频宽和陡峭程度,一般大于0.2
span=L*2; %span: 表示截断的符号范围。
sps=1; %sps: 单个符号范围的采样个数。
filiter1 = rcosdesign(beta,span,sps,'sqrt'); %根升余弦滤波器
s=upfirdn(frame_up, filiter1); % 成型滤波 发送数+span
send = reshape(s,1,size(s,1)*size(s,2));
send = [zeros(1,Ns*L),send]; %*width
%*********************信道*************************
%理想信道
receive=send;
%AWGN信道
noise=[-15,0,15];
receive_a=zeros(length(noise),length(send));
for i=1:1:length(noise)
receive_a(i,:)=awgn(send,noise(i),'measured');
end
%**********************匹配滤波*************************
receive=upfirdn(receive, filiter1); % 成型滤波 长度:发送数+span
receive1=zeros(length(noise),length(receive_a)+L*2);
for i=1:1:length(noise)
receive1(i,:)=upfirdn(receive_a(i,:), filiter1);
end
%*********************下采样*************************
%理想信道
r=downsample(receive,L); %下采样 2为偏移位数
%AWGN信道
r1=zeros(length(noise),size(receive1,2)/L); %/width
for i=1:1:length(noise)
r1(i,:)=downsample(receive1(i,:),L); %下采样
end
%**********************SC算法度量函数*************************
%理想信道
P=zeros(1,length(r));
R=zeros(1,length(r));
for d = 1:1:2*Ns
P(d) = sum(conj(r(d:d+N/2-1)).*r(d+N/2:d+N-1));
R(d) = sum(power(abs(r(d+N/2:d+N-1)),2));
end
M=power(abs(P),2)./power(abs(R),2);
[z,ro,w]=find(M>=1.0000); % M>=1的非零元素行、列的索引值、非零元素
a=max(ro);
b=min(ro);
%AWGN信道
M1=zeros(length(noise),length(r1));
P1=zeros(length(noise),length(r1));
R1=zeros(length(noise),length(r1));
for m=1:length(noise)
for d=(N/2+1):1:2*Ns
P1(m,d-N/2) = sum(conj(r1(m,d:d+N/2-1)).*r1(m,d+N/2:d+N-1));
R1(m,d-N/2) = sum(power(abs(r1(m,d+N/2:d+N-1)),2));
end
M1(m,:)=power(abs(P1(m,:)),2)./power(abs(R1(m,:)),2);
[z,g,w]=find(M1(m,:)>=max(M1(m,:)));
end
%**********************绘图******************************
% 理想信道能量差
figure(1);
d=1:1:2*Ns;
plot(d,M(d));
grid on;
axis([0,2*Ns,0,1.2]);
title('理想信道 SC算法');
xlabel('定点d');
ylabel('M(d)')
%AWGN信道
figure(2)
subplot(3,1,1)
d=1:1:2*Ns;
plot(d,M1(1,d));
grid on;
axis([0,2*Ns,0,1.2]);
title('AWGN信道 SC算法 SNR=-15dB');
xlabel('定点d');
ylabel('M1(d)');
subplot(3,1,2)
d=1:1:2*Ns;
plot(d,M1(2,d));
grid on;
axis([0,2*Ns,0,1.2]);
title('AWGN信道 SC算法 SNR=0dB');
xlabel('定点d');
ylabel('M1(d)');
subplot(3,1,3)
d=1:1:2*Ns;
plot(d,M1(3,d));
grid on;
axis([0,2*Ns,0,1.2]);
title('AWGN信道 SC算法 SNR=15dB');
xlabel('定点d');
ylabel('M1(d)');