DCT-MCM的matlab仿真

最近看了几篇基于DCT的多载波调制有关论文，看来看去主要思想都是通过添加冗余（一般是对称扩展或零填充）和预滤波器（为了CIR对称），使信道传输矩阵可以被DCT对角化，感觉没有什么新东西。 针对其中《A_Novel_Scheme_of_Multicarrier_Modulation_With_the_Discrete_Cosine_Transform》中的一部分进行仿真，在固定信道、无信道估计、BPSK映射情况下，计算信号重建的误比特率。这篇论文里的MIRA属实没有理解透彻，直接按式子仿真了。

%固定信道
%导频子载波数量Np
%重建系数数量Nr
%DCT1e长度N=N0



%% 参数设置
L=11; %CIR长度
h_ch=[1,0,0,-0.5,0,0,0,0.25,0,0,0.05];%固定的CIR
N=128; %DCT长度，子载波数量
N_simu=100; %仿真次数
Ns=20;%MCM符号数量
Np=128;%导频子载波数量
SNR=-10:1:40;%信噪比范围
ber_data=zeros(length(SNR),N_simu);%存放每次仿真的ber
ber_avg=zeros(length(SNR),1);%平均ber
N_mod=2; %64QAM或者BPSK

%% 系统仿真
for snr_iter=1:length(SNR)  %不同信噪比下实验
    snr_db=SNR(snr_iter);
    snr=10^(snr_db/10);
   
    for simu_iter=1:N_simu   %仿真次数
        %发射端
        num=0;%不同的比特数
        for n=2:Ns %MCM符号
                bit=randi([0,1],1,(N-2)*log2(N_mod));%产生01序列，这里的N-2对应论文里的N-2个子载波
                X_bit=reshape(bit,N-2,[]);
                if N_mod==2
                    %BPSK映射
                    bit_moded=2*bit-ones(size(bit));%0、1 -> -1、1
                elseif N_mod==64
                    %64QAM映射
                    bit=reshape(bit,log2(N_mod),[])';%调整数组大小
                    b2d=bi2de(bit);%每4比特进行二进制转十进制
                    bit_moded=qammod(b2d,N_mod);%QAM调制
                end

                %不进行信道估计的情况
                X=reshape(bit_moded,N-2,[]);%将调制符号调整为N-2行
                col=length(bit_moded)/(N-2);%col等于1，只有一列，一个MCM符号
                X_tmp=zeros(N,col);
                X_tmp(2:N-1,:)=X;%第2到N-1行是原始数据，第1行和第N行根据已知数据计算得到
                for i_=1:col
                    tmp1=0;
                    tmp2=0;
                    for j_=1:(N/2-1)
                        tmp1=tmp1+X_tmp(2*j_+1,i_);%式子做了一点修改，下标对应
                        tmp2=tmp2+X_tmp(2*j_,i_);
                    end
                    X_tmp(1,i_)=-2*tmp1;
                    X_tmp(N,i_)=-2*tmp2;
                end
                x=dct_matrix(N)*X_tmp;%逆DCT变换，第一行和第N行数据接近0，后面舍弃这两行数据
                x_e=reshape(x,[],1);%并串转换
                y=conv(x_e,h_ch);  %卷积，长度是Ne+2L-2
                y_noise=awgn(y,snr_db,'measured');
                h_pf=flip(h_ch);%预滤波器是h_ch的反转
                y_tilda=conv(y_noise,h_pf);%再和预滤波器卷积，长度是N+2L-2
%                 y_tilda_new=zeros(N+2*L,1);
%                 y_tilda_new(2:N+2*L-1)=y_tilda; %不是很理解第五部分y_tilda长度是N+2L，是在y_tilda两端加0吗
                %MIRA程序
                MIRA=zeros(N,N+2*(L-1));
                MIRA2=zeros(N,N+2*(L-1));
                MIRA(:,L:L-1+N)=eye(N);
                J=eye(L-1);
                J=J(:,end:-1:1);
                MIRA2(1:L-1,2:L)=J;
                MIRA2(N-(L-1)+1:N,N+(L-1):N+2*(L-1)-1)=J;
                MIRA=MIRA+MIRA2;
                y_tilda_new=MIRA*y_tilda;
                
                y_r=dct_matrix(N)*y_tilda_new;
                %BPSK解调

                y_r=y_r(2:N-1,:);
                y_r(find(y_r>0))=1;
                y_r(find(y_r<0))=0;
                num=num+sum(sum(abs(y_r-X_bit)));
                
                %64QAM解调
%                  y_r=y_r(2:N-1,:);
%                  y_r=qamdemod(y_r,N_mod,'OutputType','bit');
%                 bit2=reshape(bit,[],1);
%                 num=num+sum(sum(bit2==y_r));
            end
        %有信道估计的情况
        ber_data(snr_iter,simu_iter)=num/((N-2)*Ns*N_mod);%对于BPSK
        
    end
    ber_avg(snr_iter)=sum(ber_data(snr_iter,:))/N_simu;   
end
semilogy(SNR,ber_avg);

%dct_matrix.m
function dct_m=dct_matrix(N)
%产生DCT变换矩阵，根据式(1)(2)
dct_m=zeros(N,N);
for row=1:N
    for col=1:N
        if col==1|| col==N
            a=1/(sqrt(2*(N-1)));      %论文中是N-1，但是matlab索引从1开始
        else
            a=2/(sqrt(2*(N-1)));
        end
        dct_m(row,col)=a*cos((row-1)*(col-1)*pi/(N-1));
    end
end

结果只是一个ber曲线罢了…