function [ser,mse]=dem_DFUSC_4x2_STBC(demparams,NoisePwr,ipfid_rx1,ipfid_rx2)

PwNorm=sqrt(2048/1702);     % symbol power normalization factor
% load original data
switch demparams.qam,
    case 4,
        load orig_DFUSC_QPSK.mat;    
    case 16,
        load orig_DFUSC_QAM16.mat;
    case 64,
        load orig_DFUSC_QAM64.mat;
    otherwise,
        error('Not supported Modulation Schme!');
end
load Wk_preamble.mat;

    % get preamble modulation and index values
    preamble_index_tx1=0:4:1702;    % Tx 1 use k mod 4==0 subcarriers
    preamble_index_tx2=1:4:1702;    % Tx 2 use k mod 4==1 subcarriers
    preamble_index_tx3=2:4:1702;    % Tx 2 use k mod 4==2 subcarriers
    preamble_index_tx4=3:4:1702;    % Tx 2 use k mod 4==3 subcarriers
    preamble_mod_tx1=sqrt(2)*2*2*(0.5-Wk_preamble(1:4:end));
    preamble_mod_tx2=sqrt(2)*2*2*(0.5-Wk_preamble(2:4:end));
    preamble_mod_tx3=sqrt(2)*2*2*(0.5-Wk_preamble(3:4:end));
    preamble_mod_tx4=sqrt(2)*2*2*(0.5-Wk_preamble(4:4:end));  % Note: include DC. At No.213 element

    % ######## Perfect Channel Estimation ########
    switch demparams.chantype,
        case 1,
            load perfectCE_sui1_11.mat;
            load perfectCE_sui1_21.mat;
            load perfectCE_sui1_31.mat;
            load perfectCE_sui1_41.mat;
            load perfectCE_sui1_12.mat;
            load perfectCE_sui1_22.mat;
            load perfectCE_sui1_32.mat;
            load perfectCE_sui1_42.mat;
        case 2,
            load perfectCE_sui2_11.mat;
            load perfectCE_sui2_21.mat;
            load perfectCE_sui2_31.mat;
            load perfectCE_sui2_41.mat;
            load perfectCE_sui2_12.mat;
            load perfectCE_sui2_22.mat;
            load perfectCE_sui2_32.mat;
            load perfectCE_sui2_42.mat;
        case 3,
            load perfectCE_sui3_11.mat;
            load perfectCE_sui3_21.mat;
            load perfectCE_sui3_31.mat;
            load perfectCE_sui3_41.mat;
            load perfectCE_sui3_12.mat;
            load perfectCE_sui3_22.mat;
            load perfectCE_sui3_32.mat;
            load perfectCE_sui3_42.mat;
        case 4,
            load perfectCE_sui4_11.mat;
            load perfectCE_sui4_21.mat;
            load perfectCE_sui4_31.mat;
            load perfectCE_sui4_41.mat;
            load perfectCE_sui4_12.mat;
            load perfectCE_sui4_22.mat;
            load perfectCE_sui4_32.mat;
            load perfectCE_sui4_42.mat;
        case 5,
            load perfectCE_sui5_11.mat;
            load perfectCE_sui5_21.mat;
            load perfectCE_sui5_31.mat;
            load perfectCE_sui5_41.mat;
            load perfectCE_sui5_12.mat;
            load perfectCE_sui5_22.mat;
            load perfectCE_sui5_32.mat;
            load perfectCE_sui5_42.mat;
        case 6,
            load perfectCE_sui6_11.mat;
            load perfectCE_sui6_21.mat;
            load perfectCE_sui6_31.mat;
            load perfectCE_sui6_41.mat;
            load perfectCE_sui6_12.mat;
            load perfectCE_sui6_22.mat;
            load perfectCE_sui6_32.mat;
            load perfectCE_sui6_42.mat;
        otherwise,
            error('Unknown SUI channel type!');
    end
    perfectCE_11(:,1877:2048)=[]; perfectCE_11(:,1:173)=[];
    perfectCE_21(:,1877:2048)=[]; perfectCE_21(:,1:173)=[];
    perfectCE_31(:,1877:2048)=[]; perfectCE_31(:,1:173)=[];
    perfectCE_41(:,1877:2048)=[]; perfectCE_41(:,1:173)=[];
    perfectCE_12(:,1877:2048)=[]; perfectCE_12(:,1:173)=[];
    perfectCE_22(:,1877:2048)=[]; perfectCE_22(:,1:173)=[];
    perfectCE_32(:,1877:2048)=[]; perfectCE_32(:,1:173)=[];
    perfectCE_42(:,1877:2048)=[]; perfectCE_42(:,1:173)=[];
    % ################################################################
    
suberr=0;   % used to calculate subcarrier error rate
mse=0;
datanum=demparams.numframe*2*demparams.sympf*1702;  % data number

for frm=1:demparams.numframe
    % ############## When use Known Channel Knowledge #############
    FullChanEst_11=sqrt(2048/1702)*perfectCE_11(frm,:);
    FullChanEst_21=sqrt(2048/1702)*perfectCE_21(frm,:);
    FullChanEst_31=sqrt(2048/1702)*perfectCE_31(frm,:);
    FullChanEst_41=sqrt(2048/1702)*perfectCE_41(frm,:);
    FullChanEst_12=sqrt(2048/1702)*perfectCE_12(frm,:);
    FullChanEst_22=sqrt(2048/1702)*perfectCE_22(frm,:);
    FullChanEst_32=sqrt(2048/1702)*perfectCE_32(frm,:);
    FullChanEst_42=sqrt(2048/1702)*perfectCE_42(frm,:);
    % #############################################################
    
    % use preamble for channel estimation
       %---------- Path 11, 21, 31, 41 -------------
    [preamble_time_rx1, symlen] = fread(ipfid_rx1, 2*2048*(1+demparams.Guard), 'float32');
    if symlen~=2*2048*(1+demparams.Guard), error('fread Rx 1 input data error!'); end
    preamble_time_rx1 =  [ 1 1j ] * reshape(preamble_time_rx1, 2, symlen/2);
    symlen=symlen/2;
    % remove CP
    preamble_time_rx1=preamble_time_rx1((2048*demparams.Guard+1):end);
            % ###### Time domain ls estimation ######
      %      [FullChanEst_11,FullChanEst_21,FullChanEst_31,FullChanEst_41] ...
      %      = ls_time_4x2((preamble_time_rx1).',demparams.chantype,demparams.sampletime);
      %      FullChanEst_11=PwNorm*FullChanEst_11; FullChanEst_21=PwNorm*FullChanEst_21;
      %      FullChanEst_31=PwNorm*FullChanEst_31; FullChanEst_41=PwNorm*FullChanEst_41;
            % #######################################
    % FFT
 %   preamble_freq_rx1=fftshift(fft(preamble_time_rx1,2048)/sqrt(2048));
    % Trim Guard subcarriers
 %   preamble_freq_rx1(1877:2048)=[];  preamble_freq_rx1(1:173)=[];
            % ###### Extract preamble pilots AND Channel Estimation ######
 %           preamble_pilots_tx1=preamble_freq_rx1(preamble_index_tx1+1);
 %           preamble_pilots_tx2=preamble_freq_rx1(preamble_index_tx2+1);
 %           preamble_pilots_tx3=preamble_freq_rx1(preamble_index_tx3+1);
 %           preamble_pilots_tx4=preamble_freq_rx1(preamble_index_tx4+1);
 %           ChanEst_11=ls_freq(preamble_pilots_tx1,preamble_mod_tx1);
 %           ChanEst_21=ls_freq(preamble_pilots_tx2,preamble_mod_tx2);
 %           ChanEst_31=ls_freq(preamble_pilots_tx3,preamble_mod_tx3);
 %           ChanEst_41=ls_freq(preamble_pilots_tx4,preamble_mod_tx4);
            % Interpolate  Note: If interpolate 'Tx4', must remove DC
 %           FullChanEst_11=interpolate_4x2(preamble_index_tx1+1,ChanEst_11,'linear','Tx1');
 %           FullChanEst_21=interpolate_4x2(preamble_index_tx2+1,ChanEst_21,'linear','Tx2');
 %           FullChanEst_31=interpolate_4x2(preamble_index_tx3+1,ChanEst_31,'linear','Tx3');
 %           FullChanEst_41=interpolate_4x2(preamble_index_tx4+1,ChanEst_41,'linear','Tx4');
            % ############################################################
            
       %------------ Pathe 12, 22, 32, 42-----------
    [preamble_time_rx2, symlen] = fread(ipfid_rx2, 2*2048*(1+demparams.Guard), 'float32');
    if symlen~=2*2048*(1+demparams.Guard), error('fread Rx 2 input data error!'); end
    preamble_time_rx2 =  [ 1 1j ] * reshape(preamble_time_rx2, 2, symlen/2);
    symlen=symlen/2;
    % remove CP
    preamble_time_rx2=preamble_time_rx2((2048*demparams.Guard+1):end);
            % ###### Time domain ls estimation ######
      %      [FullChanEst_12,FullChanEst_22,FullChanEst_32,FullChanEst_42] ...
      %      = ls_time_4x2((preamble_time_rx2).',demparams.chantype,demparams.sampletime);
      %      FullChanEst_12=PwNorm*FullChanEst_12; FullChanEst_22=PwNorm*FullChanEst_22;
      %      FullChanEst_32=PwNorm*FullChanEst_32; FullChanEst_42=PwNorm*FullChanEst_42;
            % #######################################
    % FFT
 %   preamble_freq_rx2=fftshift(fft(preamble_time_rx2,2048)/sqrt(2048));
    % Trim Guard subcarriers
 %   preamble_freq_rx2(1877:2048)=[];  preamble_freq_rx2(1:173)=[];
            % ###### Extract preamble pilots AND Channel Estimation ######
 %           preamble_pilots_tx1=preamble_freq_rx2(preamble_index_tx1+1);
 %           preamble_pilots_tx2=preamble_freq_rx2(preamble_index_tx2+1);
 %           preamble_pilots_tx3=preamble_freq_rx2(preamble_index_tx3+1);
 %           preamble_pilots_tx4=preamble_freq_rx2(preamble_index_tx4+1);
 %           ChanEst_12=ls_freq(preamble_pilots_tx1,preamble_mod_tx1);
 %           ChanEst_22=ls_freq(preamble_pilots_tx2,preamble_mod_tx2);
 %           ChanEst_32=ls_freq(preamble_pilots_tx3,preamble_mod_tx3);
 %           ChanEst_42=ls_freq(preamble_pilots_tx4,preamble_mod_tx4);
            % Interpolate  Note: If interpolate 'Tx4', must remove DC
 %           FullChanEst_12=interpolate_4x2(preamble_index_tx1+1,ChanEst_12,'linear','Tx1');
 %           FullChanEst_22=interpolate_4x2(preamble_index_tx2+1,ChanEst_22,'linear','Tx2');
 %           FullChanEst_32=interpolate_4x2(preamble_index_tx3+1,ChanEst_32,'linear','Tx3');
 %           FullChanEst_42=interpolate_4x2(preamble_index_tx4+1,ChanEst_42,'linear','Tx4');
            % ############################################################
            
        % ###### MSE Calculation ######