ofdm_cyclic_spectrum.m

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

% 
% ofdm97.m
%
% 仿真主要完成OFDM系统建模,实现OFDM系统传输流程,
% 并对OFDM信号的时域波形及功率普密度进行分析.
clear;clc;
%********************** 参数设定部分 ***************************
bitrate = 64;         % 系统传输速率
Tu =1;              % 一个OFDM信号中有用信号持续时间[S],1/Tu为OFDM信号中子载波的间隔
Tp = 0;             % 一个OFDM信号中循环前缀持续时间[S]
Tg =0;              % 总保护间隔长度[S]
nf =330;             % 一帧设定的OFDM符号
ml =4;                % 映射类别:BPSK=1,QPSK=2,16QAM=4,64QAM=6
over=16;               % 过采样倍数
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
 
% %**************************************************************************
% %                           信号加窗与并串转换
% %**************************************************************************
% beta=0.1; %滚降系数
% T_cp=Tu+Tp;  %一个OFDM符号持续时间
% df=FS_cp/(T_cp); %160000000Hz：信号的采样速率
% postfix=floor(beta*(T_cp)*df); %循环后缀中的采样点数，ncp为循环前缀的采样点数
% nflen=FS_cp+postfix;  % 加窗后一个OFDM符号的采样点数
% %ifftlen = nf*FS_cp;                      %ifftlen=160*nf
% temp=zeros(nflen,nf);
% seridata=zeros(1,(ifftlen+postfix));
% %temp1=zeros(1,(ifftlen+postfix));
% ws=zeros(nflen,nf);
% %tt1 = linspace (0,beta*T_cp,postfix);
% %w1 = 0.5+0.5*cos(pi+tt1*pi./(beta*T_cp));
% %tt2=linspace(beta*T_cp,T_cp,(nflen-postfix));
% %w2=1;
% %tt3=linspace(T_cp,(T_cp+beta*T_cp),postfix);
% %w3 = 0.5+0.5*cos((tt3-T_cp)*pi./(beta*T_cp));
% for c=1:nf
%     temp(1:FS_cp,c)=paras_cp(:,c);
%     temp((FS_cp+1):nflen,c)=paras_cp((over*ncp+1):(over*ncp+postfix),c);
%     tt1 = linspace((c-1)*T_cp,(beta+c-1)*T_cp,postfix);
%     w1 = 0.5+0.5*cos(pi+tt1*pi./(beta*T_cp));
%     tt2=linspace((beta+c-1)*T_cp,c*T_cp,(nflen-postfix));
%     w2=1;
%     tt3=linspace(c*T_cp,(beta+c)*T_cp,postfix);
%     w3 = 0.5+0.5*cos((tt3-T_cp)*pi./(beta*T_cp));
%     ws(1:postfix,c)=temp(1:postfix,c).*w1.';
%     ws(postfix+1:nflen-postfix,c)=temp(postfix+1:nflen-postfix,c).*w2.';
%     ws(nflen-postfix+1:nflen,c)=temp(nflen-postfix+1:nflen,c).*w3.';   %ws=zeros(nflen,nf),加窗后OFDM符号的采样点数
% end
% for d=1:nf
%    temp1=zeros(1,(ifftlen+postfix));
%    temp1((d-1)*FS_cp+1:(d-1)*FS_cp+nflen)=ws(:,d);
%    seridata=seridata+temp1;     %得到加窗后nf个OFDM符号的串行数据采样点数
%    temp1=zeros(1,(ifftlen+postfix)); 
% end
% 
% %**************************************************************************
% %                           并串转换
% %**************************************************************************
% 
% Stx_data=seridata;
% %**************************************************************************
% %                           并串转换
% %**************************************************************************

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] = modulation(Stx_data,t,fp);
Stx_ofdm=Stx_data.*exp(j*2*pi*fp*t);
%Stx_ofdm=Stx_data;
%**************************************************************************
%                           加性高斯白噪声信道
%**************************************************************************
 snrdb=0;
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=Stx_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= 80Hz
datalen=8192;            % 选取处理数据长度
Ce=180;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算OFDM循环自相关函数

%  去掉 '%',可观察OFDM信号的循环自相关
% Ttao=80;
% TTtao=Ttao/2;
% t=[1:datalen]*Ts;
% Fs=fs/datalen;             % 频率分辨率
% af=[-511:1:512]*2*Fs;      % 循环频率搜索范围
% R_alpha=zeros(Ttao+1,datalen);    % 定义循环自相关矩阵
% for tao=-TTtao:1:TTtao       %  t延迟范围：(-30:30)
%     R_tao=redata(Ce-tao:datalen+Ce-1-tao).*conj(redata(Ce+tao:Ce+datalen+tao-1));   %计算BPSK信号的自相关函数：每一次的延迟为2tao
%     R_alpha(tao+TTtao+1,:)=1/datalen*fft(R_tao); %计算每一次延迟下的BPSK循环自相关函数
% end
% size(R_alpha);
% figure                            
% [x,y]=meshgrid([1:datalen]*Fs,(-80:2:80)*Ts);
% surf(x,y,abs(R_alpha))
% colormap([1 1 1]);                             %图形变黑
% xlabel('α [Hz]');
% ylabel('time delay [s]');
% %zlim([0 0.25])
% zlabel('Rx(α)')
% title('CAF of OFDM signal');      %  观察BPSK信号的循环自相关图
% deletf=1/Tu;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求循环谱中 f=0的截面

S_alpha=zeros(1,1024);
for i=1:10
re_x2=redata(Ce+(i-1)*datalen:Ce+datalen-1+(i-1)*datalen);
X=fft(re_x2);
Fs=fs/datalen;                      %频域分辨率
L=501;
S_alpha1=zeros(1,1024);
%plot(abs(fftshift(X)))
for m=-511:1:512    %循环频率范围
temp2=X((512+1+m):(512+1+m+L-1));  %共取L点
temp3=X((512+1-m):(512+1-m+L-1));
temp4=sum(temp2.*conj(temp3))/(501*8192);
S_alpha1(1,(m+511+1))=temp4;
end
S_alpha=abs(S_alpha1)+S_alpha;
end
S_alpha2=S_alpha/10;
alphalen=[-511:1:512]*2*Fs;
figure
plot(alphalen,abs(S_alpha2))

% 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=[0:0.5: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:10
%  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/10;
% % %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(:,1);
% 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;