sfofdm2.asv

block_mimo_ ofdm_channelestimation编码 很好的程序

function [ber]=sfofdm2(H,delays,out)

scmpar=scmparset;
NumLoop = 10;
NumSubc = 1024;
NumCP =40;

 delta_t = out.delta_t;
 f_sym = 3.84e6;                                       % sample frequency                                       
 osr = floor((1/f_sym)/scmpar.DelaySamplingInterval);  % over sample rate
 ts = 1/(f_sym*osr);                                   % sample time_step
 h_delay = round(delays./ts)+1;                        % multipath delay array: samples
 u = scmpar.NumMsElements;
 s = scmpar.NumBsElements;
 n = scmpar.NumPaths;
 h_tmp = zeros(1,NumSubc*osr);
 h = zeros(s,u,NumSubc);
 for ss = 1:s
     for uu = 1:u
             h_tmp(h_delay(1:n)) = H(uu,ss,1:n,1);
             h(ss,uu,:) = resample(h_tmp,1,osr); % downsample by osr to
     end
 end
 h11 = squeeze(h(1,1,:)); H11 = fft(h11);
 h12 = squeeze(h(1,2,:)); H12 = fft(h12).';
 h21 = squeeze(h(2,1,:)); H21 = fft(h21);
 h22 = squeeze(h(2,2,:)); H22 = fft(h22).';



% Generate the random binary stream for transmit test
BitsTx = floor(rand(1,2*NumLoop*NumSubc)*2);

% Modulate (generates QAM symbols)
% input: BitsTx(1,NumLoop*NumSubc); output: SymQAM(NumLoop,NumSubc/2)
SymQAMtmp = reshape(BitsTx,NumLoop*NumSubc,2);   %64000行2列
SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');        %64000行1列 把二进制数转换成十进制数
%00->-1-i,01->-1+i,10->1-i,11->1+i
QAMTable = [-1-i -1+i 1-i 1+i];
SymQAM = QAMTable(SymQAMtmptmp+1)/2;   %64000行1列

%A=size(SymQAM,1)*size(SymQAM,2)=1024000
% for ii = 1:2:NumLoop*NumSubc/2
%     T1_data(ii) = SymQAM(ii);
%     T2_data(ii) = SymQAM(ii+1);
%     T1_data(ii+1) = -conj(SymQAM(ii+1));
%     T2_data(ii+1) = conj(SymQAM(ii));
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%STBC编码
T1_data = zeros(NumSubc,NumLoop);
T2_data = zeros(NumSubc,NumLoop);
SymIFFTtmp = reshape(SymQAM,NumSubc,NumLoop);
for ii=1:2:NumSubc
    T1_data(ii,:)=SymIFFTtmp(ii,:);
    T2_data(ii,:)=SymIFFTtmp(ii+1,:);
    T1_data(ii+1,:)=-conj(SymIFFTtmp(ii+1,:));
    T2_data(ii+1,:)=conj(SymIFFTtmp(ii,:));
end
%    a11=size(T1_data,1) %1024 
%    a12=size(T1_data,2) %10
%    a21=size(T2_data,1)  %1024
%    a22=size(T2_data,2) %10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DFT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%小波
%%%%  WPDM ====================================
dep=3;
[t1_data,lextop1]=moduidwpt(T1_data,dep);
[t2_data,lextop2]=moduidwpt(T2_data,dep);

% a11=size(t1_data,1)
% a12=size(t1_data,2)
% a21=size(t2_data,1)
% a22=size(t2_data,2)
% % a11=2^dep%8
% T1_data=reshape(T1_data,size(T1_data,1)*size(T1_data,2)/(2^dep),2^dep);
% T2_data=reshape(T2_data,size(T2_data,1)*size(T2_data,2)/(2^dep),2^dep);
% chn1=[2^dep-1:2^(dep+1)-2]; % dep=3:[7:14];dep=2:[3:6]
% chn2=[2^dep-1:2^(dep+1)-2]; 
% dwtmode('per');    
%   lextop1= size(T1_data,1)*(2^dep); %小波变换需要的点数, 实际上是星座映射后的符号数
%   lextop2= size(T2_data,1)*(2^dep); 
%   t1 = wptree(2,dep,zeros(lextop1,1),'db4');
%   t2 = wptree(2,dep,zeros(lextop2,1),'db4');
% %  %sizet1=size(t,1)%1
% %  %sizet2=size(t,2)%1
% %  %XXXXX  
% %  
%          for wavii=1:2^dep
%             t1 = write(t1,'cfs',chn1(wavii),T1_data(:,wavii));  %分别返回节点chn（ii）的系数
%             t2 = write(t2,'cfs',chn2(wavii),T2_data(:,wavii));
%   	   end
%  	t1_data=read(wpjoin(t1),'data')';  %把小波系数的数据读出来
%       t2_data=read(wpjoin(t2),'data')'; 
% %       a11=size(t1_data,1)%10240
% %       a12=size(t1_data,2)%1
% %      %sizechant1=size(ChanT,1)%51200 25600 12800 12800
% %      %sizechant2=size(ChanT,2)%1     1     1     1
% %  	t1_data = upcoef('a',t1_data,'db4',6);%filtering
% %     t2_data = upcoef('a',t2_data,'db4',6);
% %      %sizechant1=size(ChanT,1)%3277178 1638778 819578 819578 
% %      %sizechant2=size(ChanT,2)%1       1       1      1
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%小波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DFT
%  t1_data1 = ifft(T1_data,NumSubc,1);
%   t2_data2 = ifft(T2_data,NumSubc,1);
% 
% %    a11=size(t1_data1,1)  
% %    a12=size(t1_data1,2)
% %    a21=size(t2_data2,1)  
% %    a22=size(t2_data2,2)
% 
%   
%   % Add cyclic prefix
%   NumAddPrefix = NumSubc + NumCP;
%   RowPrefix = (NumSubc - NumCP + 1):NumSubc;  % 复制
%   t1_data = [t1_data1(RowPrefix,:);t1_data1];
%   t2_data = [t2_data2(RowPrefix,:);t2_data2];
% 
% %     a11=size(t1_data,1)  
% %    a12=size(t1_data,2)
% %    a21=size(t2_data,1)  
% %    a22=size(t2_data,2)
% 
%   
%   t1_data=t1_data(:);
%   %a11=size(t1_data,1)%10640
%   %a12=size(t1_data,2)%1
%   t2_data=t2_data(:);
%   %a21=size(t2_data,1)%10640
%   %a22=size(t2_data,2)%1
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DFT
% block11 = t1_data;
% block12 = t2_data; 
% block21 = t1_data;
% block22 = t2_data; 

block11 = filter(h11,1,t1_data);
%  a11=size(block11,1)%10640
%  a12=size(block11,2)%1
block12 = filter(h12,1,t2_data); 
%  a21=size(block12,1)%10640
%  a22=size(block12,2)%1
block21 = filter(h21,1,t1_data);   
%  a31=size(block21,1)%10640
%  a32=size(block21,2)%1
block22 = filter(h22,1,t2_data);  
%  a41=size(block22,1)%10640
%  a42=size(block22,2)%1
       
% Add the AWGN
snr=0:1:25;% = SNR + 10*log10(log2(2));
M=length(snr);
ber=zeros(M,1);
for kk=1:M 
   kk
    nblock11 = awgn(block11,snr(kk),'measured');
%     a11=size(nblock11,1)%10640
%     a12=size(nblock11,2)%1
    nblock12 = awgn(block12,snr(kk),'measured');
%       a21=size(nblock12,1)%10640
%       a22=size(nblock12,2)%1
    nblock21 = awgn(block21,snr(kk),'measured');
%       a31=size(nblock21,1)%10640
%       a32=size(nblock21,2)%1
    nblock22 = awgn(block22,snr(kk),'measured');
%       a41=size(nblock22,1)%10640
%       a42=size(nblock22,2)%1


    re_block1 = nblock11 + nblock12;
%      a11=size(re_block1,1)%10640
%      a12=size(re_block1,2)%1
    re_block2 = nblock21 + nblock22;
%      a21=size(re_block2,1)%10640
%      a22=size(re_block2,2)%1
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%dft解调
%       SymDeCP1 = reshape(re_block1,NumAddPrefix,NumLoop);
%       SymDeCP2 = reshape(re_block2,NumAddPrefix,NumLoop);
%       SymDe1 = SymDeCP1((NumCP+1):NumAddPrefix,:);
%       SymDe2 = SymDeCP2((NumCP+1):NumAddPrefix,:);
% %       a11=size(SymDe1,1)%1024
% %       a12=size(SymDe1,2)%10
% %       a21=size(SymDe2,1)%1024
% %       a22=size(SymDe2,2)%10
%       N_Data1 = fft(SymDe1,NumSubc,1);
% %       a11=size(N_Data1,1)%1024
% %        a12=size(N_Data1,2)%10
%       N_Data2 = fft(SymDe2,NumSubc,1);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%小波包解调
 %%%%  WPDM ====================================
%  dep=3;
%  chn=[2^dep-1:2^(dep+1)-2];   %A Complete tree
%  dwtmode('per');   
%      re_block1= downcoef('a',re_block1,'db4',6);
%      re_block2= downcoef('a',re_block2,'db4',6);
%  	%ChanR= downcoef('a',ChanT,'db4',6);
%      re_block1= wkeep(re_block1,lextop1,'c');
%      re_block2= wkeep(re_block2,lextop2,'c');
%      %ChanR= wkeep(ChanR,lextop,'c');
%      %tr = wptree(2,dep,ChanR,'db4'); 
%      tr1 = wptree(2,dep,re_block1,'db4');
%      tr2 = wptree(2,dep,re_block2,'db4');
% 
%  		for iwavi=1:2^dep
%              N_Data1(:,iwavi)=read(tr1,'data',chn(iwavi));
%              N_Data2(:,iwavi)=read(tr2,'data',chn(iwavi));
%          end
Datarx=demodudwpt(re_block1,lextop)
Datarx=demodudwpt(re_block1,lextop)
[N_Data1,N_Data2]=demodu_wpdm(re_block1,re_block2,lextop1,lextop2);
% a11=size(N_Data1,1)%81968
% a12=size(N_Data1,2)%8
% a21=size(N_Data2,1)%81968
% a22=size(N_Data2,2)%8
             N_Data1=reshape(N_Data1,NumSubc,size(N_Data1,1)*size(N_Data1,2)/NumSubc);
             N_Data2=reshape(N_Data2,NumSubc,size(N_Data2,1)*size(N_Data2,2)/NumSubc);
% %-------------------解调————————

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%小波包解调

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DFT OVER
    R_Data = zeros(NumSubc,NumLoop);
    for jj = 1:2:NumSubc
        R_Data(jj,:) = conj(H11(jj))*N_Data1(jj,:) + H12(jj+1)*conj(N_Data1(jj+1,:))+...
                       conj(H21(jj))*N_Data2(jj,:) + H22(jj+1)*conj(N_Data2(jj+1,:));
        R_Data(jj+1,:) = conj(H12(jj))*N_Data1(jj,:) - H11(jj+1)*conj(N_Data1(jj+1,:))+...
                         conj(H22(jj))*N_Data2(jj,:) - H21(jj+1)*conj(N_Data2(jj+1,:));
    end

    R_Data=R_Data(:);
    SymDec = zeros(NumLoop*NumSubc,2);

    for m = 1:NumLoop*NumSubc
 
        Real = real(R_Data(m));
        Imag = imag(R_Data(m));
        if( abs((Real -1)) < abs((Real +1  )) )
            SymDec(m,1) = 1;
        else
            SymDec(m,1) = 0;
        end
        if( abs((Imag -1)) < abs((Imag +1  )) )
            SymDec(m,2) = 1;
        else
            SymDec(m,2) = 0;
        end
    end
    BitsRx=SymDec(:)';




% Bit Error

   [Num,Ber] = symerr(BitsTx,BitsRx);
   ber(kk)=Ber
end