From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Mon, 23 Apr 2012 01:34:01 +0000 (+0000)
Subject: rx downconverterand rate P rx filtering working in C
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=58891ed4a4a9db25aeb486a7c386ec5c403f35d9;p=freetel-svn-tracking.git

rx downconverterand rate P rx filtering working in C

git-svn-id: https://svn.code.sf.net/p/freetel/code@375 01035d8c-6547-0410-b346-abe4f91aad63
---

diff --git a/codec2-dev/octave/README_fdmdv.txt b/codec2-dev/octave/README_fdmdv.txt
index a22b65a5..541e4a81 100644
--- a/codec2-dev/octave/README_fdmdv.txt
+++ b/codec2-dev/octave/README_fdmdv.txt
@@ -7,12 +7,6 @@ sox -r 8000 -s -2 mod_dqpsk.raw -s -2 mod_dqpsk_8008hz.raw rate -h 8008
 TODO
 
 [X] Get file based mod and demod working again
-[ ] timing wraps around
-    + what is affect of bouncing back and forth over boundary?
-    + could mean we are going back and forth a symbol
-    + do we need to logic to lose or gain a frame?
-    + think so, add or remove samples, or a whole frame
-[ ] demod outputs ber (maybe after settling time)
 [ ] try interfering sine wave
     + maybe swept
     + does modem fall over?
@@ -20,9 +14,13 @@ TODO
     + make sure all estimators keep working
 [ ] test rx level sensitivity, i.e. 0 to 20dB attenuation
 [ ] try to run from shell script
-[ ] run a few tests
-[ ] start coding in C and repeat tests
 [ ] arb bit stream input to Octave mod and demod
+[ ] C port and UT framework
+    [ ] doumnent how to use
+[ ] document use of fdmdv_ut and fdmdv_demod + PathSim
+[ ] block diagram
+[ ] blog posts(s)
+[ ] Codec 2 web page update
 
 Tests
 
diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m
index 7a8f711d..afa503e6 100644
--- a/codec2-dev/octave/fdmdv.m
+++ b/codec2-dev/octave/fdmdv.m
@@ -172,7 +172,7 @@ function tx_fdm = fdm_upconvert(tx_filt)
 endfunction
 
 
-% Frequency shift each modem carrier down to Nc baseband signals
+% Frequency shift each modem carrier down to Nc+1 baseband signals
 
 function rx_baseband = fdm_downconvert(rx_fdm, nin)
   global Fs;
diff --git a/codec2-dev/octave/tfdmdv.m b/codec2-dev/octave/tfdmdv.m
index fd6e9e75..1ba7a794 100644
--- a/codec2-dev/octave/tfdmdv.m
+++ b/codec2-dev/octave/tfdmdv.m
@@ -31,6 +31,9 @@ pilot_lpf1_log = [];
 pilot_lpf2_log = [];
 S1_log = [];
 S2_log = [];
+foff_log = [];
+rx_baseband_log = [];
+rx_filt_log = [];
 
 for f=1:frames
 
@@ -53,6 +56,7 @@ for f=1:frames
   [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, M);
 
   [foff S1 S2] = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, M);
+  foff_log = [foff_log foff];
 
   pilot_baseband1_log = [pilot_baseband1_log pilot_baseband1];
   pilot_baseband2_log = [pilot_baseband2_log pilot_baseband2];
@@ -61,6 +65,11 @@ for f=1:frames
   S1_log  = [S1_log S1];
   S2_log  = [S2_log S2];
 
+  rx_baseband = fdm_downconvert(rx_fdm, M);
+  rx_baseband_log = [rx_baseband_log rx_baseband];
+  rx_filt = rx_filter(rx_baseband, M);
+  rx_filt_log = [rx_filt_log rx_filt];
+
 end
 
 % Compare to the output from the C version
@@ -136,6 +145,15 @@ plot_sig_and_error(7, 212, imag(S1_log), imag(S1_log - S1_log_c), 'S1 imag' )
 plot_sig_and_error(8, 211, real(S2_log), real(S2_log - S2_log_c), 'S2 real' )
 plot_sig_and_error(8, 212, imag(S2_log), imag(S2_log - S2_log_c), 'S2 imag' )
 
+plot_sig_and_error(9, 211, real(foff_log), real(foff_log - foff_log_c), 'Freq Offset' )
+
+c=10;
+plot_sig_and_error(10, 211, real(rx_baseband_log(c,:)), real(rx_baseband_log(c,:) - rx_baseband_log_c(c,:)), 'Rx baseband real' )
+plot_sig_and_error(10, 212, imag(rx_baseband_log(c,:)), imag(rx_baseband_log(c,:) - rx_baseband_log_c(c,:)), 'Rx baseband imag' )
+
+plot_sig_and_error(11, 211, real(rx_filt_log(c,:)), real(rx_filt_log(c,:) - rx_filt_log_c(c,:)), 'Rx filt real' )
+plot_sig_and_error(11, 212, imag(rx_filt_log(c,:)), imag(rx_filt_log(c,:) - rx_filt_log_c(c,:)), 'Rx filt imag' )
+
 % ---------------------------------------------------------------------------------------
 % AUTOMATED CHECKS ------------------------------------------
 % ---------------------------------------------------------------------------------------
@@ -169,5 +187,8 @@ check(pilot_baseband1_log, pilot_baseband1_log_c, 'pilot lpf1');
 check(pilot_baseband2_log, pilot_baseband2_log_c, 'pilot lpf2');
 check(S1_log, S1_log_c, 'S1');
 check(S2_log, S2_log_c, 'S2');
+check(foff_log, foff_log_c, 'rx_est_freq_offset');
+check(rx_baseband_log, rx_baseband_log_c, 'fdm_downconvert');
+check(rx_filt_log, rx_filt_log_c, 'fdm_downconvert');
 
 printf("\npasses: %d fails: %d\n", passes, fails);
diff --git a/codec2-dev/src/fdmdv.c b/codec2-dev/src/fdmdv.c
index add05ace..13f6c20f 100644
--- a/codec2-dev/src/fdmdv.c
+++ b/codec2-dev/src/fdmdv.c
@@ -156,6 +156,8 @@ struct FDMDV *fdmdv_create(void)
 	for(k=0; k<NFILTER; k++) {
 	    f->tx_filter_memory[c][k].real = 0.0;
 	    f->tx_filter_memory[c][k].imag = 0.0;
+	    f->rx_filter_memory[c][k].real = 0.0;
+	    f->rx_filter_memory[c][k].imag = 0.0;
 	}
 
 	/* Spread initial FDM carrier phase out as far as possible.
@@ -164,6 +166,8 @@ struct FDMDV *fdmdv_create(void)
 	
 	f->phase_tx[c].real = cos(2.0*PI*c/(NC+1));
  	f->phase_tx[c].imag = sin(2.0*PI*c/(NC+1));
+	f->phase_rx[c].real = 1.0;
+ 	f->phase_rx[c].imag = 0.0;
 
    }
     
@@ -546,7 +550,7 @@ void lpf_peak_pick(float *foff, float *max, COMP pilot_baseband[], COMP pilot_lp
     imax = 0.0;
     ix = 0;
     for(i=0; i<MPILOTFFT; i++) {
-	mag = sqrt(S[i].real*S[i].real + S[i].imag*S[i].imag);
+	mag = S[i].real*S[i].real + S[i].imag*S[i].imag;
 	if (mag > imax) {
 	    imax = mag;
 	    ix = i;
@@ -600,7 +604,7 @@ float rx_est_freq_offset(struct FDMDV *f, float rx_fdm[], int nin)
 
     /*
       Down convert latest M samples of pilot by multiplying by ideal
-      BPSK pilot signal we have generated locally.  This peak of the
+      BPSK pilot signal we have generated locally.  The peak of the
       resulting signal is sensitive to the time shift between the
       received and local version of the pilot, so we do it twice at
       different time shifts and choose the maximum.
@@ -626,3 +630,99 @@ float rx_est_freq_offset(struct FDMDV *f, float rx_fdm[], int nin)
 	
     return foff;
 }
+
+/*---------------------------------------------------------------------------*\
+                                                       
+  FUNCTION....: fdm_downconvert()	     
+  AUTHOR......: David Rowe			      
+  DATE CREATED: 22/4/2012
+
+  Frequency shift each modem carrier down to Nc+1 baseband signals.
+
+\*---------------------------------------------------------------------------*/
+
+void fdm_downconvert(COMP rx_baseband[NC+1][M+M/P], COMP rx_fdm[], COMP phase_rx[], COMP freq[], int nin)
+{
+    int  i,c;
+    COMP pilot;
+
+    /* maximum number of input samples to demod */
+
+    assert(nin < (M+M/P));
+
+    /* Nc/2 tones below centre freq */
+  
+    for (c=0; c<NC/2; c++) 
+	for (i=0; i<nin; i++) {
+	    phase_rx[c] = cmult(phase_rx[c], freq[c]);
+	    rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));
+	}
+
+    /* Nc/2 tones above centre freq */
+
+    for (c=NC/2; c<NC; c++) 
+	for (i=0; i<M; i++) {
+	    phase_rx[c] = cmult(phase_rx[c], freq[c]);
+	    rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));
+	}
+
+    /* add centre pilot tone  */
+
+    c = NC;
+    for (i=0; i<M; i++) {
+	phase_rx[c] = cmult(phase_rx[c],  freq[c]);
+	rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));
+    }
+
+}
+
+/*---------------------------------------------------------------------------*\
+                                                       
+  FUNCTION....: rx_filter()	     
+  AUTHOR......: David Rowe			      
+  DATE CREATED: 22/4/2012
+
+  Receive filter each baseband signal at oversample rate P.  Filtering at
+  rate P lowers CPU compared to rate M.
+
+  Depending on the number of input samples to the demod nin, we
+  produce P-1, P (usually), or P+1 filtered samples at rate P.  nin is
+  occasionally adjusted to compensate for timing slips due to
+  different tx and rx sample clocks.
+
+\*---------------------------------------------------------------------------*/
+
+void rx_filter(COMP rx_filt[NC+1][P+1], COMP rx_baseband[NC+1][M+M/P], COMP rx_filter_memory[NC+1][NFILTER], int nin)
+{
+    int c, i,j,k,l;
+    int n=M/P;
+
+    /* rx filter each symbol, generate P filtered output samples for
+       each symbol.  Note we keep filter memory at rate M, it's just
+       the filter output at rate P */
+
+    for(i=0, j=0; i<nin; i+=n,j++) {
+
+	/* latest input sample */
+	
+	for(c=0; c<NC+1; c++)
+	    for(k=NFILTER-n,l=i; k<NFILTER; k++,l++)	
+		rx_filter_memory[c][k] = rx_baseband[c][l];
+	
+	/* convolution (filtering) */
+
+	for(c=0; c<NC+1; c++) {
+	    rx_filt[c][j].real = 0.0; rx_filt[c][j].imag = 0.0;
+	    for(k=0; k<NFILTER; k++) 
+		rx_filt[c][j] = cadd(rx_filt[c][j], fcmult(gt_alpha5_root[k], rx_filter_memory[c][k]));
+	}
+
+	/* make room for next input sample */
+	
+	for(c=0; c<NC+1; c++)
+	    for(k=0,l=n; k<NFILTER-n; k++,l++)	
+		rx_filter_memory[c][k] = rx_filter_memory[c][l];
+    }
+
+    assert(j <= (P+1)); /* check for any over runs */
+}
diff --git a/codec2-dev/src/fdmdv_internal.h b/codec2-dev/src/fdmdv_internal.h
index 47bee78d..d2ffadee 100644
--- a/codec2-dev/src/fdmdv_internal.h
+++ b/codec2-dev/src/fdmdv_internal.h
@@ -71,6 +71,7 @@
 
 struct FDMDV {
     int  current_test_bit;
+
     int  tx_pilot_bit;
     COMP prev_tx_symbols[NC+1];
     COMP tx_filter_memory[NC+1][NFILTER];
@@ -87,6 +88,10 @@ struct FDMDV {
     COMP pilot_lpf2[NPILOTLPF];
     COMP S1[MPILOTFFT];
     COMP S2[MPILOTFFT];
+
+    COMP phase_rx[NC+1];
+    COMP rx_filter_memory[NC+1][NFILTER];
+
 };
 
 /*---------------------------------------------------------------------------*\
@@ -102,5 +107,7 @@ void generate_pilot_fdm(COMP *pilot_fdm, int *bit, float *symbol, float *filter_
 void generate_pilot_lut(COMP pilot_lut[], COMP *pilot_freq);
 float rx_est_freq_offset(struct FDMDV *f, float rx_fdm[], int nin);
 void lpf_peak_pick(float *foff, float *max, COMP pilot_baseband[], COMP pilot_lpf[], COMP S[], int nin);
+void fdm_downconvert(COMP rx_baseband[NC+1][M+M/P], COMP rx_fdm[], COMP phase_rx[], COMP freq[], int nin);
+void rx_filter(COMP rx_filt[NC+1][P+1], COMP rx_baseband[NC+1][M+M/P], COMP rx_filter_memory[NC+1][NFILTER], int nin);
 
 #endif
diff --git a/codec2-dev/unittest/tfdmdv.c b/codec2-dev/unittest/tfdmdv.c
index ae11da2d..4c82927b 100644
--- a/codec2-dev/unittest/tfdmdv.c
+++ b/codec2-dev/unittest/tfdmdv.c
@@ -4,10 +4,9 @@
   AUTHOR......: David Rowe  
   DATE CREATED: April 16 2012
                                                                              
-  Unit tests for FDMDV modem.  Combination of unit tests perfromed
-  entirely by this program and comparisons with reference Octave
-  version of the modem that require running an Octave script
-  ../octave/tfdmdv.m.
+  Tests for the C version of the FDMDV modem.  This program outputs a
+  file of Octave vectors that are loaded and compared to the
+  Octave version of thr modem by the Octave script tfmddv.m
                                                                              
 \*---------------------------------------------------------------------------*/
 
@@ -33,6 +32,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 
 #include "fdmdv_internal.h"
 #include "fdmdv.h"
@@ -40,7 +40,8 @@
 #define FRAMES 25
 
 void octave_save_int(FILE *f, char name[], int data[], int rows, int cols);
-void octave_save_complex(FILE *f, char name[], COMP data[], int rows, int cols);
+void octave_save_float(FILE *f, char name[], float data[], int rows, int cols);
+void octave_save_complex(FILE *f, char name[], COMP data[], int rows, int cols, int col_len);
 
 int main(int argc, char *argv[])
 {
@@ -49,8 +50,14 @@ int main(int argc, char *argv[])
     COMP          tx_symbols[(NC+1)];
     COMP          tx_baseband[(NC+1)][M];
     COMP          tx_fdm[M];
-    float         rx_fdm[M];
+    float         rx_fdm[M+M/P];
     float         foff;
+    COMP          foff_rect;
+    COMP          foff_phase_rect;
+    int           nin;
+    COMP          rx_fdm_fcorr[M+M/P];
+    COMP          rx_baseband[NC+1][M+M/P];
+    COMP          rx_filt[NC+1][P+1];
 
     int           tx_bits_log[FDMDV_BITS_PER_FRAME*FRAMES];
     COMP          tx_symbols_log[(NC+1)*FRAMES];
@@ -62,15 +69,24 @@ int main(int argc, char *argv[])
     COMP          pilot_lpf2_log[NPILOTLPF*FRAMES];
     COMP          S1_log[MPILOTFFT*FRAMES];
     COMP          S2_log[MPILOTFFT*FRAMES];
+    float         foff_log[FRAMES];
+    COMP          rx_baseband_log[(NC+1)][(M+M/P)*FRAMES];
+    int           rx_baseband_log_col_index;
+    COMP          rx_filt_log[NC+1][(P+1)*FRAMES];
+    int           rx_filt_log_col_index;
 
     FILE         *fout;
     int           f,c,i;
 
     fdmdv = fdmdv_create();
+    foff_phase_rect.real = 0.0; foff_phase_rect.imag = 0.0;
+
+    rx_baseband_log_col_index = 0;
+    rx_filt_log_col_index = 0;
 
     for(f=0; f<FRAMES; f++) {
 
-	/* modulator */
+	/* modulator -----------------------------------------*/
 
 	fdmdv_get_test_bits(fdmdv, tx_bits);
 	bits_to_dqpsk_symbols(tx_symbols, fdmdv->prev_tx_symbols, tx_bits, &fdmdv->tx_pilot_bit);
@@ -80,12 +96,29 @@ int main(int argc, char *argv[])
 
 	for(i=0; i<M; i++)
 	    rx_fdm[i] = tx_fdm[i].real;
+	nin = M;
+
+	/* demodulator ----------------------------------------*/
+
+	/* Freq offset estimation and correction */
+
+	foff = rx_est_freq_offset(fdmdv, rx_fdm, nin);
 
-	/* demodulator */
+	/* note this should be a C function with states in fdmdv */
 
-	foff = rx_est_freq_offset(fdmdv, rx_fdm, M);
- 
-	/* save log of outputs */
+	foff_rect.real = cos(2.0*PI*foff/FS);
+	foff_rect.imag = sin(2.0*PI*foff/FS);
+	for(i=0; i<nin; i++) {
+	    //foff_phase_rect = cmult(foff_phase_rect, foff_rect);
+	    //rx_fdm_fcorr[i] = fcmult(rx_fdm[i], foff_phase_rect);
+	    rx_fdm_fcorr[i].real = rx_fdm[i];
+	    rx_fdm_fcorr[i].imag = 0.0;
+	}
+	
+	fdm_downconvert(rx_baseband, rx_fdm_fcorr, fdmdv->phase_rx, fdmdv->freq, nin);
+	rx_filter(rx_filt, rx_baseband, fdmdv->rx_filter_memory, nin);
+
+	/* save log of outputs ------------------------------------------------------*/
 
 	memcpy(&tx_bits_log[FDMDV_BITS_PER_FRAME*f], tx_bits, sizeof(int)*FDMDV_BITS_PER_FRAME);
 	memcpy(&tx_symbols_log[(NC+1)*f], tx_symbols, sizeof(COMP)*(NC+1));
@@ -93,12 +126,33 @@ int main(int argc, char *argv[])
 	    for(i=0; i<M; i++)
 		tx_baseband_log[c][f*M+i] = tx_baseband[c][i]; 
 	memcpy(&tx_fdm_log[M*f], tx_fdm, sizeof(COMP)*M);
+
+	/* freq offset estimation */
+
 	memcpy(&pilot_baseband1_log[f*NPILOTBASEBAND], fdmdv->pilot_baseband1, sizeof(COMP)*NPILOTBASEBAND);
 	memcpy(&pilot_baseband2_log[f*NPILOTBASEBAND], fdmdv->pilot_baseband2, sizeof(COMP)*NPILOTBASEBAND);
 	memcpy(&pilot_lpf1_log[f*NPILOTLPF], fdmdv->pilot_lpf1, sizeof(COMP)*NPILOTLPF);
 	memcpy(&pilot_lpf2_log[f*NPILOTLPF], fdmdv->pilot_lpf2, sizeof(COMP)*NPILOTLPF);
 	memcpy(&S1_log[f*MPILOTFFT], fdmdv->S1, sizeof(COMP)*MPILOTFFT);
 	memcpy(&S2_log[f*MPILOTFFT], fdmdv->S2, sizeof(COMP)*MPILOTFFT);
+ 	memcpy(&foff_log[f], &foff, sizeof(float));
+
+	/* rx down conversion */
+
+	for(c=0; c<NC+1; c++) {
+	    for(i=0; i<nin; i++)
+		rx_baseband_log[c][rx_baseband_log_col_index + i] = rx_baseband[c][i]; 
+	}
+	rx_baseband_log_col_index += nin;
+
+	/* rx filtering */
+
+	for(c=0; c<NC+1; c++) {
+	    for(i=0; i<(P*M)/nin; i++)
+		rx_filt_log[c][rx_filt_log_col_index + i] = rx_filt[c][i]; 
+	}
+	rx_filt_log_col_index += (P*M)/nin;
+
     }
 
     /* dump logs to Octave file for evaluation by tfdmdv.m Octave script */
@@ -107,16 +161,19 @@ int main(int argc, char *argv[])
     assert(fout != NULL);
     fprintf(fout, "# Created by tfdmdv.c\n");
     octave_save_int(fout, "tx_bits_log_c", tx_bits_log, 1, FDMDV_BITS_PER_FRAME*FRAMES);
-    octave_save_complex(fout, "tx_symbols_log_c", tx_symbols_log, 1, (NC+1)*FRAMES);  
-    octave_save_complex(fout, "tx_baseband_log_c", (COMP*)tx_baseband_log, (NC+1), M*FRAMES);  
-    octave_save_complex(fout, "tx_fdm_log_c", (COMP*)tx_fdm_log, 1, M*FRAMES);  
-    octave_save_complex(fout, "pilot_lut_c", (COMP*)fdmdv->pilot_lut, 1, NPILOT_LUT);  
-    octave_save_complex(fout, "pilot_baseband1_log_c", pilot_baseband1_log, 1, NPILOTBASEBAND*FRAMES);  
-    octave_save_complex(fout, "pilot_baseband2_log_c", pilot_baseband2_log, 1, NPILOTBASEBAND*FRAMES);  
-    octave_save_complex(fout, "pilot_lpf1_log_c", pilot_lpf1_log, 1, NPILOTLPF*FRAMES);  
-    octave_save_complex(fout, "pilot_lpf2_log_c", pilot_lpf2_log, 1, NPILOTLPF*FRAMES);  
-    octave_save_complex(fout, "S1_log_c", S1_log, 1, MPILOTFFT*FRAMES);  
-    octave_save_complex(fout, "S2_log_c", S2_log, 1, MPILOTFFT*FRAMES);  
+    octave_save_complex(fout, "tx_symbols_log_c", tx_symbols_log, 1, (NC+1)*FRAMES, (NC+1)*FRAMES);  
+    octave_save_complex(fout, "tx_baseband_log_c", (COMP*)tx_baseband_log, (NC+1), M*FRAMES, M*FRAMES);  
+    octave_save_complex(fout, "tx_fdm_log_c", (COMP*)tx_fdm_log, 1, M*FRAMES, M*FRAMES);  
+    octave_save_complex(fout, "pilot_lut_c", (COMP*)fdmdv->pilot_lut, 1, NPILOT_LUT, NPILOT_LUT);  
+    octave_save_complex(fout, "pilot_baseband1_log_c", pilot_baseband1_log, 1, NPILOTBASEBAND*FRAMES, NPILOTBASEBAND*FRAMES);  
+    octave_save_complex(fout, "pilot_baseband2_log_c", pilot_baseband2_log, 1, NPILOTBASEBAND*FRAMES, NPILOTBASEBAND*FRAMES);  
+    octave_save_complex(fout, "pilot_lpf1_log_c", pilot_lpf1_log, 1, NPILOTLPF*FRAMES, NPILOTLPF*FRAMES);  
+    octave_save_complex(fout, "pilot_lpf2_log_c", pilot_lpf2_log, 1, NPILOTLPF*FRAMES, NPILOTLPF*FRAMES);  
+    octave_save_complex(fout, "S1_log_c", S1_log, 1, MPILOTFFT*FRAMES, MPILOTFFT*FRAMES);  
+    octave_save_complex(fout, "S2_log_c", S2_log, 1, MPILOTFFT*FRAMES, MPILOTFFT*FRAMES);  
+    octave_save_float(fout, "foff_log_c", foff_log, 1, FRAMES);  
+    octave_save_complex(fout, "rx_baseband_log_c", (COMP*)rx_baseband_log, (NC+1), rx_baseband_log_col_index, (M+M/P)*FRAMES);  
+    octave_save_complex(fout, "rx_filt_log_c", (COMP*)rx_filt_log, (NC+1), rx_filt_log_col_index, (P+1)*FRAMES);  
     fclose(fout);
 
     codec2_destroy(fdmdv);
@@ -142,18 +199,35 @@ void octave_save_int(FILE *f, char name[], int data[], int rows, int cols)
     fprintf(f, "\n\n");
 }
 
-void octave_save_complex(FILE *f, char name[], COMP data[], int rows, int cols)
+void octave_save_float(FILE *f, char name[], float data[], int rows, int cols)
 {
     int r,c;
 
     fprintf(f, "# name: %s\n", name);
-    fprintf(f, "# type: complex matrix\n");
+    fprintf(f, "# type: matrix\n");
     fprintf(f, "# rows: %d\n", rows);
     fprintf(f, "# columns: %d\n", cols);
     
     for(r=0; r<rows; r++) {
 	for(c=0; c<cols; c++)
-	    fprintf(f, " (%f,%f)", data[r*cols+c].real, data[r*cols+c].imag);
+	    fprintf(f, " %f", data[r*cols+c]);
+	fprintf(f, "\n");
+    }
+
+    fprintf(f, "\n\n");
+}
+void octave_save_complex(FILE *f, char name[], COMP data[], int rows, int cols, int col_len)
+{
+    int r,c;
+
+    fprintf(f, "# name: %s\n", name);
+    fprintf(f, "# type: complex matrix\n");
+    fprintf(f, "# rows: %d\n", rows);
+    fprintf(f, "# columns: %d\n", cols);
+    printf("rows %d cols %d col_len %d\n", rows, cols, col_len);
+    for(r=0; r<rows; r++) {
+	for(c=0; c<cols; c++)
+	    fprintf(f, " (%f,%f)", data[r*col_len+c].real, data[r*col_len+c].imag);
 	fprintf(f, "\n");
     }