Перейти к содержанию

    

Разделение суммы 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);
                }
            }           
        }

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

 

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

 

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

 

После нахождения 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 саб-пользователя? Есть ли в этом смысл, существует ли такая техника и если да, то как она называется?

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты
Есть ли в этом смысл, существует ли такая техника и если да, то как она называется?

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

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

 

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Для публикации сообщений создайте учётную запись или авторизуйтесь

Вы должны быть пользователем, чтобы оставить комментарий

Создать учетную запись

Зарегистрируйте новую учётную запись в нашем сообществе. Это очень просто!

Регистрация нового пользователя

Войти

Уже есть аккаунт? Войти в систему.

Войти
Авторизация