通信DHT变换

tic
% Select_Subcarrier 256-128-16
% ************************** 参数设置 ********************************
clear all;
clc;
nloop=3;
Number_Tx=1;
Number_Rx=1;
Number_subcarrier1=256;
CP_length1=Number_subcarrier1/4;
Mod_level=2;
Constellation_size=2^Mod_level;
minSNR=0;
maxSNR=40;
step=4;
numsymbol=3;
% ************************** 信道参数设置 ********************************
fd=133;
tstp=5*1e-8;
tau=1e-9*[0 30 150 31 370 710 1090 1730 2510];
pdb=[0 -1.5 -1.4 -3.6 -0.6 -9.1 -7.0 -12.0 -16.9];
comcase=combntns(1:8,4);
Casesize=nchoosek(8,4);
% ************************** 开始循环 ********************************
for i=1:length(Number_subcarrier1)
    Number_subcarrier=Number_subcarrier1(i);
    CP_length=CP_length1(i);
    %Casesize_1=nchoosek(Number_subcarrier,Number_subcarrier/2);
    %Sele_Bits=fix(log2(Casesize_1));
    Sele_Bits=Number_subcarrier/8*fix(log2(Casesize));
    index=1;
for SNR=minSNR:step:maxSNR
    Numerror=0;
    for a=1:nloop
% ************************** 数据映射 ********************************
        Mapdata=zeros(Number_subcarrier,numsymbol);
        randsize=randi(Casesize,1,Number_subcarrier/8);
        Randdata=randi(Constellation_size,Number_subcarrier/2,numsymbol)-1;
        Moddata=qammod(Randdata,Constellation_size);
        selectcase=selectsubcarriercase(Number_subcarrier,randsize,comcase)';
        Mapdata(selectcase,:)=Moddata;
% ************************** DHT ********************************
        Data_after_DHT=DHT_matrix*Mapdata;
        temp1 = DHT(Data_after_DHT(:,1));
        temp2 = DHT(Data_after_DHT(:,2));
        temp3 = DHT(Data_after_DHT(:,3));
% **************************IFFT ********************************
        IFFTdata1 = ifft(temp1);
        IFFTdata2 = ifft(temp2);
        IFFTdata3 = ifft(temp3);
% **************************IFFT ********************************
%        IFFTdata=ifft(Mapdata);
% ************************** 加CP ********************************
        AddCPdata=[IFFTdata1(Number_subcarrier-CP_length+1:end) IFFTdata1...
            IFFTdata2(Number_subcarrier-CP_length+1:end) IFFTdata2...
            IFFTdata3(Number_subcarrier-CP_length+1:end) IFFTdata3];
        Txsignal=reshape(AddCPdata,numsymbol*(Number_subcarrier+CP_length),1);
% ************************** 信道 ********************************
        rayChanObj=rayleighchan(tstp,fd,tau,pdb);
        FadedSignal= filter(rayChanObj,Txsignal); 
% ************************** 加噪声 ********************************
        SignalPower=sum(abs(Txsignal).^2)./(numsymbol*(Number_subcarrier+CP_length));
        Withnoisedata=awgn(FadedSignal,SNR,10*log10(SignalPower));
% ************************** 去CP ********************************
        Rxsignal=reshape(Withnoisedata,(Number_subcarrier+CP_length),numsymbol);
        MoveCPdata=Rxsignal(CP_length+1:end,:);
% ************************** FFT ********************************
        FFTdata1=fft(MoveCPdata);
        FFTdata11 = DHT(FFTdata1(:,1));
        FFTdata12 = DHT(FFTdata1(:,2));FFTdata13 = DHT(FFTdata1(:,3));
        FFTdata=[FFTdata11;FFTdata12;FFTdata13]';
%***************************信道估计 *******************************
        %***************************完美估计*******************************
        RemoveCPSignal_Cp=Rxsignal_Cp(CP_length+1:end,:);                           
        FFTSignal_cp=fft(RemoveCPSignal_Cp);  
        ChanFreqResponse2=FFTSignal_cp./Mapdata;
        %***************************实际估计*******************************
        %ChanFreqResponse_1=FFTdata./Mapdata;                                       
        %for ii=1:numsymbol
         %ChanFreqResponse2(:,ii)= ChanFreqResponse_1(:,1);                                                                         % estimated channel frequency response (firt colum is used to estimate the other colums)
        %end   
%***************************检测接收信号 *******************************
        Detecdata=FFTdata(:,1)';
        [~,Receiveddataside]=paixu(Detecdata,Number_subcarrier);
        Receiveddata=FFTdata(Receiveddataside',:);
        ChanFreqResponse=ChanFreqResponse2(Receiveddataside',:);
%***************************信号补偿 *******************************
        ChanCompenData=Receiveddata./ChanFreqResponse;
%***************************解调 *******************************   
        Demoddata=qamdemod(ChanCompenData,Constellation_size);     
%***************************计算ber *******************************
        
        Numerror=Numerror+Numerror_1;
    end
    ber=Numerror./((Sele_Bits+numsymbol*Number_subcarrier/2*Mod_level).*nloop);
    fprintf('%d\t%e\n',SNR,ber);
    BER(index)=ber;
    index=index+1;
end
semilogy(minSNR:step:maxSNR,BER,'rs-');
hold on;
end
xlabel('SNR(dB)');
ylabel('BER');
title('Spatial Modulation with OFDM 256-128-16');
grid on;
hold on;
toc