virtual.m

系统中

%虚拟子载波 条件：MJ-L>MP 
clear all
clc
close all
j=sqrt(-1); 
Q=15;P=13;e=300;P0=1;L=3;M=2;D=4;J=Q+D;%Q:子载波数 ，P:有用信号；J：加上循环前缀最后，M：联合估计符号数；D：循环前缀长度 

mse(1:9) = 0;
squ_err_h(1:9) = 0; 
%%%%%%%%%%%%%%%%%%%%%%%%%% 
for test_time = 1:50
    scale_factor = randn+randn*j;
    h(4)=exp(-3/0.6)*scale_factor;
    h(3)=exp(-2/0.6)*scale_factor;
    h(2)=exp(-1/0.6)*scale_factor;
    h(1)=exp(0/0.6)*scale_factor;
    h = h/sqrt(sum(h.*conj(h)));
    H=toeplitz([h(4);zeros(M*J-L-1,1)],[h(4) h(3) h(2) h(1) zeros(1,M*J-L-1)]); %（MJ-L）*MJ  
    dB_count = 0;
    for dB=0:5:40
        
        dB_count = dB_count + 1;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IDFT 矩阵
        w=exp(-j*2*pi/Q); 
        for i=1:Q 
            for s=1:Q 
            Wn(i,s)=1/sqrt(Q)*w^[(-1)*(i-1)*(s-1)]; 
            end  
        end  
        for i=1:Q 
            t(i)=w^[(-1)*P0*(i-1)];
        end 

        T=diag(t);%Q*Q 
        W=T*Wn(:,1:P);%Q*P 

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Tcp 
        Wp=W(Q-D+1:Q,:);%D*Q 
        Wcp=[Wp;W];%J*P 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        Im=eye(M); 
        G=kron(Im,Wcp);%MJ*MP G临时使用 发送端是M个符号
        A=H*G;%（MJ-L）*MP 
        for i1=1:e 
            %%%%%%%%%%%%%%%%%%%%%%%%%%% 
            X=rand(M*P,1);%生成MP维BPSK码 
            for m=1:M*P 
                if X(m,1)>(0.5) 
                    X(m,1)=1; 
                else
                    X(m,1)=-1;
                end 
            end 
            Y(:,i1)=A*X;%(MJ-L)*1  
        end
        Y = awgn(Y,dB,'measured');
        [s v d] = svd(Y);
        Un = s(:,M*P+1:end);

        GG = 0;
        for i = 1:M*(J-P)-L
            Un1 = Un(:,i).';
            Ui = toeplitz([Un1(1,1);zeros(L,1)],[Un1,zeros(1,L)]);
            GG = conj(Ui)*G*G'*Ui.'+GG;            
        end
        [c v d] = svd(GG);
        phase_est = angle(c(4,4)/h(1));
        h_est = c(:,4)*exp(-j*phase_est);
        h1 = [h(4) h(3) h(2) h(1)].';
        mse(dB_count)=(h1-h_est)'*(h1-h_est)/(h1'*h1)+mse(dB_count);
        squ_err_h(dB_count)=sqrt((h1-h_est)'*(h1-h_est))/sqrt(h1'*h1)/sqrt(L+1)+squ_err_h(dB_count);
    end
end

mse = mse/test_time;
squ_err_h = squ_err_h/test_time;

figure; 
semilogy(0:5:40,squ_err_h,'-*r');
grid on;
xlabel('SNR(dB)');ylabel('RMSE');

% figure;
% subplot(211);plot(abs(h),'-*r');hold on; plot(abs(hb),'-ob');grid on;
% legend('真实信道','估计信道');ylabel('幅度');
% title('盲估计算法 SNR=15dB下的结果比较');
% subplot(212);plot(angle(h),'-*r');hold on; plot(angle(hb),'-ob');grid on;
% legend('真实信道','估计信道');ylabel('相位（rad）');
% title('盲估计算法 SNR=15dB下的结果比较');
% 
% figure;
% subplot(211);plot(abs(h),'-*r');hold on; plot(abs(hb1),'-ob');grid on;
% legend('真实信道','估计信道');ylabel('幅度');
% title('半盲估计算法 SNR=15dB下的结果比较');
% subplot(212);plot(angle(h),'-*r');hold on; plot(angle(hb1),'-ob');grid on;
% legend('真实信道','估计信道');ylabel('相位(rad)');
% title('半盲估计算法 SNR=15dB下的结果比较');
%