Jump to content

    

Разделение суммы 2-х сигналов, ковариации, функция правдоподобия - суть

Задача такая. Я ее сам ставлю, не подумайте, что не разобравшись обращаюсь, но хотел бы просто небольшое Ваше мнение услышать

У нас есть 1. Генератор данных для N пользователей- 2. OFDM-модем. 3. Канал связи с импульсной характеристикой. 4. Приемник. В приемнике метод максимума функции правдоподобия и итерации для разделения пользователей (IDMA).

Полностью рабочий код приемника прилагаю.

#include <stdio.h>  
#include <stdlib.h>
#include <mex.h>
#include <math.h>
#include "define.h"

#define __RRE                 prhs[0]  // 1
#define __RIMG                prhs[1]  // 2
#define __SIGMA               prhs[2]  // 3
#define __DATALEN             prhs[3]  // 4
#define __CHIPLEN             prhs[4]  // 5
#define __PREFIX              prhs[5]  // 6
#define __NUSER               prhs[6]  // 7
#define __ITNUM               prhs[7]  // 8
#define __RAYRE               prhs[8]  // 9
#define __RAYIMG              prhs[9]  // 10
#define __SPREADLEN           prhs[10] // 11
#define __SPREADSEQ           prhs[11] // 12
#define __SCRAMBLERULEIN      prhs[12] // 13

#define __APPLLROUT           plhs[0]  // 14

int        _NUser        = NUSER;
int        _DataLen    = DATALEN;
int        _SpreadLen    = SPREADLEN;
int        _ChipLen    = CHIPLEN;
int     _Prefix     = PREFIX;
int        _ItNum        = ITNUM;
int        _Block        = BLOCK;
int        _EbNoNum    = EBNONUM;
double    _EbNoStart    = EBNOSTART;
double    _EbNoStep    = EBNOSTEP;

void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l,int il) 
{
    int it,m,is,i,j,nv,l0;
    double p,q,s,vr,vi,poddr,poddi;
    for (it=0; it<=n-1; it++)
    { 
        m=it; is=0;
        for (i=0; i<=k-1; i++)
        { 
            j=m/2; is=2*is+(m-2*j); m=j;}
            fr[it]=pr[is]; fi[it]=pi[is];
        }
        pr[0]=1.0; pi[0]=0.0;
        p=6.283185306/(1.0*n);
        pr[1]=cos(p); pi[1]=-sin(p);
        if (l!=0) 
            pi[1]=-pi[1];
        for (i=2; i<=n-1; i++)
        { 
            p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
            s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
            pr[i]=p-q; pi[i]=s-p-q;
        }
        for (it=0; it<=n-2; it=it+2)
        { 
            vr=fr[it]; vi=fi[it];
            fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];
            fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];
        }
        m=n/2; nv=2;
        for (l0=k-2; l0>=0; l0--)
        { 
            m=m/2; nv=2*nv;
            for (it=0; it<=(m-1)*nv; it=it+nv)
                for (j=0; j<=(nv/2)-1; j++)
                { 
                    p=pr[m*j]*fr[it+j+nv/2];
                    q=pi[m*j]*fi[it+j+nv/2];
                    s=pr[m*j]+pi[m*j];
                    s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
                    poddr=p-q; poddi=s-p-q;
                    fr[it+j+nv/2]=fr[it+j]-poddr;
                    fi[it+j+nv/2]=fi[it+j]-poddi;
                    fr[it+j]=fr[it+j]+poddr;
                    fi[it+j]=fi[it+j]+poddi;
                }
        }
        if (l!=0)
            for (i=0; i<=n-1; i++)
            { 
                fr[i]=fr[i]/(1.0*n);
                fi[i]=fi[i]/(1.0*n);
            }
        if (il!=0)
            for (i=0; i<=n-1; i++)
            { 
                pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
                if (fabs(fr[i])<0.000001*fabs(fi[i]))
                { 
                    if ((fi[i]*fr[i])>0) 
                        pi[i]=90.0;
                    else 
                        pi[i]=-90.0;
                }
                else
                    pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
            }
    return;
}

void Receiver( double  rRe[CONVLEN(PREFIX)],         //*in 1
               double  rImg[CONVLEN(PREFIX)],        //*in 2
               double  sigma,                        // in 3
               int     _DataLen,                     // in 4
               int     _ChipLen,                     // in 5
               int     _Prefix,                      // in 6
               int     _NUser,                       // in 7
               int     _ItNum,                       // in 8
               double  RayRe[DENULLIFY(PREFIX)],     //*in 9
               double  RayImg[DENULLIFY(PREFIX)],    //*in 10
               int     _SpreadLen,                   // in 11
               double  SpreadSeq [SPREADLEN],        //*in 12
               int     ScrambleRuleIn [FULLCODDATA], //*in 13
               double  AppLlrOut [FULLDATA])         //*in 14 
{
    int i, j, m, nuser, it, index, _tmp = 0;
    double  s2 = sigma * sigma;

    double  Chip [NUSER][CHIPLEN]         = {0}; 
    double    RCPrRe[CHIPLEN/2]             = {0}; // RCP stands for remove cyclic prefix
    double    RCPrImg[CHIPLEN/2]            = {0};
    double  FFTrRe[CHIPLEN/2]             = {0}; // Received symbols after OFDM demodulation
    double  FFTrImg[CHIPLEN/2]            = {0};
    double  CovX[NUSER][CHIPLEN/2][2][2]  = {0};
    double  MeanLRe[NUSER][CHIPLEN/2]     = {0};
    double  MeanLImg[NUSER][CHIPLEN/2]    = {0};
    double    MeanrRe[CHIPLEN/2]            = {0};
    double    MeanrImg[CHIPLEN/2]           = {0};
    double    MeanXRe[NUSER][CHIPLEN/2]     = {0}; //mean value for each chip
    double    MeanXImg[NUSER][CHIPLEN/2]    = {0};
    double  MeanLHRe[NUSER][CHIPLEN/2]    = {0};
    double  MeanLHImg[NUSER][CHIPLEN/2]   = {0};
    double  rkHRe[NUSER][CHIPLEN/2]       = {0}; //rkH means rk head
    double  rkHImg[NUSER][CHIPLEN/2]      = {0};
    double  Temp[NUSER][CHIPLEN/2][2][2]  = {0}; //for some matrix calculation
    double  CovL[NUSER][CHIPLEN/2][2][2]  = {0};
    double  Covr[CHIPLEN/2][2][2]         = {0};
    double  CovLH[NUSER][CHIPLEN/2][2][2] = {0}; //LH means L head
    double  ChipRe[NUSER][CHIPLEN/2]      = {0};
    double  ChipImg[NUSER][CHIPLEN/2]     = {0};
    double  TempChip[NUSER][CHIPLEN]      = {0};
    double    Ext[NUSER][CHIPLEN]           = {0};        
    double  R[NUSER][CHIPLEN/2][2][2]     = {0};
    double  RT[NUSER][CHIPLEN/2][2][2]    = {0}; //T means transpose

    double HRe [NUSER][CHIPLEN/2]         = {0};
    double HImg [NUSER][CHIPLEN/2]        = {0};
    
    double _HRe [CHIPLEN/2]               = {0};
    double _HImg [CHIPLEN/2]              = {0};
    double  AppLlr[NUSER][DATALEN]        = {0}; 
    int        ScrambleRule[NUSER][CHIPLEN]  = {0};

    for( nuser=0; nuser<_NUser; nuser++ ) 
        for( i=0; i<_ChipLen; i++ ) 
         ScrambleRule[nuser][i] = ScrambleRuleIn[(nuser * _ChipLen) + i];

    for( i=0; i<_ChipLen/2; i++ )
        {
            _HRe[i] = 0;
            _HImg[i] = 0;
        }

    
        for( i=0; i<_ChipLen/2; i++ )
        {
            if ( _Prefix == 0 )
            {
                j=0;
                _HRe[i] += RayRe[0]*cos(6.283185307*i*j/_ChipLen*2) + RayImg[0]*sin(6.283185307*i*j/_ChipLen*2);
                _HImg[i] += -RayRe[0]*sin(6.283185307*i*j/_ChipLen*2) + RayImg[0]*cos(6.283185307*i*j/_ChipLen*2);
            }

            if ( _Prefix > 0 )
            {
                for( j=0; j<DENULLIFY(_Prefix); j++ )
                {
                    _HRe[i] += RayRe[j]*cos(6.283185307*i*j/_ChipLen*2) + RayImg[j]*sin(6.283185307*i*j/_ChipLen*2);
                    _HImg[i] += -RayRe[j]*sin(6.283185307*i*j/_ChipLen*2) + RayImg[j]*cos(6.283185307*i*j/_ChipLen*2);
                }
            }           
        }

    for( nuser=0; nuser<_NUser; nuser++ )
        for( i=0; i<(_ChipLen/2); i++ ) 
        {
            HRe  [nuser][i] = _HRe[i];
            HImg [nuser][i] = _HImg[i];
        }
       

    for( nuser=0; nuser<_NUser; nuser++ ) 
    {
        for( i=0; i<_ChipLen/2; i++ )
        {
            R[nuser][i][0][0]  =  HRe[nuser][i];
            R[nuser][i][0][1]  = -HImg[nuser][i];
            R[nuser][i][1][0]  =  HImg[nuser][i];
            R[nuser][i][1][1]  =  HRe[nuser][i];
            RT[nuser][i][0][0] =  R[nuser][i][0][0];
            RT[nuser][i][0][1] =  R[nuser][i][1][0];
            RT[nuser][i][1][0] =  R[nuser][i][0][1];
            RT[nuser][i][1][1] =  R[nuser][i][1][1];
        }
    }


    // Remove cyclic perfix & cyclic
    for ( i=0; i<_ChipLen/2; i++ )
    {
        RCPrRe[i] = rRe[i+_Prefix];
        RCPrImg[i] = rImg[i+_Prefix];
    }

    // Get index for FFT and IFFT, where 2^index=_ChipLen/2 (расчет log2()) очень интересный способ)
    i = _ChipLen/2;
    index=0;
    while (i != 1)
    {
        i = i/2;
        index += 1;
    }    

    // FFT
    kkfft (RCPrRe,RCPrImg,_ChipLen/2,index,FFTrRe,FFTrImg,0,0);

    // Normalization
    for ( i=0; i<(_ChipLen/2); i++ )
    {
        FFTrRe[i] = FFTrRe[i]/pow(_ChipLen/2,0.5);
        FFTrImg[i] = FFTrImg[i]/pow(_ChipLen/2,0.5);
    }
        
    for( nuser=0; nuser<_NUser; nuser++ ) 
        for( i=0; i<(_ChipLen/2); i++ )
        {
            MeanXRe[nuser][i] = 0;
            MeanXImg[nuser][i] = 0;
        }

    for( i=0; i<(_ChipLen/2); i++ )
    {
        MeanrRe[i] = 0;
        MeanrImg[i] = 0;
    }
    for( i=0; i<(_ChipLen/2); i++ ) 
        for( nuser=0; nuser<_NUser; nuser++ )
        {
            MeanrRe[i] += ( HRe[nuser][i]*MeanXRe[nuser][i] - HImg[nuser][i]*MeanXImg[nuser][i] );
            MeanrImg[i] += ( HRe[nuser][i]*MeanXImg[nuser][i] + HImg[nuser][i]*MeanXRe[nuser][i] );
        }
    
    for( i=0; i<(_ChipLen/2); i++ )
    {
        Covr[i][0][0] = 0;
        Covr[i][0][1] = 0;
        Covr[i][1][0] = 0;
        Covr[i][1][1] = 0;
    }
    for( nuser=0; nuser<_NUser; nuser++ ) 
        for( i=0; i<(_ChipLen/2); i++ )
        {
            CovX[nuser][i][0][0] = 1;
            CovX[nuser][i][0][1] = 0;
            CovX[nuser][i][1][0] = 0;
            CovX[nuser][i][1][1] = 1;
        }
    for( i=0; i<(_ChipLen/2); i++ ) 
        for( nuser=0; nuser<_NUser; nuser++ )
        {
            Temp[nuser][i][0][0] = R[nuser][i][0][0]*CovX[nuser][i][0][0] + R[nuser][i][0][1]*CovX[nuser][i][1][0];
            Temp[nuser][i][0][1] = R[nuser][i][0][0]*CovX[nuser][i][0][1] + R[nuser][i][0][1]*CovX[nuser][i][1][1];
            Temp[nuser][i][1][0] = R[nuser][i][1][0]*CovX[nuser][i][0][0] + R[nuser][i][1][1]*CovX[nuser][i][1][0];
            Temp[nuser][i][1][1] = R[nuser][i][1][0]*CovX[nuser][i][0][1] + R[nuser][i][1][1]*CovX[nuser][i][1][1];

            Covr[i][0][0] += Temp[nuser][i][0][0]*RT[nuser][i][0][0] + Temp[nuser][i][0][1]*RT[nuser][i][1][0];
            Covr[i][0][1] += Temp[nuser][i][0][0]*RT[nuser][i][0][1] + Temp[nuser][i][0][1]*RT[nuser][i][1][1];
            Covr[i][1][0] += Temp[nuser][i][1][0]*RT[nuser][i][0][0] + Temp[nuser][i][1][1]*RT[nuser][i][1][0];
            Covr[i][1][1] += Temp[nuser][i][1][0]*RT[nuser][i][0][1] + Temp[nuser][i][1][1]*RT[nuser][i][1][1];
        }
    for( i=0; i<(_ChipLen/2); i++ )
    {
        Covr[i][0][0] += s2;
        Covr[i][1][1] += s2;
    }
        
    for( it=0; it<_ItNum; it++ ) 
        for( nuser=0; nuser<_NUser; nuser++ )
        {
            // Produce the LLR values for the de-interleaved chip sequence.
            for( i=0; i<(_ChipLen/2); i++ )    
            {
                //(39/u)
                MeanLRe[nuser][i] = MeanrRe[i] - HRe[nuser][i]*MeanXRe[nuser][i] + HImg[nuser][i]*MeanXImg[nuser][i];
                MeanLImg[nuser][i] = MeanrImg[i] - HRe[nuser][i]*MeanXImg[nuser][i] - HImg[nuser][i]*MeanXRe[nuser][i];
                //(38/u)
                MeanLHRe[nuser][i] = HRe[nuser][i]*MeanLRe[nuser][i] + HImg[nuser][i]*MeanLImg[nuser][i];
                MeanLHImg[nuser][i] = HRe[nuser][i]*MeanLImg[nuser][i] - HImg[nuser][i]*MeanLRe[nuser][i];
                //(35)
                rkHRe[nuser][i] = HRe[nuser][i]*FFTrRe[i] + HImg[nuser][i]*FFTrImg[i];
                rkHImg[nuser][i] = HRe[nuser][i]*FFTrImg[i] - HImg[nuser][i]*FFTrRe[i];
                //(39/d/1)
                Temp[nuser][i][0][0] = R[nuser][i][0][0]*CovX[nuser][i][0][0] + R[nuser][i][0][1]*CovX[nuser][i][1][0];
                Temp[nuser][i][0][1] = R[nuser][i][0][0]*CovX[nuser][i][0][1] + R[nuser][i][0][1]*CovX[nuser][i][1][1];
                Temp[nuser][i][1][0] = R[nuser][i][1][0]*CovX[nuser][i][0][0] + R[nuser][i][1][1]*CovX[nuser][i][1][0];
                Temp[nuser][i][1][1] = R[nuser][i][1][0]*CovX[nuser][i][0][1] + R[nuser][i][1][1]*CovX[nuser][i][1][1];
                //(39/d/2)
                CovL[nuser][i][0][0] = Covr[i][0][0] - Temp[nuser][i][0][0]*RT[nuser][i][0][0] - Temp[nuser][i][0][1]*RT[nuser][i][1][0];
                CovL[nuser][i][0][1] = Covr[i][0][1] - Temp[nuser][i][0][0]*RT[nuser][i][0][1] - Temp[nuser][i][0][1]*RT[nuser][i][1][1];
                CovL[nuser][i][1][0] = Covr[i][1][0] - Temp[nuser][i][1][0]*RT[nuser][i][0][0] - Temp[nuser][i][1][1]*RT[nuser][i][1][0];
                CovL[nuser][i][1][1] = Covr[i][1][1] - Temp[nuser][i][1][0]*RT[nuser][i][0][1] - Temp[nuser][i][1][1]*RT[nuser][i][1][1];
                //(38/d/1)
                Temp[nuser][i][0][0] = RT[nuser][i][0][0]*CovL[nuser][i][0][0] + RT[nuser][i][0][1]*CovL[nuser][i][1][0];
                Temp[nuser][i][0][1] = RT[nuser][i][0][0]*CovL[nuser][i][0][1] + RT[nuser][i][0][1]*CovL[nuser][i][1][1];
                Temp[nuser][i][1][0] = RT[nuser][i][1][0]*CovL[nuser][i][0][0] + RT[nuser][i][1][1]*CovL[nuser][i][1][0];
                Temp[nuser][i][1][1] = RT[nuser][i][1][0]*CovL[nuser][i][0][1] + RT[nuser][i][1][1]*CovL[nuser][i][1][1];
                //(38/d/2)
                CovLH[nuser][i][0][0] = Temp[nuser][i][0][0]*R[nuser][i][0][0] + Temp[nuser][i][0][1]*R[nuser][i][1][0];
                CovLH[nuser][i][0][1] = Temp[nuser][i][0][0]*R[nuser][i][0][1] + Temp[nuser][i][0][1]*R[nuser][i][1][1];
                CovLH[nuser][i][1][0] = Temp[nuser][i][1][0]*R[nuser][i][0][0] + Temp[nuser][i][1][1]*R[nuser][i][1][0];
                CovLH[nuser][i][1][1] = Temp[nuser][i][1][0]*R[nuser][i][0][1] + Temp[nuser][i][1][1]*R[nuser][i][1][1];
                //(37)
                ChipRe[nuser][i] = 2 * (pow(HRe[nuser][i],2)+pow(HImg[nuser][i],2)) * ( rkHRe[nuser][i] - MeanLHRe[nuser][i] ) / CovLH[nuser][i][0][0];
                ChipImg[nuser][i] = 2 * (pow(HRe[nuser][i],2)+pow(HImg[nuser][i],2)) * ( rkHImg[nuser][i] - MeanLHImg[nuser][i] ) / CovLH[nuser][i][1][1];
            }

        for ( i=j=m=0; i<_ChipLen; i++ )
        {
            if ( i % 2 == 0 ) TempChip[nuser][i] = ChipRe[nuser][j++];
            else if ( i % 2 == 1 ) TempChip[nuser][i] = ChipImg[nuser][m++];
        }

        for ( i=0; i<_ChipLen; i++ ) 
        {
            Chip[nuser][ScrambleRule[nuser][i]] = TempChip[nuser][i];

            if ( Chip[nuser][ScrambleRule[nuser][i]] > 20 ) Chip[nuser][ScrambleRule[nuser][i]] = 20;
            if ( Chip[nuser][ScrambleRule[nuser][i]] < -20 ) Chip[nuser][ScrambleRule[nuser][i]] = -20;
        }
                
        // De-spreading operation.        
        for( m=i=0; i<_DataLen; i++ )
        {
            AppLlr[nuser][i] = SpreadSeq[0] * Chip[nuser][m++];
            for( j=1; j<_SpreadLen; j++ ) AppLlr[nuser][i] += SpreadSeq[j] * Chip[nuser][m++];
        }

        // Data Reinterpritation
        for( j=0; j<_NUser; j++ ) 
          for( i=0; i<_DataLen; i++ ) 
             AppLlrOut[(j * _DataLen) + i] = AppLlr[j][i];
        
        // Feed the AppLlr to decoder, if there is a FEC codec in the system.

        // Spreading operation: Produce the extrinsic LLR for each chip
        for( m=i=0; i<_DataLen; i++ ) 
            for( j=0; j<_SpreadLen; j++, m++ )
                Ext[nuser][m] = SpreadSeq[j] * AppLlr[nuser][i] - Chip[nuser][m];

        // Update the statistic variable together with the interleaving process.
        for( i=0; i<(_ChipLen/2); i++ )
        {
            MeanXRe[nuser][i] = tanh( Ext[nuser][ScrambleRule[nuser][i*2]] / 2 );
            MeanXImg[nuser][i] = tanh( Ext[nuser][ScrambleRule[nuser][i*2+1]] / 2 );
            CovX[nuser][i][0][0] = 1.0 - pow(MeanXRe[nuser][i],2);
            CovX[nuser][i][1][1] = 1.0 - pow(MeanXImg[nuser][i],2);
        }

        for( i=0; i<(_ChipLen/2); i++ )
        {
            MeanrRe[i] = 0;
            MeanrImg[i] = 0;
            for( j=0; j<_NUser; j++ )
            {
                MeanrRe[i] += HRe[j][i]*MeanXRe[j][i] - HImg[j][i]*MeanXImg[j][i];
                MeanrImg[i] += HRe[j][i]*MeanXImg[j][i] + HImg[j][i]*MeanXRe[j][i];
            }

            Covr[i][0][0] = sigma*sigma;
            Covr[i][0][1] = 0;
            Covr[i][1][0] = 0;
            Covr[i][1][1] = sigma*sigma;
            for( j=0; j<_NUser; j++ )
            {
                Temp[j][i][0][0] = R[j][i][0][0]*CovX[j][i][0][0] + R[j][i][0][1]*CovX[j][i][1][0];
                Temp[j][i][0][1] = R[j][i][0][0]*CovX[j][i][0][1] + R[j][i][0][1]*CovX[j][i][1][1];
                Temp[j][i][1][0] = R[j][i][1][0]*CovX[j][i][0][0] + R[j][i][1][1]*CovX[j][i][1][0];
                Temp[j][i][1][1] = R[j][i][1][0]*CovX[j][i][0][1] + R[j][i][1][1]*CovX[j][i][1][1];

                Covr[i][0][0] += Temp[j][i][0][0]*RT[j][i][0][0] + Temp[j][i][0][1]*RT[j][i][1][0];
                Covr[i][0][1] += Temp[j][i][0][0]*RT[j][i][0][1] + Temp[j][i][0][1]*RT[j][i][1][1];
                Covr[i][1][0] += Temp[j][i][1][0]*RT[j][i][0][0] + Temp[j][i][1][1]*RT[j][i][1][0];
                Covr[i][1][1] += Temp[j][i][1][0]*RT[j][i][0][1] + Temp[j][i][1][1]*RT[j][i][1][1];
            }
        }
    }
} 

/* main function that interfaces with MATLAB */
void mexFunction(
                 int            nlhs,
                 mxArray       *plhs[],
                 int            nrhs,
                 const mxArray *prhs[] )
{
    double *rre_double;
    double *rimg_double;
    double *rayre_double;
    double *rayimg_double;    
    double *spreadseq_double;
    double *scramblerule_double;
    double *appllrout_double;        

    // Получаем информационные и служебные данные 
    rre_double = mxGetPr(__RRE);
    double rre [CONVLEN(PREFIX)];
    for (int i = 0; i < CONVLEN(PREFIX); i++) rre[i] = rre_double[i]; 
    
    rimg_double = mxGetPr(__RIMG);
    double rimg [CONVLEN(PREFIX)];
    for (int i = 0; i < CONVLEN(PREFIX); i++) rimg[i] = rimg_double[i];

    rayre_double = mxGetPr(__RAYRE);
    double rayre [DENULLIFY(PREFIX)];
    for (int i = 0; i < DENULLIFY(PREFIX); i++) rayre[i] = rayre_double[i]; 

    rayimg_double = mxGetPr(__RAYIMG);
    double rayimg [DENULLIFY(PREFIX)];
    for (int i = 0; i < DENULLIFY(PREFIX); i++) rayimg[i] = rayimg_double[i];

    spreadseq_double = mxGetPr(__SPREADSEQ);
    double spreadseq [SPREADLEN];
    for (int i = 0; i < SPREADLEN; i++) spreadseq[i] = spreadseq_double[i];

    scramblerule_double = mxGetPr(__SCRAMBLERULEIN);
    int scramblerule [FULLCODDATA];
    for (int i = 0; i < FULLCODDATA; i++) scramblerule[i] = (int) scramblerule_double[i];

    double sigma     = (double) *mxGetPr(__SIGMA); 
    int    datalen   = (int)    *mxGetPr(__DATALEN);
    int    chiplen   = (int)    *mxGetPr(__CHIPLEN);
    int    prefix    = (int)    *mxGetPr(__PREFIX);
    int    nuser     = (int)    *mxGetPr(__NUSER);
    int    itnum     = (int)    *mxGetPr(__ITNUM);
    int    spreadlen = (int)    *mxGetPr(__SPREADLEN);
    
    double appllrout [FULLDATA];

    Receiver( rre, rimg, sigma, datalen, chiplen, prefix, nuser, itnum, rayre, rayimg, spreadlen, spreadseq, scramblerule, appllrout);      
    
    /* Передаем в Матлаб данные */
    __APPLLROUT = mxCreateDoubleMatrix(FULLDATA, 1, mxREAL);
    appllrout_double = mxGetPr(__APPLLROUT);
    
    for ( int i = 0; i < FULLDATA; i++ ) appllrout_double[i] = appllrout[i];    

    return;
}

 

Этот код я хочу понять и правильно описать своими словами. Это .cpp файл в проекте, по которому создается mex файл для Matlab. В коде на вход подаются комплексные данные, представляющие собой сумму данных N-пользователей.

Например, если передавались 32 бита для каждого из 2х пользователей, то в итоге надо передать 64 бита, с ПСП это 64*16 = 1024 бит, после модуляции у нас 512 выборок, так как 2 пользователя, мы сложили 256+256 и на приемник поступает 256 выборок + длина ЦП, то есть 266, например.

Вопрос такой, там в коде есть ковариационные матрицы, итеративный метод и расчет максимума функции правдоподобия при разделении. Вам знаком этот подход? О нем можно где-то почитать?

Как бы Вы в двух словах могли это описать, если задача понятна? Хочу понять физику явления, зачем это все и как максимум правдоподобия на практике помогает разделить суперпозицию сигналов.

 

Чтобы было проще, сейчас я могу объяснить примерно так. Мы можем разделить один сигнал на несколько, используя один из 9 методов (9+), куда входят сингулярный анализ, метод независимых компонент (ICA), метод главных компонент. Однако так я подхожу к вопросу издалека. Эти методы разработаны для разделения звуковых сигналов, что используется, например, в Dolby Digital, когда нужно выделить разные музыкальные инструменты из дорожки. В каких-то из методов используются итерации, но зачем? Я еще дохожу. В ICA находится ковариационная зависимость, а разделить сигналы можно только тогда, когда есть избыток информации, классика - это у нас несколько микрофонов или приемных антенн.

Я начал разбираться и понял (еще понимаю), что в CPP, который прислал, там избыток за счет того, что мы знаем импульсную характеристику канала связи (эквалайзер с нулевым усилием (ZF)). И фактически здесь для разделения используется что-то подобное, а потом та же ковариационная матрица, расчет правдоподобия и от итерации к итерации мы находим максимум функции правдоподобия, что позволяет нам разделить сигналы. По коду я еще вижу, что сначала данные разделяются, а потом уже декодируются, а ПСП, дескать, мы знаем, и они подаются на вход Приемника.

appllr = ReceiverMx(rreout, rimgout, sigma, datalen, chiplen, prefix, nuser, itnum, rayre, rayimg, spreadlen, spreadseq, scramblerule);

Это строчка из Матлаб, scramblerule и spreadseq говорят сами за себя. И после разделения и избавления от ПСП (получаем данные для N-пользователей) мы заново на конце итерации оцениваем статистические параметры, при чем для этого используются все те же ковариационные матрицы...

 

Но в последних нескольких предложениях я еще не понял, что сказал, можете навести на мысль? Все таки как объяснить,

 

1. зачем мы находим среднее значение Mean,

2. зачем нужны ковариационные матрицы,

3. избыток какой информации в данном коде? Я прав про то, что вся фишка в том, что мы знаем комплексную импульсную характеристику канала связи H?

4. как максимум функции правдоподобия позволяет нам один сигнал разделить на 2 и более, да еще и так в итоге, что мы находим данные для каждого пользователя?

5. зачем мы комплексные коэффициенты импульсной характеристики RayRe+1i*RayImg приводим в таком гармоническом виде? (ChipLen - длина чипа, для вышеописанной частной ситуации это 1024 бита/ 2 пользователя = 512. Или же 32*16 (длина ПСП) = 512):

for( i=0; i<_ChipLen/2; i++ )
        {
            if ( _Prefix == 0 )
            {
                j=0;
                _HRe[i] += RayRe[0]*cos(6.283185307*i*j/_ChipLen*2) + RayImg[0]*sin(6.283185307*i*j/_ChipLen*2);
                _HImg[i] += -RayRe[0]*sin(6.283185307*i*j/_ChipLen*2) + RayImg[0]*cos(6.283185307*i*j/_ChipLen*2);
            }

            if ( _Prefix > 0 )
            {
                for( j=0; j<DENULLIFY(_Prefix); j++ )
                {
                    _HRe[i] += RayRe[j]*cos(6.283185307*i*j/_ChipLen*2) + RayImg[j]*sin(6.283185307*i*j/_ChipLen*2);
                    _HImg[i] += -RayRe[j]*sin(6.283185307*i*j/_ChipLen*2) + RayImg[j]*cos(6.283185307*i*j/_ChipLen*2);
                }
            }           
        }

Share this post


Link to post
Share on other sites

Если это сильно сложно, давайте по порядку теперь, общий вопрос я изложил, а теперь на понимание.

 

1. Зачем при разделении двух и более сигналов, сложенных друг с другом, используется ковариация?

2. Как участвует функция правдоподобия при разделении сигналов? Я понимаю, что функция правдоподобия позволяет понять, насколько правдоподобно, что апостериорная статистическая информация близка к тем или иным априорным распределением. Но при чем тут может быть разделение сигналов?

Share this post


Link to post
Share on other sites

Я уже почитал много чего, поэтому в точку! Но надо подробнее.

 

Подробнее хочу понять теперь следующее:

 

После нахождения e_ese = Chip:

 

for ( i=0; i<_ChipLen; i++ ) 
        {
            Chip[nuser][ScrambleRule[nuser][i]] = TempChip[nuser][i];

            if ( Chip[nuser][ScrambleRule[nuser][i]] > 20 ) Chip[nuser][ScrambleRule[nuser][i]] = 20;
            if ( Chip[nuser][ScrambleRule[nuser][i]] < -20 ) Chip[nuser][ScrambleRule[nuser][i]] = -20;
        }

 

происходит удаление ПСП, а затем снова операция расширения спектра с помощью ПСП. e_dec = Ext показывает функцию правдоподобия как отношения вероятностей того, что бит 1 и 0. Но так то e_dec должна быть равна нулю, тут же в коде взаимообратные операции Spreading/Despreading. Значит вся магия дейсствительно в ограничителях?

 

if ( Chip[nuser][ScrambleRule[nuser][i]] > 20 ) Chip[nuser][ScrambleRule[nuser][i]] = 20;
            if ( Chip[nuser][ScrambleRule[nuser][i]] < -20 ) Chip[nuser][ScrambleRule[nuser][i]] = -20;

Как лучше это понимать?

 

Кстати, а в CDMA есть возможность 2+ раза разделять пользователей (в 2+ этапа). Допустим, сначала из потока выделил 4 пользователя, а затем из каждого еще по 4 саб-пользователя? Есть ли в этом смысл, существует ли такая техника и если да, то как она называется?

Share this post


Link to post
Share on other sites
Есть ли в этом смысл, существует ли такая техника и если да, то как она называется?

В перенасыщенных CDMA делают нечто подобное, чтобы уйти от МП-приемника, который требует огромных вычислительных затрат. Сигналы неортогональные. В случае классической CDMA в этом нет смысла, потому что нет ничего проще вычисления корреляций.

Share this post


Link to post
Share on other sites

Ок. Еще упорно и нундо не могу понять одну строчку:

Ext[nuser][m] = SpreadSeq[j] * AppLlr[nuser][i] - Chip[nuser][m];

 

Я смоделировал отдельно без шумов, Ext получается 0!!!! Так как Chip и Spreaded AppLlr это одинаковые массивы. Где чего я не понял? Так то по логике турбоэквалайзера + CDMA на каждой итерации Chip должен использоваться для обновления матожидания и дисперсии, но эта строка это как-то странно.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Sign in to follow this