discuss-gnuradio
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Discuss-gnuradio] Question on gr-trellis PS and PI matricies in FSm


From: Tim Meehan
Subject: [Discuss-gnuradio] Question on gr-trellis PS and PI matricies in FSm
Date: Sun, 4 Feb 2007 06:54:36 -0800

Achilleas,

I was able to get your modifications to gr-trellis working in my
original application.  Thanks again for the changes.

It took me a bit to figure out how to include bit-stuffing and
unstuffing blocks into the gnuradio framework.  Each block has an
input rate that is different from the output rate, and is not
constant.  I am still not convinced that I have it right but it seems
to work.

I will take a look at your SISO examples when I get some time.

Tim

On 2/1/07, Achilleas Anastasopoulos <address@hidden> wrote:
Tim,

the BCJR algorithm is already implemented in gr-trellis, under the name
SISO (stands for soft in soft out).
There are two versions of it, similar to the two version of the Viterbi
algorithm...
I even have some examples on turbo equalization, and serially
concatenated convolutional codes in gnuradio-examples/python/channel-coding.

Let me know if you have any problems running and understanding
those apps.

Achilleas

Tim Meehan wrote:
> Achilleas,
>
> Thanks for the quick change.  I missed the "fsm.i" file the first time
> around, and it took me a while to figure out what was wrong.  After I
> put the correct fsm.i file in everything compiled, and it seems to
> work.
>
> I'll get back to my original application in the morning and will send
> you a quick note letting you know how it goes.
>
> Do you know if there is any interest for a BCJR algorithm?  I did I
> quick search on the discuss-gnuradio archive for "BCJR" and did not
> recieve any hits.
>
> Thanks again for all of your time
>
> Tim
>
>
>
>
> On 2/1/07, Achilleas Anastasopoulos <address@hidden> wrote:
>
>> Dear Tim,
>>
>> it was indeed a minor modification that was required.
>> I made and it is working fine now.
>>
>> It will take me some time to clean up the code
>> (as you can see the TM() method complains, but
>> the correction should be trivial there...)
>>
>>
>>
>> If you want to start working on your project
>> just copy the files I attach
>> to the gr-trellis/src/lib directory
>> execute ./generat_all.py
>> and then execute make and make install from gr-trellis
>> directory.
>>
>> I am also attaching a toy irregular FSM file.
>> put it in the directory
>> gnuradio-examples/python/channel-coding/fms_files.
>> Then run the modified test_tcm.py and run it.
>> It works.
>>
>> Please let me know if you catch any bugs that I missed...
>>
>> Achilleas
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> ============================
>> Dear Tim,,
>>
>> indeed the current implementation supports only this kind of "regular"
>> trellises.
>> I believe the intruduction of non-regular trellises (ie, ones having
>> exactly I outgoing edges from each state, but different incoming edges
>> to each state) is a good idea.
>>
>> They way I think about it, I have to substitute the PS (and PI)
>> internally generated matrix which is now of size SxI with a vector of
>> size S where each element is a  vector of size equal to the number of
>> incoming edges in the particular state. Then a couple of minor
>> modifications to the VA, and it is done...
>>
>> I am not sure how familiar you are with the code I developed in
>> gr-trellis, but I am willing to help you make the modifications.
>>
>> Best
>> Achilleas
>>
>>
>> > Message: 7
>> > Date: Thu, 1 Feb 2007 10:31:56 -0500
>> > From: "Tim Meehan" <address@hidden>
>> > Subject: [Discuss-gnuradio] Question on gr-trellis PS and PI matricies
>> >       in FSM
>> > To: address@hidden
>> > Message-ID:
>> >       <address@hidden>
>> > Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>> >
>> > Hello all,
>> >
>> > I was trying to use the Viterbi decoder in gr-trellis for maximum
>> > likelihood sequence detection of a bit-stuffed data source.  Here the
>> > FSM is the bit-stuffing operation.  A reference for this is
>> >
>> > Meehan, T.; Hicks, J.; Kragh, F., "Maximum Likelihood Sequence
>> > Detection of a Bit-stuffed Data Source," Communications, 2006. ICC
>> > '06. IEEE International Conference on , vol.4, no.pp.1514-1519, June
>> > 2006
>> >
>> > For this application the number of branches entering a state is not
>> > the same as the number leaving at a state.  The documentation for
>> > gr-trellis states
>> >
>> > "In the following we assume that for any state, the number of incoming
>> > transitions is the same as the number of outgoing transitions, ie,
>> > equal to I. All applications of interest have FSMs satisfying this
>> > requirement."
>> >
>> > Before I started digging in to modifying the code (or porting my
>> > existing decoder to the gnuradio framework) I wanted to see if anyone
>> > else was working to support different number of input branches to
>> > output branches.
>> >
>> > I think at a minimum the method "fsm::generate_PS_PI()" should check
>> > to verify that the FSM description does in fact have the same number
>> > of incoming branches and outing branches at each state.
>> >
>> > If there is interest in making gr-trellis support this type of FSM I
>> > can work on it, else I will probably just port my existing work over.
>> >
>> > Please let me know what you think.
>> >
>> > P.S.  Forgive me if I this shows up twice.  I did not see my original
>> > post and have changed my address with the list.
>> >
>> > Tim
>>
>>
>>
>> /* -*- c++ -*- */
>> /*
>>  * Copyright 2004 Free Software Foundation, Inc.
>>  *
>>  * This file is part of GNU Radio
>>  *
>>  * GNU Radio is free software; you can redistribute it and/or modify
>>  * it under the terms of the GNU General Public License as published by
>>  * the Free Software Foundation; either version 2, or (at your option)
>>  * any later version.
>>  *
>>  * GNU Radio is distributed in the hope that it will be useful,
>>  * but WITHOUT ANY WARRANTY; without even the implied warranty of
>>  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
>>  * GNU General Public License for more details.
>>  *
>>  * You should have received a copy of the GNU General Public License
>>  * along with GNU Radio; see the file COPYING.  If not, write to
>>  * the Free Software Foundation, Inc., 51 Franklin Street,
>>  * Boston, MA 02110-1301, USA.
>>  */
>>
>> #ifdef HAVE_CONFIG_H
>> #include "config.h"
>> #endif
>>
>> #include <trellis_siso_combined_f.h>
>> #include <gr_io_signature.h>
>> #include <stdexcept>
>> #include <assert.h>
>> #include <iostream>
>>
>> static const float INF = 1.0e9;
>>
>> trellis_siso_combined_f_sptr
>> trellis_make_siso_combined_f (
>>     const fsm &FSM,
>>     int K,
>>     int S0,
>>     int SK,
>>     bool POSTI,
>>     bool POSTO,
>>     trellis_siso_type_t SISO_TYPE,
>>     int D,
>>     const std::vector<float> &TABLE,
>>     trellis_metric_type_t TYPE)
>> {
>>   return trellis_siso_combined_f_sptr (new trellis_siso_combined_f
>> (FSM,K,S0,SK,POSTI,POSTO,SISO_TYPE,D,TABLE,TYPE));
>> }
>>
>> trellis_siso_combined_f::trellis_siso_combined_f (
>>     const fsm &FSM,
>>     int K,
>>     int S0,
>>     int SK,
>>     bool POSTI,
>>     bool POSTO,
>>     trellis_siso_type_t SISO_TYPE,
>>     int D,
>>     const std::vector<float> &TABLE,
>>     trellis_metric_type_t TYPE)
>>   : gr_block ("siso_combined_f",
>>                           gr_make_io_signature (1, -1, sizeof (float)),
>>                           gr_make_io_signature (1, -1, sizeof (float))),
>>   d_FSM (FSM),
>>   d_K (K),
>>   d_S0 (S0),
>>   d_SK (SK),
>>   d_POSTI (POSTI),
>>   d_POSTO (POSTO),
>>   d_SISO_TYPE (SISO_TYPE),
>>   d_D (D),
>>   d_TABLE (TABLE),
>>   d_TYPE (TYPE)//,
>>   //d_alpha(FSM.S()*(K+1)),
>>   //d_beta(FSM.S()*(K+1))
>> {
>>     int multiple;
>>     if (d_POSTI && d_POSTO)
>>         multiple = d_FSM.I()+d_FSM.O();
>>     else if(d_POSTI)
>>         multiple = d_FSM.I();
>>     else if(d_POSTO)
>>         multiple = d_FSM.O();
>>     else
>>         throw std::runtime_error ("Not both POSTI and POSTO can be
>> false.");
>>     //printf("constructor: Multiple = %d\n",multiple);
>>     set_output_multiple (d_K*multiple);
>>     //what is the meaning of relative rate for a block with 2 inputs?
>>     //set_relative_rate ( multiple / ((double) d_FSM.I()) );
>>     // it turns out that the above gives problems in the scheduler, so
>>     // let's try (assumption O>I)
>>     //set_relative_rate ( multiple / ((double) d_FSM.O()) );
>>     // I am tempted to automate like this
>>     if(d_FSM.I() <= d_D)
>>       set_relative_rate ( multiple / ((double) d_D) );
>>     else
>>       set_relative_rate ( multiple / ((double) d_FSM.I()) );
>> }
>>
>>
>> void
>> trellis_siso_combined_f::forecast (int noutput_items, gr_vector_int
>> &ninput_items_required)
>> {
>>   int multiple;
>>   if (d_POSTI && d_POSTO)
>>       multiple = d_FSM.I()+d_FSM.O();
>>   else if(d_POSTI)
>>       multiple = d_FSM.I();
>>   else if(d_POSTO)
>>       multiple = d_FSM.O();
>>   else
>>       throw std::runtime_error ("Not both POSTI and POSTO can be
>> false.");
>>   //printf("forecast: Multiple = %d\n",multiple);
>>   assert (noutput_items % (d_K*multiple) == 0);
>>   int input_required1 =  d_FSM.I() * (noutput_items/multiple) ;
>>   int input_required2 =  d_D * (noutput_items/multiple) ;
>>   //printf("forecast: Output requirements:  %d\n",noutput_items);
>>   //printf("forecast: Input requirements:  %d
>> %d\n",input_required1,input_required2);
>>   unsigned ninputs = ninput_items_required.size();
>>   assert(ninputs % 2 == 0);
>>   for (unsigned int i = 0; i < ninputs/2; i++) {
>>     ninput_items_required[2*i] = input_required1;
>>     ninput_items_required[2*i+1] = input_required2;
>>   }
>> }
>>
>> inline float min(float a, float b)
>> {
>>   return a <= b ? a : b;
>> }
>>
>> inline float min_star(float a, float b)
>> {
>>   return (a <= b ? a : b)-log(1+exp(a <= b ? a-b : b-a));
>> }
>>
>> void siso_algorithm_combined(int I, int S, int O,
>>              const std::vector<int> &NS,
>>              const std::vector<int> &OS,
>>              const std::vector< std::vector<int> > &PS,
>>              const std::vector< std::vector<int> > &PI,
>>              int K,
>>              int S0,int SK,
>>              bool POSTI, bool POSTO,
>>              float (*p2mymin)(float,float),
>>              int D,
>>              const std::vector<float> &TABLE,
>>              trellis_metric_type_t TYPE,
>>              const float *priori, const float *observations, float
>> *post//,
>>              //std::vector<float> &alpha,
>>              //std::vector<float> &beta
>>              )
>> {
>>   float norm,mm,minm;
>>   std::vector<float> alpha(S*(K+1));
>>   std::vector<float> beta(S*(K+1));
>>   float *prioro = new float[O*K];
>>
>>
>>   if(S0<0) { // initial state not specified
>>       for(int i=0;i<S;i++) alpha[0*S+i]=0;
>>   }
>>   else {
>>       for(int i=0;i<S;i++) alpha[0*S+i]=INF;
>>       alpha[0*S+S0]=0.0;
>>   }
>>
>>   for(int k=0;k<K;k++) { // forward recursion
>>       calc_metric(O, D, TABLE, &(observations[k*D]),
>> &(prioro[k*O]),TYPE); // calc metrics
>>       norm=INF;
>>       for(int j=0;j<S;j++) {
>>           minm=INF;
>>           for(int i=0;i<PS[j].size();i++) {
>>               int i0 = j*I+i;
>>
>> 
mm=alpha[k*S+PS[j][i]]+priori[k*I+PI[j][i]]+prioro[k*O+OS[PS[j][i]*I+PI[j][i]]];
>>
>>               minm=(*p2mymin)(minm,mm);
>>           }
>>           alpha[(k+1)*S+j]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int j=0;j<S;j++)
>>           alpha[(k+1)*S+j]-=norm; // normalize total metrics so they
>> do not explode
>>   }
>>
>>   if(SK<0) { // final state not specified
>>       for(int i=0;i<S;i++) beta[K*S+i]=0;
>>   }
>>   else {
>>       for(int i=0;i<S;i++) beta[K*S+i]=INF;
>>       beta[K*S+SK]=0.0;
>>   }
>>
>>   for(int k=K-1;k>=0;k--) { // backward recursion
>>       norm=INF;
>>       for(int j=0;j<S;j++) {
>>           minm=INF;
>>           for(int i=0;i<I;i++) {
>>               int i0 = j*I+i;
>>               mm=beta[(k+1)*S+NS[i0]]+priori[k*I+i]+prioro[k*O+OS[i0]];
>>               minm=(*p2mymin)(minm,mm);
>>           }
>>           beta[k*S+j]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int j=0;j<S;j++)
>>           beta[k*S+j]-=norm; // normalize total metrics so they do not
>> explode
>>   }
>>
>>
>>   if (POSTI && POSTO)
>>   {
>>     for(int k=0;k<K;k++) { // input combining
>>         norm=INF;
>>         for(int i=0;i<I;i++) {
>>             minm=INF;
>>             for(int j=0;j<S;j++) {
>>
>> mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
>>                 minm=(*p2mymin)(minm,mm);
>>             }
>>             post[k*(I+O)+i]=minm;
>>             if(minm<norm) norm=minm;
>>         }
>>         for(int i=0;i<I;i++)
>>             post[k*(I+O)+i]-=norm; // normalize metrics
>>     }
>>
>>
>>     for(int k=0;k<K;k++) { // output combining
>>         norm=INF;
>>         for(int n=0;n<O;n++) {
>>             minm=INF;
>>             for(int j=0;j<S;j++) {
>>                 for(int i=0;i<I;i++) {
>>                     mm= (n==OS[j*I+i] ?
>> alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
>>                     minm=(*p2mymin)(minm,mm);
>>                 }
>>             }
>>             post[k*(I+O)+I+n]=minm;
>>             if(minm<norm) norm=minm;
>>         }
>>         for(int n=0;n<O;n++)
>>             post[k*(I+O)+I+n]-=norm; // normalize metrics
>>     }
>>   }
>>   else if(POSTI)
>>   {
>>     for(int k=0;k<K;k++) { // input combining
>>         norm=INF;
>>         for(int i=0;i<I;i++) {
>>             minm=INF;
>>             for(int j=0;j<S;j++) {
>>
>> mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
>>                 minm=(*p2mymin)(minm,mm);
>>             }
>>             post[k*I+i]=minm;
>>             if(minm<norm) norm=minm;
>>         }
>>         for(int i=0;i<I;i++)
>>             post[k*I+i]-=norm; // normalize metrics
>>     }
>>   }
>>   else if(POSTO)
>>   {
>>     for(int k=0;k<K;k++) { // output combining
>>         norm=INF;
>>         for(int n=0;n<O;n++) {
>>             minm=INF;
>>             for(int j=0;j<S;j++) {
>>                 for(int i=0;i<I;i++) {
>>                     mm= (n==OS[j*I+i] ?
>> alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
>>                     minm=(*p2mymin)(minm,mm);
>>                 }
>>             }
>>             post[k*O+n]=minm;
>>             if(minm<norm) norm=minm;
>>         }
>>         for(int n=0;n<O;n++)
>>             post[k*O+n]-=norm; // normalize metrics
>>     }
>>   }
>>   else
>>     throw std::runtime_error ("Not both POSTI and POSTO can be false.");
>>
>>   delete [] prioro;
>>
>> }
>>
>>
>>
>>
>>
>>
>> int
>> trellis_siso_combined_f::general_work (int noutput_items,
>>                         gr_vector_int &ninput_items,
>>                         gr_vector_const_void_star &input_items,
>>                         gr_vector_void_star &output_items)
>> {
>>   assert (input_items.size() == 2*output_items.size());
>>   int nstreams = output_items.size();
>>   //printf("general_work:Streams:  %d\n",nstreams);
>>   int multiple;
>>   if (d_POSTI && d_POSTO)
>>       multiple = d_FSM.I()+d_FSM.O();
>>   else if(d_POSTI)
>>       multiple = d_FSM.I();
>>   else if(d_POSTO)
>>       multiple = d_FSM.O();
>>   else
>>       throw std::runtime_error ("Not both POSTI and POSTO can be
>> false.");
>>
>>   assert (noutput_items % (d_K*multiple) == 0);
>>   int nblocks = noutput_items / (d_K*multiple);
>>   //printf("general_work:Blocks:  %d\n",nblocks);
>>   //for(int i=0;i<ninput_items.size();i++)
>>       //printf("general_work:Input items available:
>> %d\n",ninput_items[i]);
>>
>>   float (*p2min)(float, float) = NULL;
>>   if(d_SISO_TYPE == TRELLIS_MIN_SUM)
>>     p2min = &min;
>>   else if(d_SISO_TYPE == TRELLIS_SUM_PRODUCT)
>>     p2min = &min_star;
>>
>>
>>   for (int m=0;m<nstreams;m++) {
>>     const float *in1 = (const float *) input_items[2*m];
>>     const float *in2 = (const float *) input_items[2*m+1];
>>     float *out = (float *) output_items[m];
>>     for (int n=0;n<nblocks;n++) {
>>       siso_algorithm_combined(d_FSM.I(),d_FSM.S(),d_FSM.O(),
>>         d_FSM.NS(),d_FSM.OS(),d_FSM.PS(),d_FSM.PI(),
>>         d_K,d_S0,d_SK,
>>         d_POSTI,d_POSTO,
>>         p2min,
>>         d_D,d_TABLE,d_TYPE,
>>         &(in1[n*d_K*d_FSM.I()]),&(in2[n*d_K*d_D]),
>>         &(out[n*d_K*multiple])//,
>>         //d_alpha,d_beta
>>         );
>>     }
>>   }
>>
>>   for (unsigned int i = 0; i < input_items.size()/2; i++) {
>>     consume(2*i,d_FSM.I() * noutput_items / multiple );
>>     consume(2*i+1,d_D * noutput_items / multiple );
>>   }
>>
>>   return noutput_items;
>> }
>>
>>
>> /* -*- c++ -*- */
>> /*
>>  * Copyright 2004 Free Software Foundation, Inc.
>>  *
>>  * This file is part of GNU Radio
>>  *
>>  * GNU Radio is free software; you can redistribute it and/or modify
>>  * it under the terms of the GNU General Public License as published by
>>  * the Free Software Foundation; either version 2, or (at your option)
>>  * any later version.
>>  *
>>  * GNU Radio is distributed in the hope that it will be useful,
>>  * but WITHOUT ANY WARRANTY; without even the implied warranty of
>>  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
>>  * GNU General Public License for more details.
>>  *
>>  * You should have received a copy of the GNU General Public License
>>  * along with GNU Radio; see the file COPYING.  If not, write to
>>  * the Free Software Foundation, Inc., 51 Franklin Street,
>>  * Boston, MA 02110-1301, USA.
>>  */
>>
>> #ifdef HAVE_CONFIG_H
>> #include "config.h"
>> #endif
>>
>> #include <trellis_siso_f.h>
>> #include <gr_io_signature.h>
>> #include <stdexcept>
>> #include <assert.h>
>> #include <iostream>
>>
>> static const float INF = 1.0e9;
>>
>> trellis_siso_f_sptr
>> trellis_make_siso_f (
>>     const fsm &FSM,
>>     int K,
>>     int S0,
>>     int SK,
>>     bool POSTI,
>>     bool POSTO,
>>     trellis_siso_type_t SISO_TYPE)
>> {
>>   return trellis_siso_f_sptr (new trellis_siso_f
>> (FSM,K,S0,SK,POSTI,POSTO,SISO_TYPE));
>> }
>>
>> trellis_siso_f::trellis_siso_f (
>>     const fsm &FSM,
>>     int K,
>>     int S0,
>>     int SK,
>>     bool POSTI,
>>     bool POSTO,
>>     trellis_siso_type_t SISO_TYPE)
>>   : gr_block ("siso_f",
>>                           gr_make_io_signature (1, -1, sizeof (float)),
>>                           gr_make_io_signature (1, -1, sizeof (float))),
>>   d_FSM (FSM),
>>   d_K (K),
>>   d_S0 (S0),
>>   d_SK (SK),
>>   d_POSTI (POSTI),
>>   d_POSTO (POSTO),
>>   d_SISO_TYPE (SISO_TYPE)//,
>>   //d_alpha(FSM.S()*(K+1)),
>>   //d_beta(FSM.S()*(K+1))
>> {
>>     int multiple;
>>     if (d_POSTI && d_POSTO)
>>         multiple = d_FSM.I()+d_FSM.O();
>>     else if(d_POSTI)
>>         multiple = d_FSM.I();
>>     else if(d_POSTO)
>>         multiple = d_FSM.O();
>>     else
>>         throw std::runtime_error ("Not both POSTI and POSTO can be
>> false.");
>>     //printf("constructor: Multiple = %d\n",multiple);
>>     set_output_multiple (d_K*multiple);
>>     //what is the meaning of relative rate for a block with 2 inputs?
>>     //set_relative_rate ( multiple / ((double) d_FSM.I()) );
>>     // it turns out that the above gives problems in the scheduler, so
>>     // let's try (assumption O>I)
>>     //set_relative_rate ( multiple / ((double) d_FSM.O()) );
>>     // I am tempted to automate like this
>>     if(d_FSM.I() <= d_FSM.O())
>>       set_relative_rate ( multiple / ((double) d_FSM.O()) );
>>     else
>>       set_relative_rate ( multiple / ((double) d_FSM.I()) );
>> }
>>
>>
>> void
>> trellis_siso_f::forecast (int noutput_items, gr_vector_int
>> &ninput_items_required)
>> {
>>   int multiple;
>>   if (d_POSTI && d_POSTO)
>>       multiple = d_FSM.I()+d_FSM.O();
>>   else if(d_POSTI)
>>       multiple = d_FSM.I();
>>   else if(d_POSTO)
>>       multiple = d_FSM.O();
>>   else
>>       throw std::runtime_error ("Not both POSTI and POSTO can be
>> false.");
>>   //printf("forecast: Multiple = %d\n",multiple);
>>   assert (noutput_items % (d_K*multiple) == 0);
>>   int input_required1 =  d_FSM.I() * (noutput_items/multiple) ;
>>   int input_required2 =  d_FSM.O() * (noutput_items/multiple) ;
>>   //printf("forecast: Output requirements:  %d\n",noutput_items);
>>   //printf("forecast: Input requirements:  %d
>> %d\n",input_required1,input_required2);
>>   unsigned ninputs = ninput_items_required.size();
>>   assert(ninputs % 2 == 0);
>>   for (unsigned int i = 0; i < ninputs/2; i++) {
>>     ninput_items_required[2*i] = input_required1;
>>     ninput_items_required[2*i+1] = input_required2;
>>   }
>> }
>>
>> inline float min(float a, float b)
>> {
>>   return a <= b ? a : b;
>> }
>>
>> inline float min_star(float a, float b)
>> {
>>   return (a <= b ? a : b)-log(1+exp(a <= b ? a-b : b-a));
>> }
>>
>> void siso_algorithm(int I, int S, int O,
>>              const std::vector<int> &NS,
>>              const std::vector<int> &OS,
>>              const std::vector< std::vector<int> > &PS,
>>              const std::vector< std::vector<int> > &PI,
>>              int K,
>>              int S0,int SK,
>>              bool POSTI, bool POSTO,
>>              float (*p2mymin)(float,float),
>>              const float *priori, const float *prioro, float *post//,
>>              //std::vector<float> &alpha,
>>              //std::vector<float> &beta
>>              )
>> {
>>   float norm,mm,minm;
>>   std::vector<float> alpha(S*(K+1));
>>   std::vector<float> beta(S*(K+1));
>>
>>
>>   if(S0<0) { // initial state not specified
>>       for(int i=0;i<S;i++) alpha[0*S+i]=0;
>>   }
>>   else {
>>       for(int i=0;i<S;i++) alpha[0*S+i]=INF;
>>       alpha[0*S+S0]=0.0;
>>   }
>>
>>   for(int k=0;k<K;k++) { // forward recursion
>>       norm=INF;
>>       for(int j=0;j<S;j++) {
>>           minm=INF;
>>           for(int i=0;i<PS[j].size();i++) {
>>               //int i0 = j*I+i;
>>
>> 
mm=alpha[k*S+PS[j][i]]+priori[k*I+PI[j][i]]+prioro[k*O+OS[PS[j][i]*I+PI[j][i]]];
>>
>>               minm=(*p2mymin)(minm,mm);
>>           }
>>           alpha[(k+1)*S+j]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int j=0;j<S;j++)
>>           alpha[(k+1)*S+j]-=norm; // normalize total metrics so they
>> do not explode
>>   }
>>
>>   if(SK<0) { // final state not specified
>>       for(int i=0;i<S;i++) beta[K*S+i]=0;
>>   }
>>   else {
>>       for(int i=0;i<S;i++) beta[K*S+i]=INF;
>>       beta[K*S+SK]=0.0;
>>   }
>>
>>   for(int k=K-1;k>=0;k--) { // backward recursion
>>       norm=INF;
>>       for(int j=0;j<S;j++) {
>>           minm=INF;
>>           for(int i=0;i<I;i++) {
>>               int i0 = j*I+i;
>>               mm=beta[(k+1)*S+NS[i0]]+priori[k*I+i]+prioro[k*O+OS[i0]];
>>               minm=(*p2mymin)(minm,mm);
>>           }
>>           beta[k*S+j]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int j=0;j<S;j++)
>>           beta[k*S+j]-=norm; // normalize total metrics so they do not
>> explode
>>   }
>>
>>
>> if (POSTI && POSTO)
>> {
>>   for(int k=0;k<K;k++) { // input combining
>>       norm=INF;
>>       for(int i=0;i<I;i++) {
>>           minm=INF;
>>           for(int j=0;j<S;j++) {
>>
>> mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
>>               minm=(*p2mymin)(minm,mm);
>>           }
>>           post[k*(I+O)+i]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int i=0;i<I;i++)
>>           post[k*(I+O)+i]-=norm; // normalize metrics
>>   }
>>
>>
>>   for(int k=0;k<K;k++) { // output combining
>>       norm=INF;
>>       for(int n=0;n<O;n++) {
>>           minm=INF;
>>           for(int j=0;j<S;j++) {
>>               for(int i=0;i<I;i++) {
>>                   mm= (n==OS[j*I+i] ?
>> alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
>>                   minm=(*p2mymin)(minm,mm);
>>               }
>>           }
>>           post[k*(I+O)+I+n]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int n=0;n<O;n++)
>>           post[k*(I+O)+I+n]-=norm; // normalize metrics
>>   }
>> }
>> else if(POSTI)
>> {
>>   for(int k=0;k<K;k++) { // input combining
>>       norm=INF;
>>       for(int i=0;i<I;i++) {
>>           minm=INF;
>>           for(int j=0;j<S;j++) {
>>
>> mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
>>               minm=(*p2mymin)(minm,mm);
>>           }
>>           post[k*I+i]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int i=0;i<I;i++)
>>           post[k*I+i]-=norm; // normalize metrics
>>   }
>> }
>> else if(POSTO)
>> {
>>   for(int k=0;k<K;k++) { // output combining
>>       norm=INF;
>>       for(int n=0;n<O;n++) {
>>           minm=INF;
>>           for(int j=0;j<S;j++) {
>>               for(int i=0;i<I;i++) {
>>                   mm= (n==OS[j*I+i] ?
>> alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
>>                   minm=(*p2mymin)(minm,mm);
>>               }
>>           }
>>           post[k*O+n]=minm;
>>           if(minm<norm) norm=minm;
>>       }
>>       for(int n=0;n<O;n++)
>>           post[k*O+n]-=norm; // normalize metrics
>>   }
>> }
>> else
>>     throw std::runtime_error ("Not both POSTI and POSTO can be false.");
>>
>> }
>>
>>
>>
>>
>>
>>
>> int
>> trellis_siso_f::general_work (int noutput_items,
>>                         gr_vector_int &ninput_items,
>>                         gr_vector_const_void_star &input_items,
>>                         gr_vector_void_star &output_items)
>> {
>>   assert (input_items.size() == 2*output_items.size());
>>   int nstreams = output_items.size();
>>   //printf("general_work:Streams:  %d\n",nstreams);
>>   int multiple;
>>   if (d_POSTI && d_POSTO)
>>       multiple = d_FSM.I()+d_FSM.O();
>>   else if(d_POSTI)
>>       multiple = d_FSM.I();
>>   else if(d_POSTO)
>>       multiple = d_FSM.O();
>>   else
>>       throw std::runtime_error ("Not both POSTI and POSTO can be
>> false.");
>>
>>   assert (noutput_items % (d_K*multiple) == 0);
>>   int nblocks = noutput_items / (d_K*multiple);
>>   //printf("general_work:Blocks:  %d\n",nblocks);
>>   //for(int i=0;i<ninput_items.size();i++)
>>       //printf("general_work:Input items available:
>> %d\n",ninput_items[i]);
>>
>>   float (*p2min)(float, float) = NULL;
>>   if(d_SISO_TYPE == TRELLIS_MIN_SUM)
>>     p2min = &min;
>>   else if(d_SISO_TYPE == TRELLIS_SUM_PRODUCT)
>>     p2min = &min_star;
>>
>>
>>   for (int m=0;m<nstreams;m++) {
>>     const float *in1 = (const float *) input_items[2*m];
>>     const float *in2 = (const float *) input_items[2*m+1];
>>     float *out = (float *) output_items[m];
>>     for (int n=0;n<nblocks;n++) {
>>       siso_algorithm(d_FSM.I(),d_FSM.S(),d_FSM.O(),
>>         d_FSM.NS(),d_FSM.OS(),d_FSM.PS(),d_FSM.PI(),
>>         d_K,d_S0,d_SK,
>>         d_POSTI,d_POSTO,
>>         p2min,
>>         &(in1[n*d_K*d_FSM.I()]),&(in2[n*d_K*d_FSM.O()]),
>>         &(out[n*d_K*multiple])//,
>>         //d_alpha,d_beta
>>         );
>>     }
>>   }
>>
>>   for (unsigned int i = 0; i < input_items.size()/2; i++) {
>>     consume(2*i,d_FSM.I() * noutput_items / multiple );
>>     consume(2*i+1,d_FSM.O() * noutput_items / multiple );
>>   }
>>
>>   return noutput_items;
>> }
>>
>>
>> #!/usr/bin/env python
>>
>> from gnuradio import gr
>> from gnuradio import audio
>> from gnuradio import trellis
>> from gnuradio import eng_notation
>> import math
>> import sys
>> import random
>> import fsm_utils
>>
>> def run_test (f,Kb,bitspersymbol,K,dimensionality,constellation,N0,seed):
>>     fg = gr.flow_graph ()
>>
>>
>>     # TX
>>     #packet = [0]*Kb
>>     #for i in range(Kb-1*16): # last 16 bits = 0 to drive the final
>> state to 0
>>         #packet[i] = random.randint(0, 1) # random 0s and 1s
>>     #src = gr.vector_source_s(packet,False)
>>     src = gr.lfsr_32k_source_s()
>>     src_head = gr.head (gr.sizeof_short,Kb/16) # packet size in shorts
>>     #b2s = gr.unpacked_to_packed_ss(1,gr.GR_MSB_FIRST) # pack bits in
>> shorts
>>     s2fsmi = gr.packed_to_unpacked_ss(bitspersymbol,gr.GR_MSB_FIRST) #
>> unpack shorts to symbols compatible with the FSM input cardinality
>>     enc = trellis.encoder_ss(f,0) # initial state = 0
>>     mod = gr.chunks_to_symbols_sf(constellation,dimensionality)
>>
>>     # CHANNEL
>>     add = gr.add_ff()
>>     noise = gr.noise_source_f(gr.GR_GAUSSIAN,math.sqrt(N0/2),seed)
>>
>>     # RX
>>     metrics =
>> 
trellis.metrics_f(f.O(),dimensionality,constellation,trellis.TRELLIS_EUCLIDEAN)
>> # data preprocessing to generate metrics for Viterbi
>>     va = trellis.viterbi_s(f,K,0,-1) # Put -1 if the Initial/Final
>> states are not set.
>>     fsmi2s = gr.unpacked_to_packed_ss(bitspersymbol,gr.GR_MSB_FIRST) #
>> pack FSM input symbols to shorts
>>     #s2b = gr.packed_to_unpacked_ss(1,gr.GR_MSB_FIRST) # unpack shorts
>> to bits
>>     #dst = gr.vector_sink_s();
>>     dst = gr.check_lfsr_32k_s()
>>
>>
>>     fg.connect (src,src_head,s2fsmi,enc,mod)
>>     #fg.connect (src,b2s,s2fsmi,enc,mod)
>>     fg.connect (mod,(add,0))
>>     fg.connect (noise,(add,1))
>>     fg.connect (add,metrics)
>>     fg.connect (metrics,va,fsmi2s,dst)
>>     #fg.connect (metrics,va,fsmi2s,s2b,dst)
>>
>>
>>     fg.run()
>>
>>     # A bit of cheating: run the program once and print the
>>     # final encoder state..
>>     # Then put it as the last argument in the viterbi block
>>     #print "final state = " , enc.ST()
>>
>>     ntotal = dst.ntotal ()
>>     nright = dst.nright ()
>>     runlength = dst.runlength ()
>>     #ntotal = len(packet)
>>     #if len(dst.data()) != ntotal:
>>         #print "Error: not enough data\n"
>>     #nright = 0;
>>     #for i in range(ntotal):
>>         #if packet[i]==dst.data()[i]:
>>             #nright=nright+1
>>         #else:
>>             #print "Error in ", i
>>     return (ntotal,ntotal-nright)
>>
>>
>>
>>
>> def main(args):
>>     nargs = len (args)
>>     if nargs == 3:
>>         fname=args[0]
>>         esn0_db=float(args[1]) # Es/No in dB
>>         rep=int(args[2]) # number of times the experiment is run to
>> collect enough errors
>>     else:
>>         sys.stderr.write ('usage: test_tcm.py fsm_fname Es/No_db
>> repetitions\n')
>>         sys.exit (1)
>>
>>     # system parameters
>>     f=trellis.fsm(fname) # get the FSM specification from a file
>>     Kb=1024*16  # packet size in bits (make it multiple of 16 so it
>> can be packed in a short)
>>     bitspersymbol = int(round(math.log(f.I())/math.log(2))) # bits per
>> FSM input symbol
>>     K=Kb/bitspersymbol # packet size in trellis steps
>>     modulation = fsm_utils.pam2 # see fsm_utlis.py for available
>> predefined modulations
>>     dimensionality = modulation[0]
>>     constellation = modulation[1]
>>     if len(constellation)/dimensionality != f.O():
>>         sys.stderr.write ('Incompatible FSM output cardinality and
>> modulation size.\n')
>>         sys.exit (1)
>>     # calculate average symbol energy
>>     Es = 0
>>     for i in range(len(constellation)):
>>         Es = Es + constellation[i]**2
>>     Es = Es / (len(constellation)/dimensionality)
>>     N0=Es/pow(10.0,esn0_db/10.0); # calculate noise variance
>>
>>     tot_s=0 # total number of transmitted shorts
>>     terr_s=0 # total number of shorts in error
>>     terr_p=0 # total number of packets in error
>>     for i in range(rep):
>>
>> 
(s,e)=run_test(f,Kb,bitspersymbol,K,dimensionality,constellation,N0,-long(666+i))
>> # run experiment with different seed to get different noise realizations
>>         tot_s=tot_s+s
>>         terr_s=terr_s+e
>>         terr_p=terr_p+(terr_s!=0)
>>         if ((i+1)%100==0) : # display progress
>>             print i+1,terr_p, '%.2e' %
>> ((1.0*terr_p)/(i+1)),tot_s,terr_s, '%.2e' % ((1.0*terr_s)/tot_s)
>>     # estimate of the (short or bit) error rate
>>     print rep,terr_p, '%.2e' % ((1.0*terr_p)/(i+1)),tot_s,terr_s,
>> '%.2e' % ((1.0*terr_s)/tot_s)
>>
>>
>>
>> if __name__ == '__main__':
>>     main (sys.argv[1:])
>>
>>
>> 2 2 2
>>
>> 0 0
>> 0 1
>>
>> 0 1
>> 0 1
>>
>>
>> useless irregular FSM for testing. state 0 has 3 incoming edges and state
>> 1 has 1 incoming edge.
>>
>>
>> /* -*- c++ -*- */
>> /*
>>  * Copyright 2002 Free Software Foundation, Inc.
>>  *
>>  * This file is part of GNU Radio
>>  *
>>  * GNU Radio is free software; you can redistribute it and/or modify
>>  * it under the terms of the GNU General Public License as published by
>>  * the Free Software Foundation; either version 2, or (at your option)
>>  * any later version.
>>  *
>>  * GNU Radio is distributed in the hope that it will be useful,
>>  * but WITHOUT ANY WARRANTY; without even the implied warranty of
>>  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
>>  * GNU General Public License for more details.
>>  *
>>  * You should have received a copy of the GNU General Public License
>>  * along with GNU Radio; see the file COPYING.  If not, write to
>>  * the Free Software Foundation, Inc., 51 Franklin Street,
>>  * Boston, MA 02110-1301, USA.
>>  */
>>
>> #ifndef INCLUDED_TRELLIS_FSM_H
>> #define INCLUDED_TRELLIS_FSM_H
>>
>> #include <vector>
>>
>> /*!
>>  * \brief  FSM class
>>  */
>> class fsm {
>> private:
>>   int d_I;
>>   int d_S;
>>   int d_O;
>>   std::vector<int> d_NS;
>>   std::vector<int> d_OS;
>>   std::vector< std::vector<int> > d_PS;
>>   std::vector< std::vector<int> > d_PI;
>>   std::vector<int> d_TMi;
>>   std::vector<int> d_TMl;
>>   void generate_PS_PI ();
>>   void generate_TM ();
>>   bool find_es(int es);
>> public:
>>   fsm();
>>   fsm(const fsm &FSM);
>>   fsm(int I, int S, int O, const std::vector<int> &NS, const
>> std::vector<int> &OS);
>>   fsm(const char *name);
>>   fsm(int k, int n, const std::vector<int> &G);
>>   fsm(int mod_size, int ch_length);
>>   int I () const { return d_I; }
>>   int S () const { return d_S; }
>>   int O () const { return d_O; }
>>   const std::vector<int> & NS () const { return d_NS; }
>>   const std::vector<int> & OS () const { return d_OS; }
>>   const std::vector< std::vector<int> > & PS () const { return d_PS; }
>>   const std::vector< std::vector<int> > & PI () const { return d_PI; }
>>   const std::vector<int> & TMi () const { return d_TMi; }
>>   const std::vector<int> & TMl () const { return d_TMl; }
>> };
>>
>> #endif
>>
>>
>> /* -*- c++ -*- */
>> /*
>>  * Copyright 2002 Free Software Foundation, Inc.
>>  *
>>  * This file is part of GNU Radio
>>  *
>>  * GNU Radio is free software; you can redistribute it and/or modify
>>  * it under the terms of the GNU General Public License as published by
>>  * the Free Software Foundation; either version 2, or (at your option)
>>  * any later version.
>>  *
>>  * GNU Radio is distributed in the hope that it will be useful,
>>  * but WITHOUT ANY WARRANTY; without even the implied warranty of
>>  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
>>  * GNU General Public License for more details.
>>  *
>>  * You should have received a copy of the GNU General Public License
>>  * along with GNU Radio; see the file COPYING.  If not, write to
>>  * the Free Software Foundation, Inc., 51 Franklin Street,
>>  * Boston, MA 02110-1301, USA.
>>  */
>>
>> #include <cstdio>
>> #include <stdexcept>
>> #include <cmath>
>> #include "base.h"
>> #include "fsm.h"
>>
>>
>> fsm::fsm()
>> {
>>   d_I=0;
>>   d_S=0;
>>   d_O=0;
>>   d_NS.resize(0);
>>   d_OS.resize(0);
>>   d_PS.resize(0);
>>   d_PI.resize(0);
>>   d_TMi.resize(0);
>>   d_TMl.resize(0);
>> }
>>
>> fsm::fsm(const fsm &FSM)
>> {
>>   d_I=FSM.I();
>>   d_S=FSM.S();
>>   d_O=FSM.O();
>>   d_NS=FSM.NS();
>>   d_OS=FSM.OS();
>>   d_PS=FSM.PS(); // is this going to make a deep copy?
>>   d_PI=FSM.PI();
>>   d_TMi=FSM.TMi();
>>   d_TMl=FSM.TMl();
>> }
>>
>> fsm::fsm(int I, int S, int O, const std::vector<int> &NS, const
>> std::vector<int> &OS)
>> {
>>   d_I=I;
>>   d_S=S;
>>   d_O=O;
>>   d_NS=NS;
>>   d_OS=OS;
>>
>>   generate_PS_PI();
>>   generate_TM();
>> }
>>
>> //######################################################################
>> //# Read an FSM specification from a file.
>> //# Format (hopefully will become more flexible in the future...):
>> //# I S O (in the first line)
>> //# blank line
>> //# Next state matrix (S lines, each with I integers separated by spaces)
>> //# blank line
>> //# output symbol matrix (S lines, each with I integers separated by
>> spaces)
>> //# optional comments
>> //######################################################################
>> fsm::fsm(const char *name)
>> {
>>   FILE *fsmfile;
>>
>>   if((fsmfile=fopen(name,"r"))==NULL)
>>     throw std::runtime_error ("fsm::fsm(const char *name): file open
>> error\n");
>>     //printf("file open error in fsm()\n");
>>
>>   fscanf(fsmfile,"%d %d %d\n",&d_I,&d_S,&d_O);
>>   d_NS.resize(d_I*d_S);
>>   d_OS.resize(d_I*d_S);
>>
>>   for(int i=0;i<d_S;i++) {
>>     for(int j=0;j<d_I;j++) fscanf(fsmfile,"%d",&(d_NS[i*d_I+j]));
>>   }
>>   for(int i=0;i<d_S;i++) {
>>     for(int j=0;j<d_I;j++) fscanf(fsmfile,"%d",&(d_OS[i*d_I+j]));
>>   }
>>
>>   generate_PS_PI();
>>   generate_TM();
>> }
>>
>>
>>
>>
>> //######################################################################
>> //# Automatically generate the FSM from the generator matrix
>> //# of a (n,k) binary convolutional code
>> //######################################################################
>> fsm::fsm(int k, int n, const std::vector<int> &G)
>> {
>>
>>   // calculate maximum memory requirements for each input stream
>>   std::vector<int> max_mem_x(k,-1);
>>   int max_mem = -1;
>>   for(int i=0;i<k;i++) {
>>     for(int j=0;j<n;j++) {
>>       int mem = -1;
>>       if(G[i*n+j]!=0)
>>         mem=(int)(log(G[i*n+j])/log(2.0));
>>       if(mem>max_mem_x[i])
>>         max_mem_x[i]=mem;
>>       if(mem>max_mem)
>>         max_mem=mem;
>>     }
>>   }
>>
>> //printf("max_mem_x\n");
>> //for(int j=0;j<max_mem_x.size();j++) printf("%d ",max_mem_x[j]);
>> printf("\n");
>>
>>   // calculate total memory requirements to set S
>>   int sum_max_mem = 0;
>>   for(int i=0;i<k;i++)
>>     sum_max_mem += max_mem_x[i];
>>
>> //printf("sum_max_mem = %d\n",sum_max_mem);
>>
>>   d_I=1<<k;
>>   d_S=1<<sum_max_mem;
>>   d_O=1<<n;
>>
>>   // binary representation of the G matrix
>>   std::vector<std::vector<int> > Gb(k*n);
>>   for(int j=0;j<k*n;j++) {
>>     Gb[j].resize(max_mem+1);
>>     dec2base(G[j],2,Gb[j]);
>> //printf("Gb\n");
>> //for(int m=0;m<Gb[j].size();m++) printf("%d ",Gb[j][m]); printf("\n");
>>   }
>>
>>   // alphabet size of each shift register
>>   std::vector<int> bases_x(k);
>>   for(int j=0;j<k ;j++)
>>     bases_x[j] = 1 << max_mem_x[j];
>> //printf("bases_x\n");
>> //for(int j=0;j<max_mem_x.size();j++) printf("%d ",max_mem_x[j]);
>> printf("\n");
>>
>>   d_NS.resize(d_I*d_S);
>>   d_OS.resize(d_I*d_S);
>>
>>   std::vector<int> sx(k);
>>   std::vector<int> nsx(k);
>>   std::vector<int> tx(k);
>>   std::vector<std::vector<int> > tb(k);
>>   for(int j=0;j<k;j++)
>>     tb[j].resize(max_mem+1);
>>   std::vector<int> inb(k);
>>   std::vector<int> outb(n);
>>
>>
>>   for(int s=0;s<d_S;s++) {
>>     dec2bases(s,bases_x,sx); // split s into k values, each
>> representing on of the k shift registers
>> //printf("state = %d \nstates = ",s);
>> //for(int j=0;j<sx.size();j++) printf("%d ",sx[j]); printf("\n");
>>     for(int i=0;i<d_I;i++) {
>>       dec2base(i,2,inb); // input in binary
>> //printf("input = %d \ninputs = ",i);
>> //for(int j=0;j<inb.size();j++) printf("%d ",inb[j]); printf("\n");
>>
>>       // evaluate next state
>>       for(int j=0;j<k;j++)
>>         nsx[j] = (inb[j]*bases_x[j]+sx[j])/2; // next state (for each
>> shift register) MSB first
>>       d_NS[s*d_I+i]=bases2dec(nsx,bases_x); // collect all values into
>> the new state
>>
>>       // evaluate transitions
>>       for(int j=0;j<k;j++)
>>         tx[j] = inb[j]*bases_x[j]+sx[j]; // transition (for each shift
>> register)MSB first
>>       for(int j=0;j<k;j++) {
>>         dec2base(tx[j],2,tb[j]); // transition in binary
>> //printf("transition = %d \ntransitions = ",tx[j]);
>> //for(int m=0;m<tb[j].size();m++) printf("%d ",tb[j][m]); printf("\n");
>>       }
>>
>>       // evaluate outputs
>>       for(int nn=0;nn<n;nn++) {
>>         outb[nn] = 0;
>>         for(int j=0;j<k;j++) {
>>           for(int m=0;m<max_mem+1;m++)
>>             outb[nn] = (outb[nn] + Gb[j*n+nn][m]*tb[j][m]) % 2; //
>> careful: polynomial 1+D ir represented as 110, not as 011
>> //printf("output %d equals %d\n",nn,outb[nn]);
>>         }
>>       }
>>       d_OS[s*d_I+i] = base2dec(outb,2);
>>     }
>>   }
>>
>>   generate_PS_PI();
>>   generate_TM();
>> }
>>
>>
>>
>>
>> //######################################################################
>> //# Automatically generate an FSM specification describing the
>> //# ISI for a channel
>> //# of length ch_length and a modulation of size mod_size
>> //######################################################################
>> fsm::fsm(int mod_size, int ch_length)
>> {
>>   d_I=mod_size;
>>   d_S=(int) (pow(1.0*d_I,1.0*ch_length-1)+0.5);
>>   d_O=d_S*d_I;
>>
>>   d_NS.resize(d_I*d_S);
>>   d_OS.resize(d_I*d_S);
>>
>>   for(int s=0;s<d_S;s++) {
>>     for(int i=0;i<d_I;i++) {
>>       int t=i*d_S+s;
>>       d_NS[s*d_I+i] = t/d_I;
>>       d_OS[s*d_I+i] = t;
>>     }
>>   }
>>
>>   generate_PS_PI();
>>   generate_TM();
>> }
>>
>>
>> //######################################################################
>> //# generate the PS and PI tables for later use
>> //######################################################################
>> void fsm::generate_PS_PI()
>> {
>>   d_PS.resize(d_S);
>>   d_PI.resize(d_S);
>>
>>   for(int i=0;i<d_S;i++) {
>>     d_PS[i].resize(d_I*d_S); // max possible size
>>     d_PI[i].resize(d_I*d_S);
>>     int j=0;
>>     for(int ii=0;ii<d_S;ii++) for(int jj=0;jj<d_I;jj++) {
>>       if(d_NS[ii*d_I+jj]!=i) continue;
>>       d_PS[i][j]=ii;
>>       d_PI[i][j]=jj;
>>       j++;
>>     }
>>     d_PS[i].resize(j);
>>     d_PI[i].resize(j);
>>   }
>> }
>>
>>
>> //######################################################################
>> //# generate the termination matrices TMl and TMi for later use
>> //######################################################################
>> void fsm::generate_TM()
>> {
>>   d_TMi.resize(d_S*d_S);
>>   d_TMl.resize(d_S*d_S);
>>
>>   for(int i=0;i<d_S*d_S;i++) {
>>     d_TMi[i] = -1; // no meaning
>>     d_TMl[i] = d_S; //infinity: you need at most S-1 steps
>>     if (i/d_S == i%d_S)
>>       d_TMl[i] = 0;
>>   }
>>
>>   for(int s=0;s<d_S;s++) {
>>     bool done = false;
>>     int attempts = 0;
>>     while (done == false && attempts < d_S-1) {
>>       done = find_es(s);
>>       attempts ++;
>>     }
>>     if (done == false)
>>       //throw std::runtime_error ("fsm::generate_TM(): FSM appears to
>> be disconnected\n");
>>       printf("fsm::generate_TM(): FSM appears to be disconnected\n");
>>   }
>> }
>>
>>
>> // find a path from any state to the ending state "es"
>> bool fsm::find_es(int es)
>> {
>>   bool done = true;
>>   for(int s=0;s<d_S;s++) {
>>     if(d_TMl[s*d_S+es] < d_S)
>>       continue;
>>     int minl=d_S;
>>     int mini=-1;
>>     for(int i=0;i<d_I;i++) {
>>       if( 1 + d_TMl[d_NS[s*d_I+i]*d_S+es] < minl) {
>>         minl = 1 + d_TMl[d_NS[s*d_I+i]*d_S+es];
>>         mini = i;
>>       }
>>     }
>>     if (mini != -1) {
>>       d_TMl[s*d_S+es]=minl;
>>       d_TMi[s*d_S+es]=mini;
>>     }
>>     else
>>       done = false;
>>   }
>>   return done;
>> }
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> /* -*- c++ -*- */
>> /*
>>  * Copyright 2002 Free Software Foundation, Inc.
>>  *
>>  * This file is part of GNU Radio
>>  *
>>  * GNU Radio is free software; you can redistribute it and/or modify
>>  * it under the terms of the GNU General Public License as published by
>>  * the Free Software Foundation; either version 2, or (at your option)
>>  * any later version.
>>  *
>>  * GNU Radio is distributed in the hope that it will be useful,
>>  * but WITHOUT ANY WARRANTY; without even the implied warranty of
>>  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
>>  * GNU General Public License for more details.
>>  *
>>  * You should have received a copy of the GNU General Public License
>>  * along with GNU Radio; see the file COPYING.  If not, write to
>>  * the Free Software Foundation, Inc., 51 Franklin Street,
>>  * Boston, MA 02110-1301, USA.
>>  */
>>
>> class fsm {
>> private:
>>   int d_I;
>>   int d_S;
>>   int d_O;
>>   std::vector<int> d_NS;
>>   std::vector<int> d_OS;
>>   std::vector< std::vector<int> > d_PS;
>>   std::vector< std::vector<int> > d_PI;
>>   std::vector<int> d_TMi;
>>   std::vector<int> d_TMl;
>>   void generate_PS_PI ();
>>   void generate_TM ();
>> public:
>>   fsm();
>>   fsm(const fsm &FSM);
>>   fsm(int I, int S, int O, const std::vector<int> &NS, const
>> std::vector<int> &OS);
>>   fsm(const char *name);
>>   fsm(int k, int n, const std::vector<int> &G);
>>   fsm(int mod_size, int ch_length);
>>   int I () const { return d_I; }
>>   int S () const { return d_S; }
>>   int O () const { return d_O; }
>>   const std::vector<int> & NS () const { return d_NS; }
>>   const std::vector<int> & OS () const { return d_OS; }
>>   const std::vector< std::vector<int> > & PS () const { return d_PS; }
>>   const std::vector< std::vector<int> > & PI () const { return d_PI; }
>>   const std::vector<int> & TMi () const { return d_TMi; }
>>   const std::vector<int> & TMl () const { return d_TMl; }
>> };
>>
>>
>>
>>
>>
>
>

--
_______________________________________________________
Achilleas Anastasopoulos
Associate Professor
EECS Department               Voice : (734)615-4024
UNIVERSITY OF MICHIGAN        Fax   : (734)763-8041
Ann Arbor, MI 48109-2122      E-mail: address@hidden
URL: http://www-personal.engin.umich.edu/~anastas/
_______________________________________________________





reply via email to

[Prev in Thread] Current Thread [Next in Thread]