ofdm_xunhuanpu.m

％ 仿真主要完成OFDM系统建模,实现OFDM系统传输流程, ％ 并对OFDM信号的时域波形及功率普密度进行分析.进行循环谱相关切片的分析。

% 
% ofdm97.m
%
% 仿真主要完成OFDM系统建模,实现OFDM系统传输流程,
% 并对OFDM信号的时域波形及功率普密度进行分析.
clear;clc;
%********************** 参数设定部分 ***************************
bitrate = 64;         % 系统传输速率
Tu =0.8;              % 一个OFDM信号中有用信号持续时间[S],1/Tu为OFDM信号中子载波的间隔
Tp = 0.2;             % 一个OFDM信号中循环前缀持续时间[S]
Tg =0.2;              % 总保护间隔长度[S]
nf =1500;             % 一帧设定的OFDM符号
ml =4;                % 映射类别:BPSK=1,QPSK=2,16QAM=4,64QAM=6
over=4;               % 过采样倍数
NFFT = 16;            % FFT点数
A = 1;                % 形成OFDM信号的成形脉冲幅度[V]
fp = 25;              % 载频
% fs =80;               % 采样频率
% ts = 1/fs;            % 采样间隔
%----------------------主要参数计算--------------------------------------
nb = bitrate*(Tu+Tp);   % 一个符号所需传送的bit数:64比特
numcarrier = nb/ml;     % 子载波数
numbits = nb*nf;        % 一帧的bit数
%--------------------------------------------------------------------------
tc = Tu / NFFT;            % 一个切普持续时间1.0000e-007
ncp = floor(Tp/tc);        % Number of tones of the cyclic prefix:8
n = (-ncp+1:1:NFFT);   
N = length(n);             % Total number of tones per symbol:40
%--------------------------------------------------------------------------
if ncp>0
   tc = (Tu+Tp)/N ;                        % Tone duration
end

%************************** 发射部分 **************************************

%**************************************************************************
%                            数据产生
%**************************************************************************

ge_bits=randint(1,numbits); % 产生二进制数据流

%**************************************************************************
%                             星座点映射
%**************************************************************************
S =mapping(ge_bits,ml);

%**************************************************************************
%                            串并转换
%**************************************************************************

para_S = reshape(S,numcarrier,nf);  % 串并转换
%*************************************************************************
%                           IFFT变换
%*************************************************************************
%----------------------插零满足IFFT变换长度---------------------------------
%---------------------------第一种插零方式----------------------------------
para_S0 = zeros(NFFT,nf);
up1 = floor(numcarrier/2);
down1 = numcarrier-up1;
para_S0(1:up1,:) = para_S(1:up1,:);
para_S0(NFFT-down1+1:NFFT,:) = para_S(numcarrier-down1+1:numcarrier,:);
%---------------------------第二种插零方式----------------------------------
% para_S1 = zeros(2*numcarrier,nf);
% para_S0=zeros(2*numcarrier,nf);
% for a=1:nf
% para_S0(1,a)=real(para_S(1,a));
% para_S0(17,a)=imag(para_S(1,a));
% para_S0(2:16,a)=para_S(2:end,a);
% para_S0(32:-1:18,a)=conj(para_S(2:end,a));
% end
%------------------------四倍过采样-----------------------------------------
up2 = floor(NFFT/2);
down2 = NFFT-up2;
FS = over*NFFT;                          %FS=32
FS_cp =over*(ncp+NFFT);                  %FS_cp=40
ifftlen = nf*FS_cp;                      %ifftlen=20000
paras_cp=zeros(FS_cp,nf);
 for b = 1:nf                            %循环说明一个一个发送符号
    ch = para_S0(:,b); 
    C_zp = zeros(FS,1);                   %4倍过采样中间补0(zero padding)
    C_zp(1:up2)=ch(1:up2);
    C_zp(FS-down2+1:FS)=ch(NFFT-down2+1:NFFT);
%------------------------IFFT-----------------------------------------
    C = ifft(C_zp,FS);     % IFFT of the zero-padded input
    C = sqrt(FS)*C;

%**************************************************************************
%                            插入CP
%**************************************************************************
     if ncp>0 
        C_cp=zeros(length(C)+over*ncp,1);
        C_cp(1:(over*ncp))=C(over*NFFT-(over*ncp)+1:over*NFFT);  
        C_cp(over*ncp+1:length(C_cp))=C;
     else
        C_cp=C;
     end
   paras_cp(:,b)=C_cp(:);
  
 end

%**************************************************************************
%                           并串转换
%**************************************************************************

Stx_data2=reshape(paras_cp,ifftlen,1);
Stx_data=Stx_data2.';


%--------------------------------------------------------------------------
t3 = linspace (0,nf*(Tu+Tg),ifftlen);  %时域波形
%figure
%plot(t3,real(Stx_data))
% X = xlabel('Time [S]');
% Y = ylabel('Amplitude [V]');
% title('一帧基带OFDM信号同相分量波形')

df=FS_cp/(Tu+Tp);      %10000000Hz                  
[PSD2,frequency2] = cp_PSD(Stx_data,df);     %频域波形(功率普密度)
%figure
%plot(frequency2,10*log10(PSD2))
% X = xlabel('Frequency [Hz]');
% Y = ylabel('PSD [dB]');
% title('一帧基带OFDM符号的功率普密度')
%Stx_data=zeros(1,length(Stx_data));
%t=linspace(0,nf*(Tu+Tg),length(Stx_data));
%Stx_ofdm=Stx_data.*exp(j*2*pi*fp*t);
Stx_ofdm=Stx_data;
%**************************************************************************
%                           加性高斯白噪声信道
%**************************************************************************
 snrdb=30;
Ec=sum(abs(Stx_ofdm).^2)./length(Stx_ofdm);
ecn = 10.^(snrdb/10);

 noise_var = sqrt(Ec/(2*ecn));
 [Rx_data] = add_whitenoise2(Stx_ofdm,noise_var);


%*************************************************************************
%                           接收信号检测部分
%**************************************************************************
redata=Rx_data;
redata1=redata-mean(redata);   %信号去直流处理
power_data=mean(redata1.*conj(redata1));
redata2=redata1./sqrt(power_data);   %信号功率归一化处理
[PSD3,frequency3] = cp_PSD(redata2,df);     %频域波形(功率普密度)
figure
plot(frequency3,10*log10(PSD3))
X = xlabel('Frequency [Hz]');
Y = ylabel('PSD [dB]');
title('一帧基带OFDM符号的归一化功率普密度')
T1=nf*(Tu+Tg);           % 信号总的持续时间
Smaple=length(redata2);  % 截获数据长度
Ts=T1/Smaple;            % 求得信号的采样间隔
fs=1/Ts;                 % 信号的采样频率fs= 10000000Hz
deletf=1/Tu;
datalen=5000;            % 选取处理数据长度
a=[-4:0.1:4];            % 循环频率搜索范围
Ttao=80;
Ce=Ttao+100;
TTtao=Ttao/2;             %TTtao
%T0=Tu+Tg;                  
%F0=1/T0;
%a=a*F0;                                        %OFDM的循环频率应该是Fs的整数倍处,其他地方的Rx(a)应该等于零,即所得图形的峰值应该出现在0,0.25MHz,0.5MHz,0.75MHz,1MHze...处
T=[(1:datalen)*Ts];
ai=length(a);
f=[-40:0.1:40];            % 频率搜索范围
tao1=[-TTtao:1:TTtao].*Ts.*2; %[-15:1:15]./100.*2
R_alpha=zeros(Ttao+1,length(a));
S_alpha=zeros(Ttao+1,length(f));
for i=1:5
 for tao=-TTtao:1:TTtao                            
    R_tao=redata2(Ce-tao+(i-1)*5200:Ce-tao+datalen-1+(i-1)*5200).*conj(redata2(Ce+tao+(i-1)*5200:Ce+tao+datalen-1+(i-1)*5200));
    temp(tao+TTtao+1,:)=(1/datalen)*R_tao*exp(-j*2*pi*T'*a); %循环自相关函数
 end
% R_alpha=R_alpha+temp;
temp_zz=temp.';
for ii=1:ai
    S_temp(ii,:)=temp_zz(ii,:)*exp(-j*tao1'*f*2*pi)./fs;
end
%S_temp_zz=S_temp.';
S_alpha=S_alpha+S_temp;
end
 S_alpha2=S_alpha/5;
% %size(R_alpha)
% figure                                  %显示循环自相关函数
% plot(a,abs(R_alpha2(9,:)),'b');             %经过理论分析(IEEE文章中写的)应该在第10行即延迟等于3.2us处出现峰值,因此取了延迟等于3.2e-6s的切面
% xlabel('α [Hz]');
% ylabel('Rx(α)');
% title('信号循环自相关函数');
% grid;
% tao1=(-80:2:80)*Ts;
% figure
% plot(tao1,abs(R_alpha2(:,31)),'b');
% xlabel('延迟 [s]');
% ylabel('Rx(tao)');
% title('信号循环自相关函数');
% grid;
% 
% figure                            
% [x,y]=meshgrid((-4:0.1:4)*F0,(-80:2:80)*Ts);
% surf(x,y,abs(R_alpha2))
% colormap([1 1 1]);                             %将图像变为全黑
% xlabel('α [Hz]');
% ylabel('延迟 [s]');
% zlim([0 1])
% zlabel('Rx(α)')
% title('三维循环自相关函数');
% 
% ii=length(a);
% R_alpha_zz=R_alpha.';  %求转置,横行是不同的时延，竖行是不同的循环频率
% f=[-15:1:15];  % spectrum frequence scope1*121
% tao=[-TTtao:1:TTtao].*ts.*2; %[-15:1:15]./100.*2
% S_alpha=zeros(81,length(f));
% for i=1:ii
%     S_alpha(i,:)=R_alpha_zz(i,:)*exp(-j*tao'*f*2*pi)./fs;
% end
% S_alpha_zz=S_alpha.';
[AAA,F]=meshgrid(f,a);
figure,surf(AAA,F,abs(S_alpha2))
colormap([1 1 1]);   
%取f=0的切面
data=S_alpha2(:,401);
normal_data=data/max(data);
figure
plot(a,abs(normal_data))
xlabel('频率 [Hz]');
ylabel('S0(f)');
title('OFDM谱相关函数f=0的二维切面');
grid;

%取a=0的切面
data1=S_alpha2(41,:);
normal_data1=data1/max(data1);
figure
plot(f,abs(normal_data1))
xlabel('频率 [Hz]');
ylabel('S0(f)');
title('OFDM谱相关函数α=0的二维切面');
grid;