ofdmsubspace.m

一个对OFDM信道进行盲估计的子空间算法

clear;
clc;
% ************************************************ 参数设置及初始化 ********************************************************

Fd = 1;                            % 消息序列的采样速率
Fs = 1*Fd;                         % 已调信号的采样速率
Md = 2;                            % 调制进制数,如: Md=2,Md_PSK=BPSK; Md=4,Md_PSK=QPSK

FrameNum = 50;                     % 原始数据的总帧数
N = 16;                            % 原始OFDM数据帧长度（子载波个数）
DataLen = FrameNum*N;              % 原始数据的总长度320

EbN0_Begin_dB = 0;                % 仿真初始Eb/N0值
EbN0_End_dB = 30                  % 仿真终止Eb/N0值
EbN0_Step_dB = 5;                  % 仿真Eb/N0步进值  
EbN0_dB = [ EbN0_Begin_dB : EbN0_Step_dB : EbN0_End_dB ];                     
EbN0 = 10.^(EbN0_dB/10);                                                      

Xn = zeros(1,DataLen);             % 原始数据
Xd = zeros(1,DataLen);             % 调制后的数据

CycleNum =100;                    % Monte Carlo仿真循环测试次数
ChType=1;                          % 选择信道类型

% ******************************************************* 多径信道 ************************************************************

% 多径衰落信道    
 switch ChType                                                           
     case 1
          % 多径衰落信道1：IEEE802.20 Vehicular test Channel-A 62.5ns/sample（M.1225标准）
          SampR = 200;                                    % 信号速率 200ns/sample
          ChIR = zeros(round(1090/SampR)+1,1);            % 信道冲激响应初始化
            
          ChIR(1) = 1;                                    % ChIR=[1.0000  0  0.7943  0  0.1259  0.1000]
          ChIR(round(310/SampR)+1) = 10^((-1.0)/10);      % 信道长度：6 chips
          ChIR(round(710/SampR)+1) = 10^((-9.0)/10); 
          ChIR(round(1090/SampR)+1) = 10^((-10.0)/10);
      
     case 2  
           % 多径衰落信道2：自定义的各径能量相等的多径衰落信道 
          ChIR = zeros(16,1);

          ChIR(1) = 0.4;
          ChIR(4) = 0.4;
          ChIR(7) = 0.4;
          ChIR(10) = 0.4;
          ChIR(13) = 0.4;
          ChIR(16) = 0.4;
     case 3  
           % 多径衰落信道3：自定义的多径衰落信道 
          ChIR = zeros(4,1);

          ChIR(1) = 1;
          ChIR(2) = 0.5;
          ChIR(3) = 0.3;
          ChIR(4) = 0.1;
      otherwise
        error('method参数设置错误！')
end

% 信道能量归一化
a = 1/sqrt(sum(ChIR.^2));     % 归一化因子
%ChIR=a*ChIR;          

est_ChIR=zeros(length(ChIR),length(EbN0_dB));  % 信道冲激响应估计值向量初始化  4*31

% 信道冲激响应的Toeplitz矩阵 N=16 L=4
H=toeplitz([ChIR' zeros(1,N-1)],[ChIR(1) zeros(1,N-1)]); % 由信道冲激响应生成toeplitz矩阵，(N+L-1)*N 维 19*16

% ****************************************************** 生成DFT矩阵 ************************************************************

DFT_mat=zeros(N,N);         % 初始化DFT矩阵 N*N 16*16
for k=1:N
    for n=1:N
        DFT_mat(k,n)=1/sqrt(N)*exp(-j*2*pi*(k-1)*(n-1)/N);
    end
end

% *******************************************************************************************************************************

for i = 1:length(EbN0_dB)                                               % 测试信噪比的序号 1:30                               
tic                                                                     % 程序执行时间测试
SNR = EbN0_dB(i)                                                        % 信噪比(dB)显示
sum_est_ChIR=zeros(length(ChIR),1);                                     % 信道估计值的和矩阵初始化 4*1
    for rr = 1:CycleNum                                                 % 每个SNR值测试多次，用以求均值和方差

% ******************************************************** 发射端 ****************************************************************

        Xn = randint(1,DataLen,Md);                                        % 随机产生原始信号
        Xd = dmodce(Xn,Fd,Fs,'psk',Md);                                    % 生成Md进制的PSK调制信号
  
     % 计算噪声方差
        Ps = sum(abs(Xd).^2)/DataLen;                                     % 整个传输信号的平均功率  
        Sigma = sqrt( ( abs(ChIR(1))/sqrt(N))^2*Ps/(2*EbN0(i)) );   
                                                       % 根据信噪比计算信道的噪声均方差（由于是复高斯噪声，所以再除以sqrt(2)）

% ************************************************************ 接收端 ******************************************************************
        
        sum_cor_Rsignal=zeros( N+length(ChIR)-1,N+length(ChIR)-1 );   % 初始化多个OFDM块的接收自相关矩阵的和19*19 (N+L-1)*(N+L-1)
        sum_cor_V_mat=zeros( length(ChIR),length(ChIR) );             % 4*4                                        L*L
        for n = 1:FrameNum
            Noise = normrnd(0,Sigma,1,N+length(ChIR)-1)+j*normrnd(0,Sigma,1,N+length(ChIR)-1);  % 1*19     
            Rsignal = H*DFT_mat'*Xd(1+(n-1)*N:n*N)'+Noise';        % 接收信号 r=HF'x+v
%             Rsignal = H*DFT_mat'*Xd(1+(n-1)*N:n*N)';        % 接收信号
            
            cor_Rsignal=Rsignal*Rsignal';                         % Rsignal的自相关阵   R=r*r'
            sum_cor_Rsignal=sum_cor_Rsignal+cor_Rsignal;          % 多个OFDM块的接收自相关矩阵的和

        end
            cor_Rsignal=(1/FrameNum)*(sum_cor_Rsignal);       % 对r的自相关矩阵的和求平均
            [U S V]=svd(cor_Rsignal);                         % 对Rsignal的自相关阵进行奇异值分解
            V_Hankel=hankel( [U(1:N+length(ChIR)-1,N+length(ChIR)-1)' U(1,N+length(ChIR)-1)],U(1:N,N+length(ChIR)-1));  
                                                              % 利用奇异值分解的最小特征值对应的特征向量生成henkel矩阵
            V_mat=V_Hankel(1:length(ChIR),:);                 % 取hankel矩阵的第1到第L行
            cor_V_mat=V_mat*V_mat';                           % 矩阵V_mat的自相关矩阵
            [U1 S1 V1]=svd(cor_V_mat);                        % 对cor_V_matl的自相关阵进行奇异值分解
            est_ChIR(:,i)=abs((1/a)*U1(1:length(ChIR),length(ChIR)) ) ;      % 信道冲激响应的估计值
            %est_ChIR(:,i)=est_ChIR(:,i)/est_ChIR(1)
            SE(rr)=sum((est_ChIR(:,i)-ChIR).^2)/sum((ChIR).^2);           % 每次Monte Carlo仿真得到的信道估计值与真实值的方差
    end

            MSE_ChIR(i)=(1/(CycleNum))*sum(SE);                  % 信道估计值的均方误差
            T=toc   
 
end
%  plot(1:6,est_ChIR,'r+-',1:6,ChIR,'bo-')
%  xlabel('码片');
%  ylabel('能量');
%  legend('Estimated','Accurate')
%  title('子空间算法下信道冲击响应估计 15dB')
% grid on;
%绘制方差曲线
%figure(1);
hold on;
semilogy(EbN0_dB,MSE_ChIR,'k-');  
xlabel('Eb/N0(dB)');
ylabel('MSE');
grid on;