[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Discuss-gnuradio] Question on gr-trellis PS and PI matricies in FSm
|
From: |
Achilleas Anastasopoulos |
|
Subject: |
Re: [Discuss-gnuradio] Question on gr-trellis PS and PI matricies in FSm |
|
Date: |
Thu, 01 Feb 2007 14:50:59 -0500 |
|
User-agent: |
Mozilla Thunderbird 1.0.2 (Windows/20050317) |
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 <@address@hidden>
#include <gr_io_signature.h>
#include <assert.h>
#include <iostream>
static const float INF = 1.0e9;
@SPTR_NAME@
address@hidden@ (
const fsm &FSM,
int K,
int S0,
int SK)
{
return @SPTR_NAME@ (new @NAME@ (FSM,K,S0,SK));
}
@NAME@::@NAME@ (
const fsm &FSM,
int K,
int S0,
int SK)
: gr_block ("@BASE_NAME@",
gr_make_io_signature (1, -1, sizeof (float)),
gr_make_io_signature (1, -1, sizeof (@TYPE@))),
d_FSM (FSM),
d_K (K),
d_S0 (S0),
d_SK (SK)//,
//d_trace(FSM.S()*K)
{
set_relative_rate (1.0 / ((double) d_FSM.O()));
set_output_multiple (d_K);
}
void
@NAME@::forecast (int noutput_items, gr_vector_int &ninput_items_required)
{
assert (noutput_items % d_K == 0);
int input_required = d_FSM.O() * noutput_items ;
unsigned ninputs = ninput_items_required.size();
for (unsigned int i = 0; i < ninputs; i++) {
ninput_items_required[i] = input_required;
}
}
void viterbi_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,
const float *in, @TYPE@ *out)//,
//std::vector<int> &trace)
{
std::vector<int> trace(S*K);
std::vector<float> alpha(S*2);
int alphai;
float norm,mm,minm;
int minmi;
int st;
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;
}
alphai=0;
for(int k=0;k<K;k++) {
norm=INF;
for(int j=0;j<S;j++) { // for each next state do ACS
minm=INF;
minmi=0;
for(int i=0;i<PS[j].size();i++) {
//int i0 = j*I+i;
if((mm=alpha[alphai*S+PS[j][i]]+in[k*O+OS[PS[j][i]*I+PI[j][i]]])<minm)
minm=mm,minmi=i;
}
trace[k*S+j]=minmi;
alpha[((alphai+1)%2)*S+j]=minm;
if(minm<norm) norm=minm;
}
for(int j=0;j<S;j++)
alpha[((alphai+1)%2)*S+j]-=norm; // normalize total metrics so they
do not explode
alphai=(alphai+1)%2;
}
if(SK<0) { // final state not specified
minm=INF;
minmi=0;
for(int i=0;i<S;i++)
if((mm=alpha[alphai*S+i])<minm) minm=mm,minmi=i;
st=minmi;
}
else {
st=SK;
}
for(int k=K-1;k>=0;k--) { // traceback
int i0=trace[k*S+st];
out[k]= (@TYPE@) PI[st][i0];
st=PS[st][i0];
}
}
int
@NAME@::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() == output_items.size());
int nstreams = input_items.size();
assert (noutput_items % d_K == 0);
int nblocks = noutput_items / d_K;
for (int m=0;m<nstreams;m++) {
const float *in = (const float *) input_items[m];
@TYPE@ *out = (@TYPE@ *) output_items[m];
for (int n=0;n<nblocks;n++) {
viterbi_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,&(in[n*d_K*d_FSM.O()]),&(out[n*d_K]));//,d_trace);
}
}
consume_each (d_FSM.O() * noutput_items );
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_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; }
};
/* -*- 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 <@address@hidden>
#include <gr_io_signature.h>
#include <assert.h>
#include <iostream>
static const float INF = 1.0e9;
@SPTR_NAME@
address@hidden@ (
const fsm &FSM,
int K,
int S0,
int SK,
int D,
const std::vector<float> &TABLE,
trellis_metric_type_t TYPE)
{
return @SPTR_NAME@ (new @NAME@ (FSM,K,S0,SK,D,TABLE,TYPE));
}
@NAME@::@NAME@ (
const fsm &FSM,
int K,
int S0,
int SK,
int D,
const std::vector<float> &TABLE,
trellis_metric_type_t TYPE)
: gr_block ("@BASE_NAME@",
gr_make_io_signature (1, -1, sizeof (float)),
gr_make_io_signature (1, -1, sizeof (@TYPE@))),
d_FSM (FSM),
d_K (K),
d_S0 (S0),
d_SK (SK),
d_D (D),
d_TABLE (TABLE),
d_TYPE (TYPE)//,
//d_trace(FSM.S()*K)
{
set_relative_rate (1.0 / ((double) d_D));
set_output_multiple (d_K);
}
void
@NAME@::forecast (int noutput_items, gr_vector_int &ninput_items_required)
{
assert (noutput_items % d_K == 0);
int input_required = d_D * noutput_items ;
unsigned ninputs = ninput_items_required.size();
for (unsigned int i = 0; i < ninputs; i++) {
ninput_items_required[i] = input_required;
}
}
void viterbi_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,
int D,
const std::vector<float> &TABLE,
trellis_metric_type_t TYPE,
const float *in, @TYPE@ *out)//,
//std::vector<int> &trace)
{
std::vector<int> trace(S*K);
std::vector<float> alpha(S*2);
float *metric = new float[O];
int alphai;
float norm,mm,minm;
int minmi;
int st;
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;
}
alphai=0;
for(int k=0;k<K;k++) {
calc_metric(O, D, TABLE, &(in[k*D]), metric,TYPE); // calc metrics
norm=INF;
for(int j=0;j<S;j++) { // for each next state do ACS
minm=INF;
minmi=0;
for(int i=0;i<PS[j].size();i++) {
int i0 = j*I+i;
if((mm=alpha[alphai*S+PS[j][i]]+metric[OS[PS[j][i]*I+PI[j][i]]])<minm)
minm=mm,minmi=i;
}
trace[k*S+j]=minmi;
alpha[((alphai+1)%2)*S+j]=minm;
if(minm<norm) norm=minm;
}
for(int j=0;j<S;j++)
alpha[((alphai+1)%2)*S+j]-=norm; // normalize total metrics so they
do not explode
alphai=(alphai+1)%2;
}
if(SK<0) { // final state not specified
minm=INF;
minmi=0;
for(int i=0;i<S;i++)
if((mm=alpha[alphai*S+i])<minm) minm=mm,minmi=i;
st=minmi;
}
else {
st=SK;
}
for(int k=K-1;k>=0;k--) { // traceback
int i0=trace[k*S+st];
out[k]= (@TYPE@) PI[st][i0];
st=PS[st][i0];
}
delete [] metric;
}
int
@NAME@::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() == output_items.size());
int nstreams = input_items.size();
assert (noutput_items % d_K == 0);
int nblocks = noutput_items / d_K;
for (int m=0;m<nstreams;m++) {
const float *in = (const float *) input_items[m];
@TYPE@ *out = (@TYPE@ *) output_items[m];
for (int n=0;n<nblocks;n++) {
viterbi_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_D,d_TABLE,d_TYPE,&(in[n*d_K*d_D]),&(out[n*d_K]));//,d_trace);
}
}
consume_each (d_D * noutput_items );
return noutput_items;
}