[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Commit-gnuradio] r10514 - in gnuradio/trunk/gr-trellis/src: examples li
|
From: |
anastas |
|
Subject: |
[Commit-gnuradio] r10514 - in gnuradio/trunk/gr-trellis/src: examples lib |
|
Date: |
Wed, 25 Feb 2009 17:29:51 -0700 (MST) |
Author: anastas
Date: 2009-02-25 17:29:51 -0700 (Wed, 25 Feb 2009)
New Revision: 10514
Added:
gnuradio/trunk/gr-trellis/src/examples/test_cpm.py
Modified:
gnuradio/trunk/gr-trellis/src/examples/fsm_utils.py
gnuradio/trunk/gr-trellis/src/lib/fsm.cc
gnuradio/trunk/gr-trellis/src/lib/fsm.h
gnuradio/trunk/gr-trellis/src/lib/fsm.i
gnuradio/trunk/gr-trellis/src/lib/trellis_calc_metric.cc
Log:
Added support for Continuous Phase Modulation in gr-trellis + an example
Modified: gnuradio/trunk/gr-trellis/src/examples/fsm_utils.py
===================================================================
--- gnuradio/trunk/gr-trellis/src/examples/fsm_utils.py 2009-02-25 23:38:17 UTC
(rev 10513)
+++ gnuradio/trunk/gr-trellis/src/examples/fsm_utils.py 2009-02-26 00:29:51 UTC
(rev 10514)
@@ -25,6 +25,8 @@
import math
import sys
import operator
+import numpy
+import scipy.linalg
from gnuradio import trellis
@@ -58,34 +60,7 @@
return num
-######################################################################
-# Generate a new FSM representing the concatenation of two FSMs
-######################################################################
-def fsm_concatenate(f1,f2):
- if f1.O() > f2.I():
- print "Not compatible FSMs\n"
- I=f1.I()
- S=f1.S()*f2.S()
- O=f2.O()
- nsm=list([0]*I*S)
- osm=list([0]*I*S)
- for s1 in range(f1.S()):
- for s2 in range(f2.S()):
- for i in range(f1.I()):
- ns1=f1.NS()[s1*f1.I()+i]
- o1=f1.OS()[s1*f1.I()+i]
- ns2=f2.NS()[s2*f2.I()+o1]
- o2=f2.OS()[s2*f2.I()+o1]
- s=s1*f2.S()+s2
- ns=ns1*f2.S()+ns2
- nsm[s*I+i]=ns
- osm[s*I+i]=o2
-
-
- f=trellis.fsm(I,S,O,nsm,osm)
- return f
-
######################################################################
# Generate a new FSM representing n stages through the original FSM
######################################################################
@@ -143,10 +118,80 @@
return (1,lookup)
+
+
+
+
+######################################################################
+# Automatically generate the signals appropriate for CPM
+# decomposition.
+# This decomposition is based on the paper by B. Rimoldi
+# "A decomposition approach to CPM", IEEE Trans. Info Theory, March 1988
+# See also my own notes at http://www.eecs.umich.edu/~anastas/docs/cpm.pdf
+######################################################################
+def make_cpm_signals(K,P,M,L,q,frac):
+
+ Q=numpy.size(q)/L
+ h=(1.0*K)/P
+ f0=-h*(M-1)/2
+ dt=0.0; # maybe start at t=0.5
+ t=(dt+numpy.arange(0,Q))/Q
+ qq=numpy.zeros(Q)
+ for m in range(L):
+ qq=qq + q[m*Q:m*Q+Q]
+ w=math.pi*h*(M-1)*t-2*math.pi*h*(M-1)*qq+math.pi*h*(L-1)*(M-1)
+ X=(M**L)*P
+ PSI=numpy.empty((X,Q))
+ for x in range(X):
+ xv=dec2base(x/P,M,L)
+ xv=numpy.append(xv, x%P)
+ qq1=numpy.zeros(Q)
+ for m in range(L):
+ qq1=qq1+xv[m]*q[m*Q:m*Q+Q]
+ psi=2*math.pi*h*xv[-1]+4*math.pi*h*qq1+w
+ #print psi
+ PSI[x]=psi
+ PSI = numpy.transpose(PSI)
+ SS=numpy.exp(1j*PSI) # contains all signals as columns
+ #print SS
+
+ # Now we need to orthogonalize the signals
+ F = scipy.linalg.orth(SS) # find an orthonormal basis for SS
+ #print numpy.dot(numpy.transpose(F.conjugate()),F) # check for
orthonormality
+ S = numpy.dot(numpy.transpose(F.conjugate()),SS)
+ #print F
+ #print S
+ # We only want to keep those dimensions that contain most
+ # of the energy of the overall constellation (eg, frac=0.9 ==> 90%)
+ # evaluate mean energy in each dimension
+ E=numpy.sum(numpy.absolute(S)**2,axis=1)/Q
+ E=E/numpy.sum(E)
+ #print E
+ Es = -numpy.sort(-E)
+ Esi = numpy.argsort(-E)
+ #print Es
+ #print Esi
+ Ecum=numpy.cumsum(Es)
+ #print Ecum
+ v0=numpy.searchsorted(Ecum,frac)
+ N = v0+1
+ #print v0
+ #print Esi[0:v0+1]
+ Ff=numpy.transpose(numpy.transpose(F)[Esi[0:v0+1]])
+ #print Ff
+ Sf = S[Esi[0:v0+1]]
+ #print Sf
+
+ return (f0,SS,S,F,Sf,Ff,N)
+ #return f0
+
+
+
+
######################################################################
# A list of common modulations.
# Format: (dimensionality,constellation)
@@ -194,20 +239,23 @@
if __name__ == '__main__':
f1=trellis.fsm('fsm_files/awgn1o2_4.fsm')
#f2=trellis.fsm('fsm_files/awgn2o3_4.fsm')
- print f1.I(), f1.S(), f1.O()
- print f1.NS()
- print f1.OS()
+ #print f1.I(), f1.S(), f1.O()
+ #print f1.NS()
+ #print f1.OS()
#print f2.I(), f2.S(), f2.O()
#print f2.NS()
#print f2.OS()
##f1.write_trellis_svg('f1.svg',4)
#f2.write_trellis_svg('f2.svg',4)
#f=fsm_concatenate(f1,f2)
- f=fsm_radix(f1,2)
+ #f=fsm_radix(f1,2)
- print "----------\n"
- print f.I(), f.S(), f.O()
- print f.NS()
- print f.OS()
+ #print "----------\n"
+ #print f.I(), f.S(), f.O()
+ #print f.NS()
+ #print f.OS()
#f.write_trellis_svg('f.svg',4)
+ q=numpy.arange(0,8)/(2.0*8)
+ (f0,SS,S,F,Sf,Ff,N) = make_cpm_signals(1,2,2,1,q,0.99)
+
Added: gnuradio/trunk/gr-trellis/src/examples/test_cpm.py
===================================================================
--- gnuradio/trunk/gr-trellis/src/examples/test_cpm.py
(rev 0)
+++ gnuradio/trunk/gr-trellis/src/examples/test_cpm.py 2009-02-26 00:29:51 UTC
(rev 10514)
@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+##################################################
+# Gnuradio Python Flow Graph
+# Title: CPM test
+# Author: Achilleas Anastasopoulos
+# Description: gnuradio flow graph
+# Generated: Thu Feb 19 23:16:23 2009
+##################################################
+
+from gnuradio import gr
+from gnuradio import trellis
+from gnuradio.gr import firdes
+from grc_gnuradio import blks2 as grc_blks2
+import math
+import numpy
+import scipy.stats
+import fsm_utils
+from gnuradio import trellis
+
+def run_test(seed,blocksize):
+ tb = gr.top_block()
+
+ ##################################################
+ # Variables
+ ##################################################
+ M = 2
+ K = 1
+ P = 2
+ h = (1.0*K)/P
+ L = 3
+ Q = 4
+ frac = 0.99
+ f = trellis.fsm(P,M,L)
+
+ # CPFSK signals
+ #p = numpy.ones(Q)/(2.0)
+ #q = numpy.cumsum(p)/(1.0*Q)
+
+ # GMSK signals
+ BT=0.3;
+ tt=numpy.arange(0,L*Q)/(1.0*Q)-L/2.0;
+ #print tt
+
p=(0.5*scipy.stats.erfc(2*math.pi*BT*(tt-0.5)/math.sqrt(math.log(2.0))/math.sqrt(2.0))-0.5*scipy.stats.erfc(2*math.pi*BT*(tt+0.5)/math.sqrt(math.log(2.0))/math.sqrt(2.0)))/2.0;
+ p=p/sum(p)*Q/2.0;
+ #print p
+ q=numpy.cumsum(p)/Q;
+ q=q/q[-1]/2.0;
+ #print q
+
+ (f0T,SS,S,F,Sf,Ff,N) = fsm_utils.make_cpm_signals(K,P,M,L,q,frac)
+ #print N
+ #print Ff
+ Ffa = numpy.insert(Ff,Q,numpy.zeros(N),axis=0)
+ #print Ffa
+ MF = numpy.fliplr(numpy.transpose(Ffa))
+ #print MF
+ E = numpy.sum(numpy.abs(Sf)**2,axis=0)
+ Es = numpy.sum(E)/f.O()
+ #print Es
+
+ constellation = numpy.reshape(numpy.transpose(Sf),N*f.O())
+ #print Ff
+ #print Sf
+ #print constellation
+ #print numpy.max(numpy.abs(SS - numpy.dot(Ff , Sf)))
+
+ EsN0_db = 10.0
+ N0 = Es * 10.0**(-(1.0*EsN0_db)/10.0)
+ #N0 = 0.0
+ #print N0
+ head = 4
+ tail = 4
+ numpy.random.seed(seed*666)
+ data = numpy.random.randint(0, M, head+blocksize+tail+1)
+ #data = numpy.zeros(blocksize+1+head+tail,'int')
+ for i in range(head):
+ data[i]=0
+ for i in range(tail+1):
+ data[-i]=0
+
+
+
+ ##################################################
+ # Blocks
+ ##################################################
+ random_source_x_0 = gr.vector_source_b(data, False)
+ gr_chunks_to_symbols_xx_0 = gr.chunks_to_symbols_bf((-1, 1), 1)
+ gr_interp_fir_filter_xxx_0 = gr.interp_fir_filter_fff(Q, p)
+ gr_frequency_modulator_fc_0 =
gr.frequency_modulator_fc(2*math.pi*h*(1.0/Q))
+
+ gr_add_vxx_0 = gr.add_vcc(1)
+ gr_noise_source_x_0 = gr.noise_source_c(gr.GR_GAUSSIAN, (N0/2.0)**0.5,
-long(seed))
+
+ gr_multiply_vxx_0 = gr.multiply_vcc(1)
+ gr_sig_source_x_0 = gr.sig_source_c(Q, gr.GR_COS_WAVE, -f0T, 1, 0)
+ # only works for N=2, do it manually for N>2...
+ gr_fir_filter_xxx_0_0 = gr.fir_filter_ccc(Q, MF[0].conjugate())
+ gr_fir_filter_xxx_0_0_0 = gr.fir_filter_ccc(Q, MF[1].conjugate())
+ gr_streams_to_stream_0 = gr.streams_to_stream(gr.sizeof_gr_complex*1, N)
+ gr_skiphead_0 = gr.skiphead(gr.sizeof_gr_complex*1, N*(1+0))
+ viterbi = trellis.viterbi_combined_cb(f, head+blocksize+tail, 0, -1, N,
constellation, trellis.TRELLIS_EUCLIDEAN)
+
+ gr_vector_sink_x_0 = gr.vector_sink_b()
+
+ ##################################################
+ # Connections
+ ##################################################
+ tb.connect((random_source_x_0, 0), (gr_chunks_to_symbols_xx_0, 0))
+ tb.connect((gr_chunks_to_symbols_xx_0, 0), (gr_interp_fir_filter_xxx_0,
0))
+ tb.connect((gr_interp_fir_filter_xxx_0, 0),
(gr_frequency_modulator_fc_0, 0))
+ tb.connect((gr_frequency_modulator_fc_0, 0), (gr_add_vxx_0, 0))
+ tb.connect((gr_noise_source_x_0, 0), (gr_add_vxx_0, 1))
+ tb.connect((gr_add_vxx_0, 0), (gr_multiply_vxx_0, 0))
+ tb.connect((gr_sig_source_x_0, 0), (gr_multiply_vxx_0, 1))
+ tb.connect((gr_multiply_vxx_0, 0), (gr_fir_filter_xxx_0_0, 0))
+ tb.connect((gr_multiply_vxx_0, 0), (gr_fir_filter_xxx_0_0_0, 0))
+ tb.connect((gr_fir_filter_xxx_0_0, 0), (gr_streams_to_stream_0, 0))
+ tb.connect((gr_fir_filter_xxx_0_0_0, 0), (gr_streams_to_stream_0, 1))
+ tb.connect((gr_streams_to_stream_0, 0), (gr_skiphead_0, 0))
+ tb.connect((gr_skiphead_0, 0), (viterbi, 0))
+ tb.connect((viterbi, 0), (gr_vector_sink_x_0, 0))
+
+
+ tb.run()
+ dataest = gr_vector_sink_x_0.data()
+ #print data
+ #print numpy.array(dataest)
+ perr = 0
+ err = 0
+ for i in range(blocksize):
+ if data[head+i] != dataest[head+i]:
+ #print i
+ err += 1
+ if err != 0 :
+ perr = 1
+ return (err,perr)
+
+if __name__ == '__main__':
+ blocksize = 1000
+ ss=0
+ ee=0
+ for i in range(10000):
+ (s,e) = run_test(i,blocksize)
+ ss += s
+ ee += e
+ if (i+1) % 100 == 0:
+ print i+1,ss,ee,(1.0*ss)/(i+1)/(1.0*blocksize),(1.0*ee)/(i+1)
+ print i+1,ss,ee,(1.0*ss)/(i+1)/(1.0*blocksize),(1.0*ee)/(i+1)
Property changes on: gnuradio/trunk/gr-trellis/src/examples/test_cpm.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: gnuradio/trunk/gr-trellis/src/lib/fsm.cc
===================================================================
--- gnuradio/trunk/gr-trellis/src/lib/fsm.cc 2009-02-25 23:38:17 UTC (rev
10513)
+++ gnuradio/trunk/gr-trellis/src/lib/fsm.cc 2009-02-26 00:29:51 UTC (rev
10514)
@@ -170,7 +170,7 @@
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
+ dec2bases(s,bases_x,sx); // split s into k values, each representing one
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++) {
@@ -239,8 +239,55 @@
}
+
+
//######################################################################
//# Automatically generate an FSM specification describing the
+//# the trellis for a CPM with h=K/P (relatively prime),
+//# alphabet size M, and frequency pulse duration L symbols
+//#
+//# This FSM is based on the paper by B. Rimoldi
+//# "A decomposition approach to CPM", IEEE Trans. Info Theory, March 1988
+//# See also my own notes at http://www.eecs.umich.edu/~anastas/docs/cpm.pdf
+//######################################################################
+fsm::fsm(int P, int M, int L)
+{
+ d_I=M;
+ d_S=(int)(pow(1.0*M,1.0*L-1)+0.5)*P;
+ d_O=(int)(pow(1.0*M,1.0*L)+0.5)*P;
+
+ d_NS.resize(d_I*d_S);
+ d_OS.resize(d_I*d_S);
+ int nv;
+ for(int s=0;s<d_S;s++) {
+ for(int i=0;i<d_I;i++) {
+ int s1=s/P;
+ int v=s%P;
+ int ns1= (i*(int)(pow(1.0*M,1.0*(L-1))+0.5)+s1)/M;
+ if (L==1)
+ nv=(i+v)%P;
+ else
+ nv=(s1%M+v)%P;
+ d_NS[s*d_I+i] = ns1*P+nv;
+ d_OS[s*d_I+i] = i*d_S+s;
+ }
+ }
+
+ generate_PS_PI();
+ generate_TM();
+}
+
+
+
+
+
+
+
+
+
+
+//######################################################################
+//# Automatically generate an FSM specification describing the
//# the joint trellis of fsm1 and fsm2
//######################################################################
fsm::fsm(const fsm &FSM1, const fsm &FSM2)
@@ -419,7 +466,7 @@
//######################################################################
-//# Write trellis specification to a text files,
+//# Write trellis specification to a text file,
//# in the same format used when reading FSM files
//######################################################################
void fsm::write_fsm_txt(std::string filename)
Modified: gnuradio/trunk/gr-trellis/src/lib/fsm.h
===================================================================
--- gnuradio/trunk/gr-trellis/src/lib/fsm.h 2009-02-25 23:38:17 UTC (rev
10513)
+++ gnuradio/trunk/gr-trellis/src/lib/fsm.h 2009-02-26 00:29:51 UTC (rev
10514)
@@ -50,6 +50,7 @@
fsm(const char *name);
fsm(int k, int n, const std::vector<int> &G);
fsm(int mod_size, int ch_length);
+ fsm(int P, int M, int L);
fsm(const fsm &FSM1, const fsm &FSM2);
int I () const { return d_I; }
int S () const { return d_S; }
Modified: gnuradio/trunk/gr-trellis/src/lib/fsm.i
===================================================================
--- gnuradio/trunk/gr-trellis/src/lib/fsm.i 2009-02-25 23:38:17 UTC (rev
10513)
+++ gnuradio/trunk/gr-trellis/src/lib/fsm.i 2009-02-26 00:29:51 UTC (rev
10514)
@@ -40,6 +40,7 @@
fsm(const char *name);
fsm(int k, int n, const std::vector<int> &G);
fsm(int mod_size, int ch_length);
+ fsm(int P, int M, int L);
fsm(const fsm &FSM1, const fsm &FSM2);
int I () const { return d_I; }
int S () const { return d_S; }
Modified: gnuradio/trunk/gr-trellis/src/lib/trellis_calc_metric.cc
===================================================================
--- gnuradio/trunk/gr-trellis/src/lib/trellis_calc_metric.cc 2009-02-25
23:38:17 UTC (rev 10513)
+++ gnuradio/trunk/gr-trellis/src/lib/trellis_calc_metric.cc 2009-02-26
00:29:51 UTC (rev 10514)
@@ -31,6 +31,7 @@
{
float minm = FLT_MAX;
int minmi = 0;
+
switch (type){
case TRELLIS_EUCLIDEAN:
@@ -213,6 +214,7 @@
float minm = FLT_MAX;
int minmi = 0;
+
switch (type){
case TRELLIS_EUCLIDEAN:
for(int o=0;o<O;o++) {
@@ -222,6 +224,7 @@
metric[o]+=s.real()*s.real()+s.imag()*s.imag();
}
}
+ break;
case TRELLIS_HARD_SYMBOL:
for(int o=0;o<O;o++) {
metric[o]=0.0;
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Commit-gnuradio] r10514 - in gnuradio/trunk/gr-trellis/src: examples lib,
anastas <=