park.txt

OFDM同步中的PARK经典算法

% Simulation program to realize SAC algorithm
% (under five path rayleigh fading and awgn channel)
%

% clear;%clc;

%********************** preparation part ***************************

para=256;   % Number of parallel channel to transmit (points)
fftlen=256; % FFT length
noc=256;    % Number of carrier
nd=1;       % Number of information OFDM symbol for one loop
ml=2;       % Modulation level : QPSK
sr=5*10^6/fftlen;  % Symbol rate
br=sr.*ml;  % Bit rate per carrier
gilen=16;   % Length of guard interval (points)
fftlen2=fftlen+gilen;
v1=0.4;      % frequency offset normalized by subcarrier interval

% snr=3;
for snr=1:1:6,      % signal to noise ratio
ebn0=5*snr-10*log10(2);    % Eb/N0
snr=10^(ebn0/10);
%******************* Fading initialization ********************
% If you use fading function "sefade", you can initialize all of parameters.
% Otherwise you can comment out the following initialization.
% The detailed explanation of all of valiables are mentioned in Program 2-8.

% Time resolution

tstp=1/sr/(fftlen+gilen); 

% Arrival time for each multipath normalized by tstp
% If you would like to simulate under one path fading model, you have only to set 
% direct wave.

itau = [0 3 5 9 12];

% Mean power for each multipath normalized by direct wave.
% If you would like to simulate under one path fading model, you have only to set 
% direct wave.
dlvl = -20*log10(exp(-itau));

% Number of waves to generate fading for each multipath.
% In normal case, more than six waves are needed to generate Rayleigh fading
n0=[6 6 6 6 6];

% Initial Phase of delayed wave
% In this simulation four-path Rayleigh fading are considered.
%th1=[0 0 0 0 0];

% Number of fading counter to skip 
itnd0=nd*(fftlen+gilen)*10;

% Initial value of fading counter
% In this simulation one-path Rayleigh fading are considered.
% Therefore one fading counter are needed.
  
itnd1=1000:1000:5000;

% Number of directwave + Number of delayed wave
% In this simulation one-path Rayleigh fading are considered
now1=5;        

% Maximum Doppler frequency [Hz]
% You can insert your favorite value
fd=48.1;       

% You can decide two mode to simulate fading by changing the variable flat
% flat     : flat fading or not 
% (1->flat (only amplitude is fluctuated),0->nomal(phase and amplitude are fluctutated)
flat =1;

%************************** main loop part **************************

nloop=100;  % Number of simulation loops
for iii=1:nloop

%************************** transmitter *********************************

%************************** Data generation **************************** 

seldata_b=rand(1,para*nd*ml)>0.5;  %  rand : built in function
seldata_d=rand(1,para*nd*ml)>0.5;  %  rand : built in function

%****************** Serial to parallel conversion ***********************

paradata_b=reshape(seldata_b,para,nd*ml); %  reshape : built in function
paradata_d=reshape(seldata_d,para,nd*ml); %  reshape : built in function

%************************** QPSK modulation ***************************** 

[ich_b,qch_b]=qpskmod(paradata_b,para,nd,ml);
ich_b(2:2:fftlen)=0;

[ich_d,qch_d]=qpskmod(paradata_d,para,nd,ml);
ich_d=ich_d/sqrt(2);
qch_d=qch_d/sqrt(2);
%******************* IFFT ************************

x_b=ich_b;
y_b=sqrt(2)*ifft(x_b);      %  ifft : built in function
ich2_b=real(y_b);   %  real : built in function
qch2_b=imag(y_b);   %  imag : built in function

x_d=ich_d+qch_d.*i;
y_d=ifft(x_d);      %  ifft : built in function
ich2_d=real(y_d);   %  real : built in function
qch2_d=imag(y_d);   %  imag : built in function

%********* Gurad interval insertion **********
[ich3_b,qch3_b]= giins(ich2_b,qch2_b,fftlen,gilen,nd);
[ich3_d,qch3_d]= giins(ich2_d,qch2_d,fftlen,gilen,nd);

%********* Attenuation Calculation *********

spow_b=sum(ich3_b.^2+qch3_b.^2)/nd./para;  %  sum : built in function
attn_b=sqrt(0.5*spow_b*sr/br*10.^(-ebn0/10));

spow_d=sum(ich3_d.^2+qch3_d.^2)/nd./para;  %  sum : built in function
attn_d=sqrt(0.5*spow_d*sr/br*10.^(-ebn0/10));

%********************** Fading channel **********************
% [ifade_b,qfade_b]=sefade(ich3_b,qch3_b,itau,dlvl,th1,n0,itnd1,now1,length(ich3_b),tstp,fd,flat);
% [ifade_d,qfade_d]=sefade(ich3_d,qch3_d,itau,dlvl,th1,n0,itnd1,now1,length(ich3_b),tstp,fd,flat);
% engy1(iii)=sum(abs(ich3_b+i*qch3_b).^2); engy2(iii)=sum(abs(ifade_b+i*qfade_b).^2);
% ich3_b=ifade_b; qch3_b=qfade_b;
% ich3_d=ifade_d; qch3_d=qfade_d;
  
% ich4_b=ich3_b*ray_chan(iii); qch4_b=qch3_b*ray_chan(iii);
% ich4_d=ich3_d*ray_chan(iii); qch4_d=qch3_d*ray_chan(iii);

[ich4_b,qch4_b]=fade(ich3_b,qch3_b,length(ich3_b),tstp,fd,n0(1),itnd1(1),flat);%fade(idata,qdata,nsamp,tstp,fd,no,counter,flat)
[ich4_d,qch4_d]=fade(ich3_d,qch3_d,length(ich3_d),tstp,fd,n0(1),itnd1(1),flat);
% engy1(iii)=20*log10(sum(abs(ich3_b+i*qch3_b).^2)); 
% engy2(iii)=20*log10(sum(abs(ifade_b+i*qfade_b).^2));
% ich4_b=ifade_b; qch4_b=qfade_b;
% ich4_d=ifade_d; qch4_d=qfade_d;

% Updata fading counter
itnd1 = itnd1+ itnd0;

%********************** Add frequency offset************************
ich5_b=exp(j*2*pi*(-gilen:fftlen-1)*v1/fftlen).*ich4_b; 
qch5_b=exp(j*2*pi*(-gilen:fftlen-1)*v1/fftlen).*qch4_b;
ich5_d=exp(j*2*pi*(-gilen:fftlen-1)*v1/fftlen).*ich4_d; 
qch5_d=exp(j*2*pi*(-gilen:fftlen-1)*v1/fftlen).*qch4_d;


%***************** AWGN addition ********* 

[ich6_b,qch6_b]=comb(ich5_b,qch5_b,attn_b); % only white gaussian noise is considered
[ich6_d,qch6_d]=comb(ich5_d,qch5_d,attn_d);

%***************************  Receiver  *****************************
r_all_b=[ich6_d+i*qch6_d ich6_b+i*qch6_b ich6_d+i*qch6_d];

d=-gilen-fftlen/2:fftlen/2-1;
P_b=zeros(1,length(d)); R_b=zeros(1,length(d));
for m=1:fftlen/2-1,
    P_b=P_b+r_all_b(fftlen2+gilen+fftlen/2+1+d-m).*r_all_b(fftlen2+gilen+fftlen/2+1+d+m);
    R_b=R_b+abs(r_all_b(fftlen2+gilen+fftlen/2+1+d+m)).^2;
end;
P1_b=abs(P_b).^2; R1_b=abs(R_b).^2;
result_b=P1_b./R1_b;
index_b=find(result_b==max(result_b))-gilen-fftlen/2-1;
figure(1);
plot(result_b);
end;
mean_timing_b(snr+1)=mean(index_b)
var_timing_b(snr+1)=var(index_b)
error_b(snr+1)=length(find(index_b~=0))/iii
end;