From d1af45a588f547c351ed722c0f21688d51c8144b Mon Sep 17 00:00:00 2001
From: ha7ilm <retzlerandras@gmail.com>
Date: Sun, 30 Apr 2017 17:51:35 +0200
Subject: [PATCH] Added matched filter

---
 csdr.c    | 70 +++++++++++++++++++++++++++++++++------
 libcsdr.c | 99 +++++++++++++++++++++++++++++++++++++++++++------------
 libcsdr.h | 14 ++++++++
 3 files changed, 152 insertions(+), 31 deletions(-)

diff --git a/csdr.c b/csdr.c
index 52f1c66..0c415ee 100755
--- a/csdr.c
+++ b/csdr.c
@@ -3030,26 +3030,76 @@ int main(int argc, char *argv[])
         fwrite(zeros, sizeof(float), n_zero_samples, stdout);
         clone_(the_bufsize);
     }
-/*
-    if(!strcmp(argv[1], "matched_filter_cc")) //<matched_filter_type> <samples_per_symbol> <num_taps>
+
+    int matched_filter_which = 0;
+    if(
+            (!strcmp(argv[1], "firdes_matched_filter_f") && (matched_filter_which = 1)) ||
+            (!strcmp(argv[1], "matched_filter_cc") && (matched_filter_which = 2))
+    ) //(RRC <samples_per_symbol> <num_taps> <beta> | COSINE <samples_per_symbol>)
     {
+        if(argc<=2) return badsyntax("required parameter <matched_filter_type> is missing.");
+        matched_filter_type_t type = matched_filter_get_type_from_string(argv[2]);
         int samples_per_symbol = 0;
-        if(argc<=2) return badsyntax("required parameter <samples_per_symbol> is missing.");
-        sscanf(argv[2],"%d",&samples_per_symbol);
+
+        if(argc<=3) return badsyntax("required parameter <samples_per_symbol> is missing.");
+        sscanf(argv[3],"%d",&samples_per_symbol);
+
         int num_taps = 0;
-        if(argc<=2) return badsyntax("required parameter <num_taps> is missing.");
-        sscanf(argv[2],"%d",&num_taps);
+        if(argc<=4 && type!=MATCHED_FILTER_COSINE) 
+            return badsyntax("required parameter <num_taps> is missing.");
+        sscanf(argv[4],"%d",&num_taps);
+
+        float beta = 0;
+        if(argc<=5 && type==MATCHED_FILTER_RRC) 
+            return badsyntax("required parameter <beta> is missing.");
+        sscanf(argv[5],"%f",&beta);
+
+        if(type==MATCHED_FILTER_COSINE) num_taps = (2*samples_per_symbol)+1;
+        float* taps = (float*)malloc(sizeof(float)*num_taps);
+        switch(type)
+        {
+        case MATCHED_FILTER_RRC:
+            firdes_rrc_f(taps, num_taps, samples_per_symbol, beta);
+            break;
+        case MATCHED_FILTER_COSINE:
+            firdes_cosine_f(taps, num_taps, samples_per_symbol);
+            break;
+        }
+
+        if(!sendbufsize(initialize_buffers())) return -2;
+
+        if(matched_filter_which==1)
+        {
+            for(int i=0;i<num_taps;i++) printf("%f ", taps[i]);
+            return 0;
+        }
+
+        for(;;)
+        {
+            FEOF_CHECK;
+            FREAD_C;
+            apply_real_fir_cc((complexf*)input_buffer, (complexf*)output_buffer, the_bufsize, taps, num_taps);
+            FWRITE_C;
+            TRY_YIELD;
+        }
+    }
+
+    if(!strcmp(argv[1], "generic_slicer_f_u8")) //<n_symbols>
+    {
+        int n_symbols = 0;
+        if(argc<=2) return badsyntax("required parameter <n_symbols> is missing.");
+        sscanf(argv[2],"%d",&n_symbols);
         if(!sendbufsize(initialize_buffers())) return -2;
         for(;;)
         {
             FEOF_CHECK;
-            FREAD_C; 
-
-            FWRITE_C; 
+            if(!FREAD_R) break;
+            generic_slicer_f_u8(input_buffer, (unsigned char*)output_buffer, the_bufsize, n_symbols);
+            FWRITE_U8;
             TRY_YIELD;
         }
+        return 0;
     }
-*/
 
     if(!strcmp(argv[1],"none"))
     {
diff --git a/libcsdr.c b/libcsdr.c
index f2ea3e8..89e4b10 100755
--- a/libcsdr.c
+++ b/libcsdr.c
@@ -1721,6 +1721,33 @@ void serial_line_decoder_f_u8(serial_line_t* s, float* input, unsigned char* out
     DEBUG_SERIAL_LINE_DECODER && fprintf(stderr, "sld: >> output_size = %d  (+%d)\n", s->output_size, s->input_used);
 }
 
+void generic_slicer_f_u8(float* input, unsigned char* output, int input_size, int n_symbols)
+{
+    float symbol_distance = 2.0/(n_symbols+1);
+    for(int i=0;i<input_size;i++) 
+        for(int j=0;j<n_symbols;i++)
+        {
+            float symbol_center = -1+j*symbol_distance;
+            float symbol_low_limit = symbol_center-(symbol_distance/2);
+            float symbol_high_limit = symbol_center+(symbol_distance/2);
+            if(j==0)
+            {
+                if(input[i]<symbol_high_limit) output[i]=j;
+                break;
+            }
+            else if (j==n_symbols-1)
+            {
+                if(input[i]>=symbol_low_limit) output[i]=j;
+                break;
+            }
+            else 
+            {
+                if(input[i]>=symbol_low_limit && input[i]<symbol_high_limit) output[i]=j;
+                break;
+            }
+        }
+}
+
 void binary_slicer_f_u8(float* input, unsigned char* output, int input_size)
 {
     for(int i=0;i<input_size;i++) output[i] = input[i] > 0;
@@ -2169,6 +2196,24 @@ int apply_fir_cc(complexf* input, complexf* output, int input_size, complexf* ta
     return i;
 }
 
+
+int apply_real_fir_cc(complexf* input, complexf* output, int input_size, float* taps, int taps_length)
+{
+    int i;
+    for(i=0; i<input_size-taps_length+1; i++)
+    {
+        float acci = 0, accq = 0;
+        for(int ti=0;ti<taps_length;ti++)
+        {
+            acci += iof(input,i)*taps[ti];
+            accq += qof(input,i)*taps[ti];
+        }
+        iof(output,i)=acci;
+        qof(output,i)=accq;
+    }
+    return i;
+}
+
 float normalized_timing_variance_u32_f(unsigned* input, float* temp, int input_size, int samples_per_symbol, int initial_sample_offset, int debug_print)
 {
     float *ndiff_rad = temp;
@@ -2310,32 +2355,44 @@ void get_random_gaussian_samples_c(complexf* output, int output_size, FILE* stat
     }
 }
 
-/*
-void get_awgn_samples_c(complexf* output, int output_size, FILE* status)
-{
-    int* pioutput = (int*)output;
-    int cnt = 0;
-    for(int i=0;i<output_size;i++)
-    {
-        do
-        {
-            fread(pioutput+2*i, sizeof(complexf), 2, status);
-            iof(output, i)=((float)pioutput[2*i])/((float)INT_MAX); 
-            qof(output, i)=((float)pioutput[2*i+1])/((float)INT_MAX); 
-        }
-        while(sqrt(iof(output, i)*iof(output, i)+qof(output, i)*qof(output, i))>1);
-        iof(output, i)=(1/0.82)*iof(output, i);
-        qof(output, i)=(1/0.82)*qof(output, i);
-    }
-}
-*/
-
-
 int deinit_get_random_samples_f(FILE* status)
 {
     return fclose(status);
 }
 
+int firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol)
+{
+    //needs a taps_length 2 × samples_per_symbol + 1
+    int middle_i=taps_length/2;
+    for(int i=0;i<samples_per_symbol;i++) taps[middle_i+i]=taps[middle_i-i]=(1+cos(PI*i/(float)samples_per_symbol))/2;
+}
+
+int firdes_rrc_f(float* taps, int taps_length, int samples_per_symbol, float beta)
+{
+    //needs an odd taps_length
+    int middle_i=taps_length/2;
+    taps[middle_i]=(1/samples_per_symbol)*(1+beta*(4/PI-1));
+    for(int i=1;i<taps_length/2;i++) 
+    {
+        if(i==samples_per_symbol/(4*beta)) 
+            taps[middle_i+i]=taps[middle_i-i]=(beta/(samples_per_symbol*sqrt(2)))*((1+(2/PI))*sin(PI/(4*beta))+(1-(2/PI)*cos(PI/(4*beta))));
+        else
+            taps[middle_i+i]=taps[middle_i-i]=(1/samples_per_symbol)*
+                (sin(PI*(i/(float)samples_per_symbol)*(1-beta)) + 4*beta*(i/(float)samples_per_symbol)*cos(PI*(i/(float)samples_per_symbol)*(1+beta)))/
+                (PI*(i/(float)samples_per_symbol)*(1-powf(4*beta*(i/(float)samples_per_symbol),2)));
+    }
+}
+
+#define MMATCHEDFILT_GAS(NAME) \
+    if(!strcmp( #NAME , input )) return MATCHED_FILTER_ ## NAME;
+
+matched_filter_type_t matched_filter_get_type_from_string(char* input)
+{
+    MMATCHEDFILT_GAS(RRC);
+    MMATCHEDFILT_GAS(COSINE);
+    return MATCHED_FILTER_DEFAULT;
+}
+
 float* add_ff(float* input1, float* input2, float* output, int input_size)
 {
     for(int i=0;i<input_size;i++) output[i]=input1[i]+input2[i];
diff --git a/libcsdr.h b/libcsdr.h
index bd25bac..3aa5ecc 100644
--- a/libcsdr.h
+++ b/libcsdr.h
@@ -370,3 +370,17 @@ int deinit_get_random_samples_f(FILE* status);
 float* add_ff(float* input1, float* input2, float* output, int input_size);
 float total_logpower_cf(complexf* input, int input_size);
 float normalized_timing_variance_u32_f(unsigned* input, float* temp, int input_size, int samples_per_symbol, int initial_sample_offset, int debug_print);
+
+typedef enum matched_filter_type_e
+{
+    MATCHED_FILTER_RRC, 
+    MATCHED_FILTER_COSINE 
+} matched_filter_type_t;
+
+#define MATCHED_FILTER_DEFAULT MATCHED_FILTER_RRC
+
+int firdes_cosine_f(float* taps, int taps_length, int samples_per_symbol);
+int firdes_rrc_f(float* taps, int taps_length, int samples_per_symbol, float beta);
+matched_filter_type_t matched_filter_get_type_from_string(char* input);
+int apply_real_fir_cc(complexf* input, complexf* output, int input_size, float* taps, int taps_length);
+void generic_slicer_f_u8(float* input, unsigned char* output, int input_size, int n_symbols);