From 554836d1f3d42d46f19682cad4fc498c7c121276 Mon Sep 17 00:00:00 2001
From: Antoine Grenier <antoine.grenier@tuni.fi>
Date: Wed, 12 Apr 2023 11:12:01 +0300
Subject: [PATCH] [Convert Channel class to C implementation #44] Continue
 acquisition convertion, not working.

---
 core/c_functions/channel_l1ca.c | 97 ++++++++++++++++++++++++++++-----
 core/c_functions/fft8g.h        | 33 +++++++++++
 2 files changed, 115 insertions(+), 15 deletions(-)

diff --git a/core/c_functions/channel_l1ca.c b/core/c_functions/channel_l1ca.c
index 57322e3..1efdce6 100644
--- a/core/c_functions/channel_l1ca.c
+++ b/core/c_functions/channel_l1ca.c
@@ -2,6 +2,8 @@
 #include <stdlib.h>
 #include <complex.h>
 
+#include "fft8g.h"
+
 /// CONSTANTS
 // GPS' definition of Pi
 #define PI 3.1415926535898
@@ -38,6 +40,7 @@ struct st_Channel_L1CA {
     // Satellite variables
     int  satelliteID;    // PRN number
     int* code;           // PRN code sequence (size 1023 for GPS L1CA)
+    complex double* codeFFT;
 
     // Numerically Controled Oscilator (NCO)
     double NCO_code;
@@ -53,7 +56,14 @@ struct st_Channel_L1CA {
     double codeCounter;
 
     // Acquisition parameters
-    
+    complex double* cohCorrelationMap;
+    double* correlationMap;
+    int acq_nbCohIntegration;
+    int acq_nbNoncohIntegration;
+    double acq_frequencyRange;
+    double acq_frequencySteps;
+    int cohIntegrationCounter;
+    int noncohIntegrationCounter;
 
     // Tracking parameters
     double* correlatorsSpacing;
@@ -106,8 +116,45 @@ void ProcessHandler(st_Channel_L1CA channelVar){
 
 // --------------------------------------------------------------------------------------------------------------------
 
-void RunAcquisition(){
+void RunAcquisition(st_Channel_L1CA channelVar){
+    
+    int nbBins = channelVar.acq_frequencyRange*2 / channelVar.acq_frequencySteps + 1;
+    complex double _correlationMap[nbBins][channelVar.size];
+
+    // Perform PCPS
+    PCPS(
+        channelVar.rfData, channelVar.codeFFT, channelVar.size, channelVar.acq_frequencyRange, 
+        channelVar.acq_frequencySteps, channelVar.samplingFrequency, _correlationMap
+    );
+
+    // Perform coherent / non-coherent integration
+    for(int i=0; i < nbBins; i++){
+        for(int j=0; j < channelVar.size; j++){
+            
+            // Coherent integration
+            channelVar.cohCorrelationMap[i][j] += _correlationMap[i][j];
+
+            if (channelVar.cohIntegrationCounter+1 == channelVar.acq_nbCohIntegration){
+                // Non-coherent integration
+                channelVar.correlationMap[i][j] += abs(channelVar.cohCorrelationMap[i][j]);
+
+                // Reset variables
+                channelVar.cohCorrelationMap[i][j] = 0.0;
+                channelVar.cohIntegrationCounter = 0;
+            }
+            else{
+                channelVar.cohIntegrationCounter++;
+            }
+        }
+    }
+
+    // Perform non-coherent integration
+    if(channelVar.cohIntegrationCounter >= channelVar.acq_nbCohIntegration){
+
+    }
+
 
+    return;
 }
 
 // --------------------------------------------------------------------------------------------------------------------
@@ -158,30 +205,50 @@ void RunTracking(st_Channel_L1CA channelVar){
 
 // --------------------------------------------------------------------------------------------------------------------
 
-void PCPS(complex double* rfData, double* code, size_t size, double freqRange, double freqstep, int cohIntegration, 
-          int nonCohIntegration, double samplesPerCode, double samplingFreq){
+void PCPS(complex double* rfData, complex double* codeFFT, size_t size, double freqRange, double freqStep, 
+          double samplingFreq, complex double* r_correlationMap){
     
-    size_t size_phasePoints = cohIntegration * samplesPerCode];
-    double* phasePoints[size_phasePoints];    
+    int nbBins = freqRange*2 / freqStep + 1;
+    double phasePoints[size];
+    complex double signal[size];  
+    double freq = 0.0;
     
     // Compute the phase points
-    for(int i=0; i < cohIntegration * samplesPerCode; i++){
-        phasePoints[i] = i * 2 * PI / samplingFreq
+    for(int i=0; i < size; i++){
+        phasePoints[i] = i * 2 * PI / samplingFreq;
     }
 
     // Compute frequency bins
-    int i = 0;
-    double freq = 0.0
-    while(freq < freqRange){
-        freq = -freqRange + freqStep * i
+    for(int i = 0; i < nbBins; i++){
+        freq = -freqRange + freqStep * i;
 
-        // Generate replica
-        for(int i=0; i <size_phasePoint< i++){
+        // Generate replica and mix
+        for(int i=0; i <size; i++){
             
+            // Generate replica
+            signal[i] = cexp(I * freq * phasePoints[i]);
+
+            // Mix with carrier
+            signal[i] *= rfData[i];
         }
-    }
 
+        // Perform FFT (in-place) 
+        fft(signal, size, 0);
 
+        // Perform correlation with code
+        // TODO check if size corrects
+        for(int i=0; i < size; i++){
+            signal[i] *= codeFFT[i];
+        }
+
+        // Perform IFFT (in-place)
+        fft(signal, size, 1);
+
+        // Save results
+        r_correlationMap[i] = *signal;
+    }
+
+    return;
 }
 
 // --------------------------------------------------------------------------------------------------------------------
diff --git a/core/c_functions/fft8g.h b/core/c_functions/fft8g.h
index 314713a..14521de 100644
--- a/core/c_functions/fft8g.h
+++ b/core/c_functions/fft8g.h
@@ -1189,3 +1189,36 @@ void rftbsub(int n, double *a, int nc, double *c)
     }
     a[m + 1] = -a[m + 1];
 }
+
+/// ===================================================================================================================
+// Added by Hans
+/*
+@brief Pad an input signal to a power of two length and perform the FFT/IFFT on it.
+@param signal Array with a real-valued signal (will be over-written with complex values).
+@param size   Size of the array.
+@param doInv  Whether to perform forward or inverse FFT.
+@return void.
+*/
+void fft(double* signal, const size_t size, const char doInv)
+{
+    // Compute the size to pad the signal to
+    size_t paddedSize = 1;
+    while (paddedSize < size)
+        paddedSize <<= 1;
+
+    // Pad the signal with zeroes
+    double paddedSignal[paddedSize];
+    for (size_t i = 0; i < size; i++)
+        paddedSignal[i] = signal[i];
+    for (size_t i = size; i < paddedSize; i++)
+        paddedSignal[i] = signal[i];
+
+    // Perform the FFT in place
+    rdft(paddedSize, doInv == 0 ? 1 : -1, paddedSignal, ip, w);
+
+    // Return the result in the same array
+    for (size_t i = 0; i < size; i++)
+        signal[i] = paddedSignal[i];
+
+    return;
+}
\ No newline at end of file