%###### 基于DVB-T標準,COFDM調制系統的利用導頻信號進行符號粗同步 ######
%###### ?仿真條件：68個OFDM符號，64QAM調制，萊斯信道 ?######
% IFFT_bin_length: IFFT和FFT的點數
% carrier_count: 子載波個數
% bits_per_symbol: 每個載波含有的比特數
% symbols_per_carrier: 參與估計的總的OFDM符號數
% Ｘ：欲發送的二進制比特流

clear all;
close all;
rng(2);
IFFT_bin_length=2048;
carrier_count=1705;
bits_per_symbol=4;
symbols_per_carrier=68; % 總的OFDM符號數
LI=12; ? ? ? ? ? ? ? ? ?% 導頻之間的間隔
GI=2048/4; ? ? ? ? ? ? ?% 保護間隔長度
N_number=carrier_count*symbols_per_carrier*bits_per_symbol;
% 格雷碼特性設定
EbNo = 4.5:.5:7; linEbNo = 10.^(EbNo(:).*0.1);
ModulateIndex = 16; codeRate = 1/2; constlen = 7; k = log2(ModulateIndex); codegen = [171 133];
tblen = 32; ? ? % traceback length
trellis = poly2trellis(constlen, codegen);
dspec = distspec(trellis, 7);
expVitBER = bercoding(EbNo, 'conv', 'hard', codeRate, dspec);
%--------------------------------------------------------------------------
%矩陣初始化
XX=zeros(1,N_number);
dif_bit=zeros(1,N_number);
dif_bit1=zeros(1,N_number);
dif_bit2=zeros(1,N_number);
dif_bit3=zeros(1,N_number);
X=rand(1,N_number); ? ? ? ? ? %產生二進制隨機序列（非0即1）
X=randint(carrier_count*symbols_per_carrier,1,16);
%--------------------------------------------------------------------------
%% 16QAM調制 %%
numSymb = 1512*4*2;
hStr = RandStream('mt19937ar', 'Seed', 654321);
msg_orig = randi(hStr, [0 1], numSymb, 1);
msg_enc = convenc(msg_orig, trellis);
hMod = modem.qammod('M', ModulateIndex, 'PhaseOffset', 0, ...
? ? 'SymbolOrder', 'Gray', 'InputType', 'Bit');
msg_tx = modulate(hMod, msg_enc);

X1=msg_tx;

%--------------------------------------------------------------------------
%% 串并轉換 %%
X2=reshape(X1,1512,4).';

%--------------------------------------------------------------------------
%產生隨機導頻信號
Kmin=0;
Kmax=1704;
pilot(:,1)=Kmin; ? ? ? ? ? ? ? %使第一列是導頻
for l=0:symbols_per_carrier-1
? ? p=0;
? ? k=0;
? ? while k<=Kmax-12
? ? ? ?k=Kmin+3*mod(l,4)+12*p;
? ? ? ?pilot(l+1,p+2)=k; ? ? ? %分散導頻位置
? ? ? ?p=p+1;
? ? end
end
pilot(:,p+2)=Kmax; ? ? ? ? ? ? %使最后一列也是導頻
pilot=pilot+1;

% 16QAM的平均幅度
A_avg=(3*sqrt(2)*4+sqrt(2)*4+sqrt(10)*8)/16;
Burst=1*A_avg*4/3;
%--------------------------------------------------------------------------
%% 插入分散導頻 %%
g=bin2dec('100000000101'); ? ? %偽隨機二進制序列生成多項式:x11+x2+1
state=bin2dec('11111111111'); ?%偽隨機二進制序列生成寄存器初始狀態
N=2^11-1; ? ? ? ? ? ? ? ? ? ? ?%生成二進制序列長度
train=zeros(symbols_per_carrier,carrier_count);
train_sym=zeros(symbols_per_carrier,carrier_count);
for i=1:l+1
? ? m=round(rand(1,carrier_count)); ? %改成隨機數，不用上列PRBS生成法。簡化
? ? train(i,:)=m(1:carrier_count);
end

for i=1:l+1
? ? train_sym(i,pilot(i,:))=Burst*2.*(1/2-train(i,pilot(i,:))); %分散導頻值
end
signal=1:carrier_count;
X3(:,signal)=0;
for i=1:l+1
? ? X3(i,pilot(i,:))=train_sym(i,pilot(i,:)); ? ? ? ?%插入分散導頻
end
X3_SPCP=X3(1:4,:); ?%保留原始插入分散導頻
X3_SPCP(2:4,1)=0;X3_SPCP(2:4,1705)=0;
ScPilotX=X3(1:4,:); ? ? ?
%% 插入連續導頻 %%
CP=[0 48 54 87 141 156 192 201 255 279 282 333 432 450 483 525 531 618,...
? ? 636 714 759 765 780 804 873 888 918 939 942 969 984 1050 1101 1107,...
? ? 1110 1137 1140 1146 1206 1269 1323 1377 1491 1683 1704]+1; ? ? ? ? ? ? ? ? ? ?%45個連續導頻位置

for i=1:length(CP)
? ? train_sym(:,CP(i))=Burst*2.*(1/2-train(:,CP(i)));%連續導頻值
end
for i=1:length(CP)
? ? X3(:,CP(i))=train_sym(:,CP(i)); ? ? ? ? ? ? ? ? ?%插入連續導頻
end

%% 插入傳數參數信令 %%
TPS=[34 50 209 346 413 569 595 688 790 901 1073 1219 1262 1286 1469 1594 1687]+1; %17個傳數參數信令位置

for i=1:length(TPS)
? ? train_sym(:,TPS(i))=randint;
? ? train_sym(:,TPS(i))=A_avg*2*(1/2-train_sym(:,TPS(i))); %傳數參數信令值?
end
for i=1:length(TPS)
? ? X3(:,TPS(i))=train_sym(:,TPS(i)); ? ? ? ? ? ? ? ?%插入傳數參數信令 ? ??
end
%--------------------------------------------------------------------------
%% 插入數據
Data_index=zeros(4,1705);
X_data=X2;
X_out=X3(1:4,:);
for i=1:4
? ? m=1;
? ? for j1=1:1705
? ? ? ? if abs(X_out(i,j1))<0.1
? ? ? ? ? ? X_out(i,j1)=X_data(i,m);
? ? ? ? ? ? Data_index(i,j1)=3; ?% 記憶有效數據點處為3，導頻和TPS為0，為接收端提取數據用。
? ? ? ? ? ? m=m+1; ?%只有1512個數據，最后一個m多加了1.
? ? ? ? end
? ? end
? ? M(i)=m; ? ?
end
Data_index=[Data_index;Data_index;Data_index;Data_index];
? %組成68個OFDM符號，4個FODM符號為一循環
? ? X41=[X_out;X_out;X_out;X_out;X_out;X_out;X_out;X_out];
? ? X42=[X41;X41;X_out];
? ? CP_pilot=X42(1:6,:);
?% ------------------------------------------------------------------------
%% IFFT變換 %%
IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);
IFFT_modulation(:,signal)=X42;?
X4=ifft(IFFT_modulation,IFFT_bin_length,2); %2為行運算，1為列運算。
% 下面的方法與前一行運算等同。復數矩陣行列調換按如下方式進行，不能用reshape
%ifft_x=(IFFT_modulation).';
%X41=ifft(ifft_x,2048);
%X4=(X41).';

%--------------------------------------------------------------------------
%% 加循環前綴保護間隔 %%
X10=zeros(68,2560);
for j1=1:68
? ? X10(j1,1:GI)=X4(j1,2048-512+1:end);
? ? X10(j1,GI+1:end)=X4(j1,1:end);
end

%--------------------------------------------------------------------------
%% 并串轉換 %% ?reshape按列重排，X6是68x2560，必須先轉成為2560x68，再重排
X7=reshape(X10.',1,symbols_per_carrier*(IFFT_bin_length+GI));

% %--------------------------------------------------------------------------
% %% 加入基于DVB-T的瑞利多徑衰落信道(固定接受模式) %%
% 加衰減信道
??
chan = rayleighchan; ? ? ?% Static Rician channel
chan.PathDelays=[0 0.05 0.4 1.45 2.3 2.8]*1e-6; % 長時延
chan.AvgPathGaindB=[-2.8 0 -3.8 -0.1 -2.6 -1.3];?
? ? ? ? %chan.PathDelays=[0 5 14 35 54 75]*1e-6; % 長時延
? ? ? ? %chan.AvgPathGaindB=[0 -9 -22 -25 -27 -28];?
? ? ? ? %chan.PathDelays=[0 0.05 0.4 1.45 2.3 2.8]*1e-6; %短時延
? ? ? ? %chan.AvgPathGaindB=[-2.8 0 -3.8 -0.1 -2.6 -1.3];
chan.InputSamplePeriod=1e-7; ? ?
fadedSig = filter(chan, X7); ? ?% Apply channel effects
Tx_data = awgn(fadedSig,16,'measured'); ? ? ? ? ? % Add Gaussian noise. ??

Y7=[zeros(1,400) Tx_data(2560*3+1:end)]; % Tx_data(2560*2+1:end)] **********************************************************************信道參數

%% ----------最大似然法，確定粗定時偏差 --------------------------
r=Y7; ? ? ? ? ? %任意選取一段數據進行估計
N=IFFT_bin_length;
G=GI;
T=1:3*N+2*G;
for i=1:length(T) ? ? ? ? ? ?% FFT窗移動的距離
? ? data1=r(i:G+i-1); ? ? ? ?% 取GI長個數據存data1
? ? data2=r(N+i:N+G+i-1); ? ?% 取相距IFFT_bin_length的GI長個數據存data2
? ? s(i)=data1*data2'; ? ? ? % 求互相關值 ? ?
? ? cor(i)=abs(s(i));?
% ? ? cor(i)=abs(real(data1))*(abs(real(data2)))'+abs(imag(data1))*(abs(imag(data2)))';
end

% ---------------極大似然作圖----------------------------------------------
figure(1)
subplot(211);
plot(T,cor);
xlabel('載波數');
ylabel('相關值');
grid on;

[Max Max_i]=max(cor(1:N)); ? ? ? ? ? ? ? ? ?% Max_i對應一個符號的第一個數據
start_place=Max_i;
for i=1:length(T)
? ? deltaF(i)=1/(2*pi)*angle(conj(s(i)));
? ? r1=r(i:G+i-1);
? ? r2=r(i+N:G+i+N-1);
? ? e(i)=1/(2*pi)*atan(imag(r1*r2')/(real(r1*r2')));
end
% deltaF=1/(2*pi)*angle(conj(s(Max_i))); ? ? ?% 分數頻偏估計值
% r1=r(Max_i:G+Max_i-1);
% r2=r(Max_i+N:G+Max_i+N-1);
% e=1/(2*pi)*atan(imag(r1*r2')/(real(r1*r2')));
% disp(start_place);disp(e);
subplot(212); ?
plot(T,e);
xlabel('載波數');
ylabel('粗頻偏估計值');
%axis(1,length(T),-0.5,0.5);
grid on;
% Max_i=400;
r=r(Max_i-30:Max_i-30+12*2560);

%% -------------------------- FFT變換 ---------------------------
N=IFFT_bin_length;
Ng=GI;
Y8=r; ?% 加延時，使截取點落入循環前綴之內。
r=Y8;
r_ofdmin=reshape(r(1:12*2560),2560,12).';
r_nocp=r_ofdmin(:,GI+1:end);
r_Sym=fft(r_nocp,2048,2);
%figure(2)
scatterplot(r_Sym(1:1705));
title('存在定時誤差的FFT窗口取樣符號');

%% ------------------細定時估計------------------------------------
% 找出分散導頻位置
P=12;
Y_1=r_Sym(1,1:end);
Y_5=r_Sym(5,1:end);
% 找出導頻位置點，逐點相關。
for i=1:P*4;
? ? cc=0;
? ? for p=0:142
? ? ? ? cc=cc+Y_1(12*p+i)*conj(Y_5(12*p+i));
? ? end
? ? c(i)=abs(cc);
end
% 另一種找出導頻點的方法，更好且簡單，確定四種導頻方式其中的一種。
for k1=0:3
? ? dd=0;
? ? for p=1:142
? ? ? ? dd=dd+Y_1(12*p+1+3*k1)*conj(Y_5(12*p+1+3*k1));
? ? end
? ? d(k1+1)=abs(dd);
end

[cMax Max_ip]=max(d(1:4));
figure(3)
subplot(2,1,1)
plot(c);
subplot(2,1,2)
stem(d);
% 確定細定時偏差，假定定時偏差在（-85，+85）之間。
Y_SP=r_Sym(1,:); %取第1個OFDM符號。
sum=0;K1=-N/(2*pi*P);
SP_X=ScPilotX(Max_ip,:); ?% 取出模式Max_ip指定的分散導頻
SP_X=[SP_X, ?zeros(1,2*P)]; %補點，防止下標溢出。
SP_start=3*(Max_ip-1)+1;
for i=0:142
? ? sum=sum+conj(Y_SP(i*P+SP_start))*(SP_X(i*P+SP_start))...
? ? ? ? *conj(conj(Y_SP((i+1)*P+SP_start))*(SP_X((i+1)*P+SP_start)));
? ? %sum=sum+(Y_SP(i*P+SP_start))*conj(Y_SP((i+1)*P+SP_start));
end
err_avg=K1*angle(sum)/1;
Time_err=round(err_avg);
%Time_err=200; ?%保證取樣點正確，為星座圖便于理解的例外設定

% 修正定時誤差產生的旋轉量
? ? k=0:2048-1;
? ? X_modify=r_Sym(1:12,:);
? ? % Time_err=-16;
for i1=1:12
? ? X_modify1(i1,:)=X_modify(i1,:).*exp(j*2*pi*(Time_err)*k/N);
end
%figure(4)
scatterplot(X_modify1(Max_ip,:));
title('經定時誤差修正的符號');
% ---------------------------------------------------------

%% -------------------- 信道估計 ?--------------------------
% 簡化方式，只能對第一個OFDM符號做估計，且認為是分散導頻模式0.
% 只保留分散導頻和連續導頻 ?
r_chestimation=X_modify1(:,1:1705);
if (Max_ip==1)
? ? First_ip=1;
end
if (Max_ip==2)
? ? First_ip=4;
end
if (Max_ip==3)
? ? First_ip=3;
end
if (Max_ip==4)
? ? First_ip=2;
end
if (Max_ip>4)
? ? First_ip=1;
end
% First_ip=4;
r_chestimation=X_modify1(First_ip:First_ip+8-1,:);
X_modify2=r_chestimation;

r_chestimation(:,TPS)=0;
for m=1:8
? ? for k=1:1705
? ? ? ? if (abs(Data_index(m,k))>0.5)
? ? ? ? ? ? r_chestimation(m,k)=0;
? ? ? ? end
? ? end
end
r_chestimation_sum=zeros(1,1705); ?% *****************************************************************
X3_SPCP_sum=zeros(1,1705);
X3_SPCP12=[X3_SPCP;X3_SPCP;X3_SPCP];
CP_pilot12=[CP_pilot(1:4,:); CP_pilot(1:4,:); CP_pilot(1:4,:)];
X3_SPCP12(:,CP)=CP_pilot12(:,CP); ?% 發射端導頻，數據載波為0；
% 將四個符號的導頻合并到一個符號中
%for m=1:4
% ? ?r_chestimation_sum=r_chestimation_sum+r_chestimation(m,1:1705);
% ? ?X3_SPCP_sum=X3_SPCP_sum+X3_SPCP(m,1:1705);
%end

% 連續導頻用第一個符號的連續導頻置換。修改為Max_ip置換。

%CP_pilot
for i=1:8
? ? for k=1:1705;
? ? ? ? if(abs(X3_SPCP12(i,k))>0.001)
? ? ? ? ? ? Hp(i,k)=r_chestimation(i,k)./X3_SPCP12(i,k);
? ? ? ? end
? ? end
end
Hp_CP(:,CP)=Hp(:,CP); %保留連續導頻的信道估計待用

% 插值
% 列方向（時間）插值
?i=1;
for n=0:568-1
? ? while (i<=5) ? ?
? ? ? ? Hp(i+1,3*n+1)=(Hp(i,3*n+1)+3*Hp(i+4,3*n+1))/4;
? ? ? ? Hp(i+2,3*n+1)=(Hp(i,3*n+1)+Hp(i+4,3*n+1))/2;
? ? ? ? Hp(i+3,3*n+1)=(3*Hp(i,3*n+1)+Hp(i+4,3*n+1))/4;
? ? ? ? break ??
? ? end
? ? i=i+1;
? ? if i==5
? ? ? ? i=1;
? ? end
end
% 第一列連續導頻再替換回去。
Hp1=Hp;
Hp1(:,CP)=Hp_CP(:,CP);
% 行方向（頻率）插值，對第5個OFDM符號的信道估計
for n=0:568-1
? ? Hp1(5,3*n+2)=2*Hp1(5,3*n+1)/3+Hp1(5,3*n+4)/3;
? ? ?Hp1(5,3*n+3)=Hp1(5,3*n+1)/3+2*Hp1(5,3*n+4)/3;
end
% Hp1=Hp;?

%for k=0:568-1;
% ? ?Hp(1,3*k+2)=(Hp(1,3*k+1)*2/3)+(Hp(1,3*(k+1)+1))*1/3;
% ? ?Hp(1,3*k+3)=(Hp(1,3*k+1)*1/3)+(Hp(1,3*(k+1)+1))*2/3;
%end
% 信道估計修正接收的符號，第5個符號
% X_modify2(5,1:1705)=X_modify1(5,1:1705)./Hp1(5,1:1705);
X_modify3(5,1:1705)=X_modify2(5,1:1705)./Hp1(5,1:1705);
% figure(5)
?scatterplot(X_modify3(5,1:1705));
title('經信道估計修正的符號');
xlim([-5 5]); ylim([-5 5]);
% --------------統計超過幅度的符號點--------------------------------
recod=zeros(2,1705);sum1=0;
for k=1:1705
? ? if (abs(real(X_modify2(1,k)))>7)||(abs(imag(X_modify2(1,k)))>7)
? ? ? ? recod(1,k)=X_modify2(1,k);
? ? ? ? recod(2,k)=k;
? ? ? ? sum1=sum1+1;
? ? end
end

%% ------------------------------------------------------------- ? ?
? ?
% 從前1705個子載波中提取1512個數據符號，去除導頻和TPS。 ??
k2=1;
for k=1:1705
? ? if (abs(Data_index(5,k))>0.001)
? ? ? ? S_data(k2)=X_modify3(5,k);
? ? ? ? k2=k2+1;
? ? end
end
%figure(6)
scatterplot(S_data);
?title('抽取的數據符號');
?xlim([-5 5]); ylim([-5 5]);

%% -------------解調和viterbi解碼 --------------------------------
hDemod = modem.qamdemod('M', ModulateIndex, 'PhaseOffset', 0, ...
? ? 'SymbolOrder', 'Gray', 'OutputType', 'Bit');
% scatterplot(msg_rx_int);
msg_demod = demodulate(hDemod, S_data.');
msg_dec = vitdec(msg_demod, trellis, tblen, 'cont', 'hard');
[nChnlErrs BERChnl] = biterr(msg_enc(1:end/4), msg_demod);
[nCodErrs BERCoded] = biterr(msg_orig(1:end/4-tblen), msg_dec(1+tblen:end));
Time_err
Max_ip
Max_i
nChnlErrs
nCodErrs

?系統誤碼率仿真結果如下：

?

總結

以上是生活随笔為你收集整理的基于DVB-T标准,COFDM调制系统的利用导频信号进行符号粗同步的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。