API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Nov 20th, 2017
90
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 46.09 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <fftw3.h>
  3. #include <cmath>
  4. #include "print_utilities.h"
  5. #include <string>
  6. #include <stdio.h>
  7. #include <time.h>
  8. #include <math.h>
  9.  
  10. #define             REAL 0
  11. #define             IMAG 1
  12. #define             PI 3.14159265
  13. using namespace     std;
  14.  
  15. struct input_params{
  16.         double      f0;
  17.         double      chirp_rate;
  18.         double      Fs;
  19.         double      Tmod;
  20.         double      va;
  21.         double*     xa;
  22.         double*     xa_offset;
  23.         int         Nfft_range;
  24.         double*     range_ax;
  25.         double*     azimuth_ax;
  26.         int         azimuth_ax_size;
  27.         double*     Win_range;
  28.         double*     Win_azimuth;
  29.         int         rg_min_index;
  30.         int         rg_max_index;
  31. };
  32.  
  33.  
  34. void stampa_matrice_complessa(fftw_complex** matrice, int n, int m){
  35.     for(int i = 0; i < n; i++)
  36.         for(int j = 0; j < m; j++){
  37.             char spazio = (!((j+1)%m))?'\n':'\t';
  38.             cout << matrice[i][j][REAL] << " " << matrice[i][j][IMAG] << spazio;
  39.         }
  40. }
  41.  
  42. void stampa_matrice_reale(double** matrice, int n, int m){
  43.     for(int i = 0; i < n; i++)
  44.         for(int j = 0; j < m; j++){
  45.             char spazio = (!((j+1)%m))?'\n':'\t';
  46.             cout << matrice[i][j] << spazio;
  47.         }
  48. }
  49.  
  50. void inizializza_vettore(double* v, int n){
  51.     for(int i = 0; i < n; i++)
  52.         v[i] = 0;
  53. }
  54.  
  55. double* linspace(double start, double end, int n){
  56.     double* vect = new double[n];
  57.     double step = (end-start)/(n-1);
  58.     for(int i = 0; i < n; i++)
  59.         vect[i] = (step*i)+start;
  60.     return vect;
  61. }
  62.  
  63. void stampa_vettore(double* v, int size){
  64.     for(int i = 0; i < size; i++)
  65.         cout << v[i] << "\t";
  66.     cout << "\n";
  67. }
  68.  
  69. double** create_real_matrix(int n, int m){
  70.     double** matrix = new double*[n];
  71.     for(int i = 0; i < n; i++){
  72.         matrix[i] = new double[m];
  73.         for(int j = 0; j < m; j++)
  74.             matrix[i][j] = 0;
  75.     }
  76.  
  77.     return matrix;
  78. }
  79.  
  80. double* create_real_vector(int n){
  81.     double* v = new double[n];  
  82.     for(int i = 0; i  < n; i++)
  83.         v[i] = 0;
  84.     return v;
  85. }
  86.  
  87. void stampa_tensore(double*** tensore, int n, int m, int p){
  88.     for(int i = 0; i < n; i++)
  89.         for(int j = 0; j < m; j++)
  90.             for(int k = 0; k < p; k++)
  91.                 cout << tensore[i][j][k] << endl;
  92. }
  93.  
  94. fftw_complex*** zeros_3D_complex(int n, int m, int p){
  95.     fftw_complex*** tensore = new fftw_complex**[n];
  96.     for(int i = 0; i < n; i++){
  97.         tensore[i] = new fftw_complex*[m];
  98.         for(int j = 0; j < m; j++){
  99.             tensore[i][j] = new fftw_complex[p];
  100.             for(int k = 0; k < p; k++){
  101.                 tensore[i][j][k][REAL] = 0;
  102.                 tensore[i][j][k][IMAG] = 0;
  103.             }
  104.         }
  105.     }
  106.     return tensore;
  107. }
  108.  
  109. //fftw_complex*** zeros_3D_complex(int L, int M, int N){
  110. //  fftw_complex*** data;
  111. //
  112. //  data = (fftw_complex***) fftw_malloc(sizeof(fftw_complex) * L*M*N);
  113. //
  114. //  for(int l = 0; l < L; l++)
  115. //      for(int m = 0; m < M; m++)
  116. //          for(int n = 0; n < N; n++){
  117. //              data[l][m][n][REAL] = 0;
  118. //              data[l][m][n][IMAG] = 0;
  119. //          }
  120. //  return data;
  121. //}
  122.  
  123. double*** real_zeros(int n, int m, int p){
  124.     double*** tensore = new double**[n];
  125.     for(int i = 0; i < n; i++){
  126.         tensore[i] = new double*[m];
  127.         for(int j = 0; j < m; j++){
  128.             tensore[i][j] = new double[p];
  129.             for(int k = 0; k < p; k++)
  130.                 tensore[i][j][k] = 0;
  131.         }
  132.     }
  133.     return tensore;
  134. }
  135.  
  136. fftw_complex** create_complex_matrix(int rows, int columns){
  137.      
  138.     fftw_complex** matrix = new fftw_complex*[rows];
  139.     for(int i=0; i < rows; i++){
  140.         matrix[i] = new fftw_complex[columns];
  141.         for(int j = 0; j < columns; j++){
  142.             matrix[i][j][REAL] = 0;
  143.             matrix[i][j][IMAG] = 0;
  144.         }
  145.     }      
  146.  
  147.     return matrix;
  148.  
  149. }
  150.  
  151. void fft(fftw_complex* in, fftw_complex* out, int n){
  152.  
  153.     fftw_plan plan = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
  154.     fftw_execute(plan);
  155.     fftw_destroy_plan(plan);
  156.     fftw_cleanup();
  157.  
  158. }
  159.  
  160. void ifft(fftw_complex* in, fftw_complex* out, double n){
  161.  
  162.     fftw_plan plan = fftw_plan_dft_1d(n, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
  163.     fftw_execute(plan);
  164.     fftw_destroy_plan(plan);
  165.     fftw_cleanup();
  166.  
  167.     for(int i = 0; i < n; i++){
  168.         out[i][REAL] /= n;
  169.         out[i][IMAG] /= n;
  170.     }
  171.  
  172. }
  173.  
  174. void fftshift(fftw_complex* in, int n){
  175.      
  176.     fftw_complex support[n];
  177.     int mid = (n/2);
  178.     for(int i = mid; i < n; i++){
  179.         support[i][REAL] = in[i-mid][REAL];
  180.         support[i][IMAG] = in[i-mid][IMAG];
  181.     }
  182.     for(int i = 0; i < mid; i++){
  183.         if(n%2 == 1){
  184.             support[i][REAL] = in[i+mid+1][REAL];
  185.             support[i][IMAG] = in[i+mid+1][IMAG];
  186.         }
  187.         else{
  188.             support[i][REAL] = in[i+mid][REAL];
  189.             support[i][IMAG] = in[i+mid][IMAG];
  190.         }
  191.     }
  192.     for(int i = 0; i < n; i++){
  193.         in[i][REAL] = support[i][REAL];
  194.         in[i][IMAG] = support[i][IMAG];
  195.     }
  196.  
  197. }
  198.  
  199. void ifftshift(fftw_complex* in, int n){
  200.      
  201.     fftw_complex support[n];
  202.     int mid = (n/2);
  203.     for(int i = mid; i < n; i++){
  204.         support[i-mid][REAL] = in[i][REAL];
  205.         support[i-mid][IMAG] = in[i][IMAG];
  206.     }
  207.     for(int i = 0; i < mid; i++){
  208.         if(n%2){
  209.             support[i+mid+1][REAL] = in[i][REAL];
  210.             support[i+mid+1][IMAG] = in[i][IMAG];
  211.         }
  212.         else{
  213.             support[i+mid][REAL] = in[i][REAL];
  214.             support[i+mid][IMAG] = in[i][IMAG];
  215.         }
  216.     }
  217.     for(int i = 0; i < n; i++){
  218.         in[i][REAL] = support[i][REAL];
  219.         in[i][IMAG] = support[i][IMAG];
  220.     }
  221.  
  222. }
  223.  
  224. void copy_vector(fftw_complex* source, fftw_complex* destination, int n){
  225.      
  226.     for(int i = 0; i < n; i++){
  227.         destination[i][REAL] = source[i][REAL];
  228.         destination[i][IMAG] = source[i][IMAG];
  229.     }
  230.  
  231. }
  232.  
  233. void repmat_scalar_vector(double rep, double* destination, int n){
  234.     for(int i = 0; i < n; i++)
  235.         destination[i] = rep;
  236. };
  237.  
  238. void repmat_colvector_matrix(double* rep, double** destination, int n, int m){
  239.     for(int i = 0; i < m; i++)
  240.         for(int j = 0; j < n; j++)
  241.             destination[j][i] = rep[j];
  242. }
  243.  
  244. void repmat_rowvector_matrix(double* rep, double** destination, int n, int m){
  245.     for(int i = 0; i < n; i++)
  246.         for(int j = 0; j < m; j++)
  247.             destination[i][j] = rep[j];
  248. }
  249.  
  250. void repmat_rowvector_tensor(double* rep, double*** tensor, int n, int m, int p){
  251.     for(int i = 0; i < n; i++)
  252.         for(int j = 0; j < m; j++)
  253.             for(int k = 0; k < p; k++)
  254.                 tensor[i][j][k] = rep[j];
  255. }
  256.  
  257. void scalar_times_vector(double* source, double* target, double scalar, int n){
  258.     for(int i = 0; i < n; i++)
  259.         target[i] = scalar*source[i];
  260. }
  261.  
  262. void vector_time_matrix_scalar(double* vector, double** matrix, int n, int m){
  263.     for(int i = 0; i < n; i++)
  264.         for(int j = 0; j < m; j++)
  265.             matrix[i][j] *= vector[i];
  266. }
  267.  
  268. void vector_time_complex_matrix_scalar(double* vector, fftw_complex** matrix, int n, int m){
  269.     for(int i = 0; i < n; i++)
  270.         for(int j = 0; j < m; j++){
  271.             matrix[i][j][REAL] *= vector[i];
  272.             matrix[i][j][IMAG] *= vector[i];
  273.         }
  274. }
  275.  
  276. void invvector_time_matrix_scalar(double* vector, double** matrix, int n, int m){
  277.     for(int i = 0; i < n; i++)
  278.         for(int j = 0; j < m; j++)
  279.             matrix[i][j] *= 1/vector[i];
  280. }
  281.  
  282. double complex_exponential_real_part(double phi){
  283.     return cos(phi);
  284. }
  285.  
  286. double complex_exponential_imag_part(double phi){
  287.     return sin(phi);
  288. }
  289.  
  290. void matrix_scalar_square(double** matrix, int n, int m){
  291.     for(int i = 0; i < n; i++)
  292.         for(int j = 0; j < m; j++)
  293.             matrix[i][j] *= matrix[i][j];
  294. }
  295.  
  296. void matrix_scalar_cube(double** matrix, int n, int m){
  297.     for(int i = 0; i < n; i++)
  298.         for(int j = 0; j < m; j++)
  299.             matrix[i][j] *= (matrix[i][j]*matrix[i][j]);
  300. }
  301.  
  302. void vector_scalar_square(double* target, int n){
  303.     for(int i = 0; i < n; i++)
  304.         target[i] *= target[i];
  305. }
  306.  
  307. void matrix_times_scalar(double** matrix, double scalar, int n, int m){
  308.     for(int i = 0; i < n; i++)
  309.         for(int j = 0; j < m; j++)
  310.             matrix[i][j] *= scalar;
  311. }
  312.  
  313. void matrix_minus_scalar(double** matrix, double scalar, int n, int m){
  314.     for(int i = 0; i < n; i++)
  315.         for(int j = 0; j < m; j++)
  316.             matrix[i][j] -= scalar;
  317. }
  318.  
  319. void complex_matrix_exponential(fftw_complex** target, double** source, int n, int m){
  320.     for(int i = 0; i < n; i++){
  321.         for(int j = 0; j < m; j++){
  322.             target[i][j][REAL] = complex_exponential_real_part(source[i][j]);
  323.             target[i][j][IMAG] = complex_exponential_imag_part(source[i][j]);
  324.         }
  325.     }
  326. }
  327.  
  328. void complex_matrix_member_per_member(fftw_complex** matrix1, fftw_complex** matrix2, fftw_complex** target, int n, int m){
  329.     double a, b, c, d;
  330.     for(int i = 0; i < n; i++){
  331.         for(int j = 0; j < m; j++){
  332.             a = matrix1[i][j][REAL];
  333.             b = matrix1[i][j][IMAG];
  334.             c = matrix2[i][j][REAL];
  335.             d = matrix2[i][j][IMAG];
  336.             target[i][j][REAL] = (a*c) - (b*d);
  337.             target[i][j][IMAG] = (b*c) + (a*d);
  338.         }
  339.     }
  340. }
  341.  
  342. void matrix_member_per_member(double** matrix1, double** matrix2, double** target, int n, int m){
  343.     for(int i = 0; i < n; i++)
  344.         for(int j = 0; j < m; j++)
  345.             target[i][j] = matrix1[i][j]*matrix2[i][j];
  346. }
  347.  
  348. void matrix_ratio_member_per_member(double** matrix1, double** matrix2, double** target, int n, int m){
  349.     for(int i = 0; i < n; i++)
  350.         for(int j = 0; j < m; j++)
  351.             target[i][j] = matrix1[i][j] / matrix2[i][j];
  352. }
  353.  
  354. void real_matrix_times_complex_matrix_scalar(double** rmatrix, fftw_complex** cmatrix, fftw_complex** target, int n, int m){
  355.     for(int i = 0; i < n; i++)
  356.         for(int j = 0; j < m; j++){
  357.             target[i][j][REAL] = cmatrix[i][j][REAL]*rmatrix[i][j];
  358.             target[i][j][IMAG] = cmatrix[i][j][IMAG]*rmatrix[i][j];
  359.         }
  360. }
  361.  
  362. void fast_time_windowing(double* window, fftw_complex** cmatrix, fftw_complex** target, int M, int N){
  363.     for(int m = 0; m < M; m++)
  364.         for(int n = 0; n < N; n++){
  365.             target[m][n][REAL] = cmatrix[m][n][REAL]*window[m];
  366.             target[m][n][IMAG] = cmatrix[m][n][IMAG]*window[m];
  367.         }
  368. }
  369.  
  370. void slow_time_windowing_cube(double* window, fftw_complex*** cmatrix, fftw_complex*** target, int L, int M, int N){
  371.     for(int l = 0; l < L; l++)
  372.         for(int m = 0; m < M; m++)
  373.             for(int n = 0; n < N; n++){
  374.                 target[l][m][n][REAL] = cmatrix[l][m][n][REAL]*window[m];
  375.                 target[l][m][n][IMAG] = cmatrix[l][m][n][IMAG]*window[m];
  376.             }
  377. }
  378.  
  379. void matrix_minus_matrix(double** matrix1, double** matrix2, double** target, int n, int m){
  380.     for(int i = 0; i < n; i++)
  381.         for(int j = 0; j < m; j++)
  382.             target[i][j] = matrix1[i][j] - matrix2[i][j];
  383. }
  384.  
  385. void matrix_scalar_sqrt(double** matrix, double** target, int n, int m){
  386.     for(int i = 0; i < n; i++)
  387.         for(int j = 0; j < m; j++)
  388.             target[i][j] = sqrt(matrix[i][j]);
  389. }
  390.  
  391. void matrix_sind(double** matrix, int n, int m){
  392.     for(int i = 0; i < n; i++)
  393.         for(int j = 0; j < m; j++)
  394.             matrix[i][j] = (sin(matrix[i][j]*PI/180));
  395. }
  396.  
  397. void complex_matrix_times_complex_matrix(fftw_complex** matrix1, fftw_complex** matrix2, fftw_complex** target, int n1, int m1, int n2, int m2){
  398.     double a, b, c, d;
  399.     for(int i = 0; i < n1; i++)
  400.         for(int j = 0; j < m2; j++)
  401.             for(int k = 0; k < m1; k++){
  402.                 a = matrix1[i][k][REAL];
  403.                 b = matrix1[i][k][IMAG];
  404.                 c = matrix2[k][j][REAL];
  405.                 d = matrix2[k][j][IMAG];
  406.                 target[i][j][REAL] += (a*c) - (b*d);
  407.                 target[i][j][IMAG] += (b*c) + (a*d);
  408.             }
  409. }
  410.  
  411. void complex_matrix_transposition(fftw_complex** source, fftw_complex** target, int n, int m){
  412.     for(int i = 0; i < m; i++)
  413.         for(int j = 0; j < n; j++){
  414.             target[i][j][REAL] = source[j][i][REAL];
  415.             target[i][j][IMAG] = source[j][i][IMAG];
  416.         }
  417. }
  418.  
  419. fftw_complex*** Frequency_Scaling_Algorithm(input_params ip, fftw_complex*** RawData, int Nmod, int Nsweeps, int Nch, int* Nrange_out){
  420.  
  421.     double              pi = 3.14159265;
  422.     int                 c = 3e8;
  423.     fftw_complex***     RangeData = zeros_3D_complex(ip.Nfft_range, Nsweeps, Nch);
  424.     double*             fast_time = linspace(0, ip.Tmod, Nmod);
  425.     double*             fr_ax = linspace(0, ip.Fs, Nmod);
  426.     double*             fr_ax_pad = linspace(0, ip.Fs, ip.Nfft_range);
  427.     double              r0 = (ip.range_ax[ip.rg_min_index]+ip.range_ax[ip.rg_max_index])/2;
  428.     fftw_complex**      TimeDopplerData_appo = create_complex_matrix(Nmod, Nsweeps);
  429.     fftw_complex        in[Nsweeps];
  430.     fftw_complex        out[Nsweeps];
  431.  
  432.     /////////////////////////////////////////////////////////////////////////
  433.     //print_1D_array_real("Fast time", fast_time, Nmod);                                    //OK
  434.     //print_1D_array_real("Range frequency axis", fr_ax, Nmod);                             //OK
  435.     //print_1D_array_real("Range frequency axis padded", fr_ax_pad, ip.Nfft_range);         //OK
  436.     //cout << "\n";
  437.     /////////////////////////////////////////////////////////////////////////
  438.  
  439.     for(int ch = 0; ch < Nch; ch++){
  440.  
  441.         /*Step 1: Slow Time FFT*/
  442.         for(int i = 0; i < Nmod; i++){
  443.             for(int j = 0; j < Nsweeps; j++){
  444.                 in[j][REAL] = RawData[i][j][ch][REAL];
  445.                 in[j][IMAG] = RawData[i][j][ch][IMAG];
  446.             }
  447.             fft(in, out, Nsweeps);
  448.  
  449.             /////////////////////////////////////////////////////////////////////////
  450.             //print_1D_array_complex("Slow Time FFT in row", in, Nsweeps);                  //OK
  451.             //print_1D_array_complex("Slow Time FFT out row", out, Nsweeps);                //OK
  452.             /////////////////////////////////////////////////////////////////////////
  453.  
  454.             fftshift(out, Nsweeps);
  455.  
  456.             ////////////////////////////////////////////////////////////////////////
  457.             //print_1D_array_complex("Slow Time FFT shift out row", out, Nsweeps);           //OK
  458.             /////////////////////////////////////////////////////////////////////////
  459.  
  460.             copy_vector(out, TimeDopplerData_appo[i], Nsweeps);
  461.         }
  462.  
  463.         /////////////////////////////////////////////////////////////////////////
  464. //      print_2D_array_complex("Step 1 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);      //OK
  465. //      cout << "\n";
  466.         /////////////////////////////////////////////////////////////////////////
  467.  
  468.         /*Step 2: frequency scaling*/
  469.         double          va_r[Nsweeps];
  470.         double          fd_a[Nsweeps];
  471.         double          beta[Nsweeps];
  472.         double**        matrix1 = create_real_matrix(Nmod, Nsweeps);
  473.         double**        matrix2 = create_real_matrix(Nmod, Nsweeps);
  474.         double          beta_n[Nsweeps];
  475.         fftw_complex**  H_FS = create_complex_matrix(Nmod, Nsweeps);
  476.         for(int i = 0; i < Nsweeps; i++)
  477.             va_r[i] = sin(atan((ip.xa[i]+ip.xa_offset[ch])/sqrt((r0*r0)+(ip.xa_offset[ch]*ip.xa_offset[ch]))))*ip.va;
  478.         scalar_times_vector(va_r, fd_a, -2*ip.f0/c, Nsweeps);
  479.  
  480.         /////////////////////////////////////////////////////////////////////////
  481.         //print_1D_array_real("va_r", va_r, Nsweeps);                                   //OK
  482.         //print_1D_array_real("fd_a", fd_a, Nsweeps);                                   //OK
  483.         //cout << "\n";
  484.         /////////////////////////////////////////////////////////////////////////
  485.  
  486.         scalar_times_vector(fd_a, beta, c/ip.f0/ip.va/2, Nsweeps);
  487.  
  488.         /////////////////////////////////////////////////////////////////////////
  489.         //print_1D_array_real("beta", beta, Nsweeps);                                   //OK
  490.         //cout << "\n";
  491.         /////////////////////////////////////////////////////////////////////////
  492.  
  493.         for(int i = 0; i < Nsweeps; i++)
  494.             beta[i] = sqrt(1 - (beta[i]*beta[i]));
  495.         for(int i = 0; i < Nsweeps; i++)
  496.             beta_n[i] = 1 - beta[i];
  497.         repmat_rowvector_matrix(beta_n, matrix1, Nmod, Nsweeps);
  498.         repmat_colvector_matrix(fast_time, matrix2, Nmod, Nsweeps);
  499.         matrix_scalar_square(matrix2, Nmod, Nsweeps);
  500.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  501.         matrix_times_scalar(matrix1, 3.1416*ip.chirp_rate, Nmod, Nsweeps);
  502.         complex_matrix_exponential(H_FS, matrix1, Nmod, Nsweeps);
  503.         complex_matrix_member_per_member(TimeDopplerData_appo, H_FS, TimeDopplerData_appo, Nmod, Nsweeps);
  504.  
  505.         /////////////////////////////////////////////////////////////////////////
  506. //      print_2D_array_complex("Step 2 H_FS", H_FS, Nmod, Nsweeps);                     //OK
  507. //      cout << "\n";
  508. //      print_2D_array_complex("Step 2 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);                     //OK
  509. //      cout << "\n";
  510.         /////////////////////////////////////////////////////////////////////////
  511.  
  512.         /*step 3: fast-time FFT*/
  513.         fftw_complex**  RangeDopplerData_appo = create_complex_matrix(Nmod, Nsweeps);
  514.         fftw_complex    in2[Nmod];
  515.         fftw_complex    out2[Nmod];
  516.         for(int j = 0; j < Nsweeps; j++){
  517.             for(int i = 0; i < Nmod; i++){                                   //Errore: Nsweeps al posto di Nmod
  518.                 in2[i][REAL] = TimeDopplerData_appo[i][j][REAL];            //Errore; H_FS al posto di TimeDopplerData_appo
  519.                 in2[i][IMAG] = TimeDopplerData_appo[i][j][IMAG];
  520.             }
  521.             fft(in2, out2, Nmod);
  522.             for(int i = 0; i < Nmod; i++){
  523.                 RangeDopplerData_appo[i][j][REAL] = out2[i][REAL];
  524.                 RangeDopplerData_appo[i][j][IMAG] = out2[i][IMAG];
  525.             }
  526.         }
  527.  
  528.         /////////////////////////////////////////////////////////////////////////
  529. //      print_2D_array_complex("Step 3 RangeDopplerData_appo", RangeDopplerData_appo, Nmod, Nsweeps);                       //OK
  530. //      cout << "\n";
  531.         /////////////////////////////////////////////////////////////////////////
  532.  
  533.         /*Step 4: RVP correction*/
  534.         fftw_complex**  H_RVP = create_complex_matrix(Nmod, Nsweeps);
  535.         double          fr_ax_n[Nmod];
  536.         for(int i = 0; i < Nmod; i++)
  537.             fr_ax_n[i] = fr_ax[i];                                                  //a cosa serve fr_ax_n?? mi apre uguale a fr_ax
  538.  
  539.         /////////////////////////////////////////////////////////////////////////
  540. //      print_1D_array_real("fr_ax", fr_ax, Nmod);                                  //OK
  541. //      print_1D_array_real("fr_ax_n", fr_ax_n, Nmod);                                  //OK
  542. //      cout << "\n";
  543.         /////////////////////////////////////////////////////////////////////////
  544.  
  545.         vector_scalar_square(fr_ax_n, Nmod);
  546.         repmat_colvector_matrix(fr_ax_n, matrix1, Nmod, Nsweeps);
  547.         repmat_rowvector_matrix(beta, matrix2, Nmod, Nsweeps);
  548.         matrix_times_scalar(matrix2, ip.chirp_rate, Nmod, Nsweeps);
  549.         matrix_ratio_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  550.         matrix_times_scalar(matrix1, -pi, Nmod, Nsweeps);
  551.         complex_matrix_exponential(H_RVP, matrix1, Nmod, Nsweeps);
  552.         complex_matrix_member_per_member(RangeDopplerData_appo, H_RVP, RangeDopplerData_appo, Nmod, Nsweeps);      
  553.  
  554.         /////////////////////////////////////////////////////////////////////////
  555. //      print_2D_array_complex("Step 4 RangeDopplerData_appo", RangeDopplerData_appo, Nmod, Nsweeps);                       //OK
  556. //      cout << "\n";
  557.         /////////////////////////////////////////////////////////////////////////
  558.  
  559.         /*Step 5: fast-time ifft*/
  560.         for(int j = 0; j < Nsweeps; j++){
  561.             for(int i = 0; i < Nmod; i++){                                   //Errore: Nsweeps al posto di Nmod
  562.                 in2[i][REAL] = RangeDopplerData_appo[i][j][REAL];
  563.                 in2[i][IMAG] = RangeDopplerData_appo[i][j][IMAG];
  564.             }
  565.             ifft(in2, out2, Nmod);
  566.             for(int i = 0; i < Nmod; i++){
  567.                 TimeDopplerData_appo[i][j][REAL] = out2[i][REAL];
  568.                 TimeDopplerData_appo[i][j][IMAG] = out2[i][IMAG];
  569.             }
  570.         }
  571.  
  572.         /////////////////////////////////////////////////////////////////////////
  573. //      print_2D_array_complex("Step 5 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);                     //OK
  574. //      cout << "\n";
  575.         /////////////////////////////////////////////////////////////////////////
  576.  
  577.         /*Step 6: inverse frequency scaling*/
  578.         fftw_complex** H_IFS = create_complex_matrix(Nmod, Nsweeps);
  579.         for(int i = 0; i < Nsweeps; i++)
  580.             beta_n[i] = beta[i]*(beta[i]-1);
  581.         repmat_rowvector_matrix(beta_n, matrix2, Nmod, Nsweeps);
  582.         repmat_colvector_matrix(fast_time, matrix1, Nmod, Nsweeps);
  583.         matrix_scalar_square(matrix1, Nmod, Nsweeps);
  584.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  585.         matrix_times_scalar(matrix1, ip.chirp_rate*pi, Nmod, Nsweeps);
  586.         complex_matrix_exponential(H_IFS, matrix1, Nmod, Nsweeps);
  587.         complex_matrix_member_per_member(TimeDopplerData_appo, H_IFS, TimeDopplerData_appo, Nmod, Nsweeps);
  588.  
  589.         /////////////////////////////////////////////////////////////////////////
  590. //      print_2D_array_complex("Step 6 H_IFS", H_IFS, Nmod, Nsweeps);                       //OK
  591. //      cout << "\n";
  592. //      print_2D_array_complex("Step 6 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);                     //OK
  593. //      cout << "\n";
  594.         /////////////////////////////////////////////////////////////////////////
  595.  
  596.         /*Step 7: Doppler frequency correction*/
  597.         fftw_complex** H_DFC = create_complex_matrix(Nmod, Nsweeps);
  598.         for(int i = 0; i < Nsweeps; i++)
  599.             beta_n[i] = fd_a[i]*beta[i];
  600.         repmat_rowvector_matrix(beta_n, matrix2, Nmod, Nsweeps);
  601.         repmat_colvector_matrix(fast_time, matrix1, Nmod, Nsweeps);
  602.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  603.         matrix_times_scalar(matrix1, -2*pi, Nmod, Nsweeps);
  604.         complex_matrix_exponential(H_DFC, matrix1, Nmod, Nsweeps);
  605.         complex_matrix_member_per_member(TimeDopplerData_appo, H_DFC, TimeDopplerData_appo, Nmod, Nsweeps);
  606.  
  607.         /////////////////////////////////////////////////////////////////////////
  608. //      print_2D_array_complex("Step 7 H_DFC", H_DFC, Nmod, Nsweeps);                       //OK
  609. //      cout << "\n";
  610. //      print_2D_array_complex("Step 7 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);                     //OK
  611. //      cout << "\n";
  612.         /////////////////////////////////////////////////////////////////////////
  613.  
  614.         /*Step 8: secondary range compression*/
  615.         double rc = r0;
  616.         fftw_complex** H_SRC = create_complex_matrix(Nmod, Nsweeps);
  617.         fftw_complex** H_SRC1 = create_complex_matrix(Nmod, Nsweeps);
  618.         fftw_complex** H_SRC2 = create_complex_matrix(Nmod, Nsweeps);      
  619.         /* part 1*/
  620.         for(int i = 0; i < Nsweeps; i++)
  621.             beta_n[i] = (beta[i]*beta[i]-1)/(beta[i]*beta[i]*beta[i]);
  622.         repmat_rowvector_matrix(beta, matrix1, Nmod, Nsweeps);
  623.         repmat_colvector_matrix(fast_time, matrix2, Nmod, Nsweeps);
  624.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  625.         matrix_minus_scalar(matrix1, 2*rc/c, Nmod, Nsweeps);
  626.         matrix_scalar_square(matrix1, Nmod, Nsweeps);
  627.         repmat_rowvector_matrix(beta_n, matrix2, Nmod, Nsweeps);
  628.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  629.         matrix_times_scalar(matrix1, 2*pi*rc*ip.chirp_rate*ip.chirp_rate/c/ip.f0, Nmod, Nsweeps);
  630.         complex_matrix_exponential(H_SRC1, matrix1, Nmod, Nsweeps);
  631.         /*part 2*/
  632.         for(int i = 0; i < Nsweeps; i++)
  633.             beta_n[i] /= (beta[i]*beta[i]);
  634.         repmat_rowvector_matrix(beta, matrix1, Nmod, Nsweeps);
  635.         repmat_colvector_matrix(fast_time, matrix2, Nmod, Nsweeps);
  636.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  637.         matrix_minus_scalar(matrix1, 2*rc/c, Nmod, Nsweeps);
  638.         matrix_scalar_cube(matrix1, Nmod, Nsweeps);
  639.         repmat_rowvector_matrix(beta_n, matrix2, Nmod, Nsweeps);
  640.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  641.         matrix_times_scalar(matrix1, -(2*pi*rc*ip.chirp_rate*ip.chirp_rate*ip.chirp_rate/ip.f0/ip.f0/ip.f0), Nmod, Nsweeps);
  642.         complex_matrix_exponential(H_SRC2, matrix1, Nmod, Nsweeps);
  643.         complex_matrix_member_per_member(H_SRC1, H_SRC2, H_SRC, Nmod, Nsweeps);
  644.         complex_matrix_member_per_member(TimeDopplerData_appo, H_SRC, TimeDopplerData_appo, Nmod, Nsweeps);
  645.          
  646.         /////////////////////////////////////////////////////////////////////////
  647. //      print_2D_array_complex("Step 8 H_DFC", H_SRC, Nmod, Nsweeps);                       //OK
  648. //      cout << "\n";
  649. //      print_2D_array_complex("Step 8 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);                     //OK
  650. //      cout << "\n";
  651.         /////////////////////////////////////////////////////////////////////////
  652.  
  653.  
  654.         /*Step 9: bulk range shift*/
  655.         fftw_complex** H_BS = create_complex_matrix(Nmod, Nsweeps);
  656.         for(int i = 0; i < Nsweeps; i++)
  657.             beta_n[i] = (1/beta[i]) -1;
  658.         repmat_rowvector_matrix(beta, matrix1, Nmod, Nsweeps);
  659.         repmat_colvector_matrix(fast_time, matrix2, Nmod, Nsweeps);
  660.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  661.         matrix_minus_scalar(matrix1, 2*rc/c, Nmod, Nsweeps);
  662.         repmat_rowvector_matrix(beta_n, matrix2, Nmod, Nsweeps);
  663.         matrix_member_per_member(matrix1, matrix2, matrix1, Nmod, Nsweeps);
  664.         matrix_times_scalar(matrix1, 4*pi*ip.chirp_rate/c*rc, Nmod, Nsweeps);
  665.         complex_matrix_exponential(H_BS, matrix1, Nmod, Nsweeps);
  666.         complex_matrix_member_per_member(TimeDopplerData_appo, H_BS, TimeDopplerData_appo, Nmod, Nsweeps);
  667.  
  668.         /////////////////////////////////////////////////////////////////////////
  669. //      print_2D_array_complex("Step 9 H_BS", H_BS, Nmod, Nsweeps);                     //OK
  670. //      cout << "\n";
  671. //      print_2D_array_complex("Step 9 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);                     //OK
  672. //      cout << "\n";
  673.         /////////////////////////////////////////////////////////////////////////
  674.  
  675.         /*Step 10: range windowing*/
  676.         fftw_complex**  RawData_filt = create_complex_matrix(Nmod, Nsweeps);
  677.         fast_time_windowing(ip.Win_range, TimeDopplerData_appo, RawData_filt, Nmod, Nsweeps);
  678.          
  679.         /////////////////////////////////////////////////////////////////////////
  680. //      print_2D_array_complex("Step 10 TimeDopplerData_appo", TimeDopplerData_appo, Nmod, Nsweeps);                        //OK
  681. //      cout << "\n";
  682. //      print_2D_array_complex("Step 10 RawData_filt", RawData_filt, Nmod, Nsweeps);                        //OK
  683. //      cout << "\n";
  684.         /////////////////////////////////////////////////////////////////////////
  685.  
  686.         /*Step 11: final range compressione*/
  687.         fftw_complex**  RangeDopplerData_filt = create_complex_matrix(ip.Nfft_range, Nsweeps);
  688.         fftw_complex    in3[ip.Nfft_range];
  689.         fftw_complex    out3[ip.Nfft_range];
  690.         for(int j = 0; j < Nsweeps; j++){
  691.             for(int i = 0; i < ip.Nfft_range; i++){
  692.                 if(i < Nmod){
  693.                     in3[i][REAL] = RawData_filt[i][j][REAL];            //Errore: usato in2 e out2 al posto di in3 e out3
  694.                     in3[i][IMAG] = RawData_filt[i][j][IMAG];
  695.                 }else{
  696.                     in3[i][REAL] = 0;
  697.                     in3[i][IMAG] = 0;
  698.                 }                  
  699.             }
  700.             fft(in3, out3, ip.Nfft_range);
  701.             for(int i = 0; i < ip.Nfft_range; i++){
  702.                 RangeDopplerData_filt[i][j][REAL] = out3[i][REAL];
  703.                 RangeDopplerData_filt[i][j][IMAG] = out3[i][IMAG];
  704.             }
  705.         }
  706.  
  707.         /////////////////////////////////////////////////////////////////////////
  708. //      print_2D_array_complex("Step 11 RangeDopplerData_filt", RangeDopplerData_filt, ip.Nfft_range, Nsweeps);                     //OK
  709. //      cout << "\n";
  710.         /////////////////////////////////////////////////////////////////////////
  711.  
  712.         /*Step 12: phase preservation*/
  713.         fftw_complex**  H_PPC = create_complex_matrix(ip.Nfft_range, Nsweeps);
  714.         double**        matrix3 = create_real_matrix(ip.Nfft_range, Nsweeps);
  715.         double**        matrix4 = create_real_matrix(ip.Nfft_range, Nsweeps);
  716.         repmat_colvector_matrix(fr_ax_pad, matrix3, ip.Nfft_range, Nsweeps);
  717.         repmat_rowvector_matrix(beta, matrix4, ip.Nfft_range, Nsweeps);
  718.         matrix_ratio_member_per_member(matrix3, matrix4, matrix3, ip.Nfft_range, Nsweeps);
  719.         matrix_times_scalar(matrix3, (4*pi*rc/c), ip.Nfft_range, Nsweeps);
  720.         complex_matrix_exponential(H_PPC, matrix3, ip.Nfft_range, Nsweeps);
  721.         complex_matrix_member_per_member(RangeDopplerData_filt, H_PPC, RangeDopplerData_filt, ip.Nfft_range, Nsweeps);
  722.  
  723.         /////////////////////////////////////////////////////////////////////////
  724. //      print_2D_array_complex("Step 12 H_PPC", H_PPC, ip.Nfft_range, Nsweeps);                     //OK
  725. //      cout << "\n";
  726. //      print_2D_array_complex("Step 12 RangeDopplerData_filt", RangeDopplerData_filt, ip.Nfft_range, Nsweeps);                     //OK
  727. //      cout << "\n";
  728.         /////////////////////////////////////////////////////////////////////////
  729.  
  730.         /*Step 13: final azimuth IFFT*/
  731.         fftw_complex**  RangeData_filt = create_complex_matrix(ip.Nfft_range, Nsweeps);
  732.         for(int i = 0; i < ip.Nfft_range; i++)
  733.             ifftshift(RangeDopplerData_filt[i], Nsweeps);
  734.         for(int i = 0; i < ip.Nfft_range; i++){
  735.             for(int j = 0; j < Nsweeps; j++){
  736.                 in[j][REAL] = RangeDopplerData_filt[i][j][REAL];
  737.                 in[j][IMAG] = RangeDopplerData_filt[i][j][IMAG];
  738.             }              
  739.             ifft(in, out, Nsweeps);                                     //Errore: usato fft al posto ifft
  740.             //fftshift(out, Nsweeps);
  741.             copy_vector(out, RangeData_filt[i], Nsweeps);  
  742.         }
  743.  
  744.         /////////////////////////////////////////////////////////////////////////
  745. //      print_2D_array_complex("Step 13 RangeData_filt", RangeData_filt, ip.Nfft_range, Nsweeps);                       //OK
  746. //      cout << "\n";
  747.         /////////////////////////////////////////////////////////////////////////
  748.  
  749.  
  750.         /*Step 14: azimuth compression*/
  751.         double*         xa_n = create_real_vector(Nsweeps);
  752.         double**        RR_stolt = create_real_matrix(ip.Nfft_range, Nsweeps);
  753.         fftw_complex**  H_AMF = create_complex_matrix(ip.Nfft_range, Nsweeps);
  754.         for(int i = 0; i < Nsweeps; i++)
  755.             xa_n[i] = ip.xa[i] + ip.xa_offset[ch];
  756.  
  757.         /////////////////////////////////////////////////////////////////////////
  758. //      print_1D_array_real("Step 14 xa_n", xa_n, Nsweeps);         //OK
  759. //      cout << "\n";
  760.         /////////////////////////////////////////////////////////////////////////
  761.  
  762.         repmat_colvector_matrix(ip.range_ax, matrix3, ip.Nfft_range, Nsweeps);
  763.         repmat_rowvector_matrix(xa_n, matrix4, ip.Nfft_range, Nsweeps);
  764.  
  765.         /////////////////////////////////////////////////////////////////////////
  766. //      print_2D_array_real("Step 14 matrix3", matrix3, ip.Nfft_range, Nsweeps);            //OK
  767. //      cout << "\n";
  768. //      print_2D_array_real("Step 14 matrix4", matrix4, ip.Nfft_range, Nsweeps);            //OK
  769. //      cout << "\n";
  770.         /////////////////////////////////////////////////////////////////////////
  771.  
  772.         matrix_scalar_square(matrix3, ip.Nfft_range, Nsweeps);
  773.         matrix_scalar_square(matrix4, ip.Nfft_range, Nsweeps);
  774.  
  775.         /////////////////////////////////////////////////////////////////////////
  776. //      print_2D_array_real("Step 14 matrix3", matrix3, ip.Nfft_range, Nsweeps);            //OK
  777. //      cout << "\n";
  778. //      print_2D_array_real("Step 14 matrix4", matrix4, ip.Nfft_range, Nsweeps);            //OK
  779. //      cout << "\n";
  780.         /////////////////////////////////////////////////////////////////////////
  781.  
  782.         matrix_minus_matrix(matrix3, matrix4, matrix4, ip.Nfft_range, Nsweeps);
  783.         matrix_scalar_sqrt(matrix4, RR_stolt, ip.Nfft_range, Nsweeps);
  784.  
  785.         /////////////////////////////////////////////////////////////////////////
  786. //      print_2D_array_real("Step 14 matrix4", matrix4, ip.Nfft_range, Nsweeps);            //OK
  787. //      cout << "\n";
  788. //      print_2D_array_real("Step 14 RR_stolt", RR_stolt, ip.Nfft_range, Nsweeps);          //OK
  789. //      cout << "\n";
  790.         /////////////////////////////////////////////////////////////////////////
  791.  
  792.         matrix_times_scalar(RR_stolt, 4*pi*ip.f0*rc/c, ip.Nfft_range, Nsweeps);
  793.         complex_matrix_exponential(H_AMF, RR_stolt, ip.Nfft_range, Nsweeps);
  794.         complex_matrix_member_per_member(RangeData_filt, H_AMF, RangeData_filt, ip.Nfft_range, Nsweeps);
  795.          
  796.         /////////////////////////////////////////////////////////////////////////
  797. //      print_2D_array_real("Step 14 RR_stolt", RR_stolt, ip.Nfft_range, Nsweeps);                      //OK
  798. //      cout << "\n";
  799. //      print_2D_array_complex("Step 14 H_AMF", H_AMF, ip.Nfft_range, Nsweeps);                     //OK
  800. //      cout << "\n";
  801. //      print_2D_array_complex("Step 14 RangeData_filt", RangeData_filt, ip.Nfft_range, Nsweeps);                       //OK
  802. //      cout << "\n";
  803.         /////////////////////////////////////////////////////////////////////////
  804.  
  805.         /*Finish*/
  806.         for(int i = 0; i < ip.Nfft_range; i++)
  807.             for(int j = 0; j < Nsweeps; j++){
  808.                 RangeData[i][j][ch][REAL] = RangeData_filt[i][j][REAL];
  809.                 RangeData[i][j][ch][IMAG] = RangeData_filt[i][j][IMAG];
  810.             }
  811.  
  812.     }
  813.  
  814.     /*Monostatic equivalent correction*/
  815.     double**        matrix5 = create_real_matrix(ip.Nfft_range, Nsweeps);
  816.     double**        eq_mono_corr = create_real_matrix(ip.Nfft_range, Nsweeps);
  817.     fftw_complex**  H_eq_mono_corr = create_complex_matrix(ip.Nfft_range, Nsweeps);
  818.     fftw_complex**  matrix_supp = create_complex_matrix(ip.Nfft_range, Nsweeps);
  819.  
  820.     for(int ch = 0; ch < Nch; ch++){
  821.         repmat_colvector_matrix(ip.range_ax, matrix5, ip.Nfft_range, Nsweeps);
  822.         matrix_scalar_square(matrix5, ip.Nfft_range, Nsweeps);
  823.         matrix_minus_scalar(matrix5, (ip.xa_offset[ch]*ip.xa_offset[ch]), ip.Nfft_range, Nsweeps);
  824.         matrix_scalar_sqrt(matrix5, eq_mono_corr, ip.Nfft_range, Nsweeps);
  825.         matrix_times_scalar(eq_mono_corr, 4*pi*ip.f0/c, ip.Nfft_range, Nsweeps);
  826.         complex_matrix_exponential(H_eq_mono_corr, eq_mono_corr, ip.Nfft_range, Nsweeps);
  827.         for(int i = 0; i < ip.Nfft_range; i++)
  828.             for(int j = 0; j < Nsweeps; j++){
  829.                 matrix_supp[i][j][REAL] = RangeData[i][j][ch][REAL];
  830.                 matrix_supp[i][j][IMAG] = RangeData[i][j][ch][IMAG];
  831.             }
  832.         complex_matrix_member_per_member(matrix_supp, H_eq_mono_corr, matrix_supp, ip.Nfft_range, Nsweeps);
  833.         for(int i = 0; i < ip.Nfft_range; i++)
  834.             for(int j = 0; j < Nsweeps; j++){
  835.                 RangeData[i][j][ch][REAL] = matrix_supp[i][j][REAL];
  836.                 RangeData[i][j][ch][IMAG] = matrix_supp[i][j][IMAG];
  837.             }  
  838.     }
  839.  
  840.     /////////////////////////////////////////////////////////////////////////
  841. //  print_3D_array_complex("RangeData after monostatic range correction", RangeData, ip.Nfft_range, Nsweeps, Nch);          //OK
  842. //  cout << "\n";
  843.     /////////////////////////////////////////////////////////////////////////
  844.  
  845.  
  846.     /*Range gating*/  //credo fosse tutto sbagliato
  847.     int                 Nrange = ip.rg_max_index - ip.rg_min_index + 1;
  848.     fftw_complex***     RangeData_n = zeros_3D_complex(Nrange, Nsweeps, Nch);
  849.     int                 i_n = 0;
  850.     for(int i = ip.rg_min_index; i <= ip.rg_max_index; i++){
  851.         for(int j = 0; j < Nsweeps; j++)
  852.             for(int ch = 0; ch < Nch; ch++){
  853.                     RangeData_n[i_n][j][ch][REAL] = RangeData[i][j][ch][REAL];
  854.                     RangeData_n[i_n][j][ch][IMAG] = RangeData[i][j][ch][IMAG];
  855. //                  cout << "i=" << i << "i_n=" << i_n  << "j=" << j << "ch=" << ch << "\n";
  856.             }
  857.         i_n++;
  858.     }
  859.     *Nrange_out = Nrange;
  860.  
  861.     /////////////////////////////////////////////////////////////////////////
  862. //  print_3D_array_complex("RangeData before range gating", RangeData, ip.Nfft_range, Nsweeps, Nch);            //OK
  863. //  cout << "\n";
  864. //  print_3D_array_complex("RangeData after range gating", RangeData_n, Nrange, Nsweeps, Nch);          //OK
  865. //  cout << "\n";
  866.     /////////////////////////////////////////////////////////////////////////
  867.  
  868.  
  869.  
  870.     /*Slow-time windowing*/
  871. //  double*** matrix6;
  872. //  matrix6 = real_zeros(Nrange, Nsweeps, Nch);
  873. //  repmat_rowvector_tensor(ip.Win_azimuth, matrix6, Nrange, Nsweeps, Nch); //operazione evitabile
  874. //  for(int i = 0; i < Nrange; i++)
  875. //      for(int j = 0; j < Nsweeps; j++)
  876. //          for(int k = 0; k < Nch; k++){
  877. //              RangeData_n[i][j][k][REAL] *= matrix6[i][j][k];
  878. //              RangeData_n[i][j][k][IMAG] *= matrix6[i][j][k];
  879. //          }
  880.  
  881.     slow_time_windowing_cube(ip.Win_azimuth, RangeData_n, RangeData_n, Nrange, Nsweeps, Nch);
  882.  
  883.     /////////////////////////////////////////////////////////////////////////
  884. //  print_3D_array_complex("RangeData after slow-time windowing", RangeData_n, Nrange, Nsweeps, Nch);           //OK
  885. //  cout << "\n";
  886.     /////////////////////////////////////////////////////////////////////////
  887.  
  888.     /*Beamforming (DFT)*/
  889.     fftw_complex***     IMM = zeros_3D_complex(Nrange, ip.azimuth_ax_size, Nch);
  890.     fftw_complex**      W = create_complex_matrix(ip.azimuth_ax_size, Nsweeps);
  891.     fftw_complex**      W_transp = create_complex_matrix(Nsweeps, ip.azimuth_ax_size);
  892.     double**            X_a = create_real_matrix(ip.azimuth_ax_size, Nsweeps);
  893.     double**            AZ_ax = create_real_matrix(ip.azimuth_ax_size, Nsweeps);
  894.     repmat_rowvector_matrix(ip.xa, X_a, ip.azimuth_ax_size, Nsweeps);
  895.     repmat_colvector_matrix(ip.azimuth_ax, AZ_ax, ip.azimuth_ax_size, Nsweeps);
  896.     matrix_sind(AZ_ax, ip.azimuth_ax_size, Nsweeps);
  897.     matrix_member_per_member(X_a, AZ_ax, X_a, ip.azimuth_ax_size, Nsweeps);
  898.     matrix_times_scalar(X_a, 4*pi*ip.f0/c, ip.azimuth_ax_size, Nsweeps);
  899.     complex_matrix_exponential(W, X_a, ip.azimuth_ax_size, Nsweeps);
  900.     complex_matrix_transposition(W, W_transp, ip.azimuth_ax_size, Nsweeps);
  901.     fftw_complex** matrix7 = create_complex_matrix(Nrange, Nsweeps);
  902.     fftw_complex** matrix8 = create_complex_matrix(Nrange, ip.azimuth_ax_size);
  903.     for(int ch = 0; ch < Nch; ch++){
  904.         for(int i = 0; i < Nrange; i++)
  905.             for(int j = 0; j < Nsweeps; j++){
  906.                 matrix7[i][j][REAL] = RangeData_n[i][j][ch][REAL];
  907.                 matrix7[i][j][IMAG] = RangeData_n[i][j][ch][IMAG];
  908.             }
  909.         complex_matrix_times_complex_matrix(matrix7, W_transp, matrix8, Nrange, Nsweeps, Nsweeps, ip.azimuth_ax_size);
  910.         for(int i = 0; i < Nrange; i++)
  911.             for(int j = 0; j < ip.azimuth_ax_size; j++){
  912.                 IMM[i][j][ch][REAL] = matrix8[i][j][REAL];
  913.                 IMM[i][j][ch][IMAG] = matrix8[i][j][IMAG];
  914.             }
  915.     }
  916.  
  917.     /////////////////////////////////////////////////////////////////////////
  918. //  print_3D_array_complex("IMM", IMM, Nrange, Nsweeps, Nch);           //OK
  919. //  cout << "\n";
  920.     /////////////////////////////////////////////////////////////////////////
  921.  
  922.     return IMM;
  923.  
  924. }
  925.  
  926. int read_from_file_2D_array_complex(fftw_complex** data, char* filename, int L, int M){
  927.  
  928.     FILE*               pFile  = NULL;
  929.     unsigned long int   lSize  = 0;
  930.     size_t              result = 0;
  931.     long int            i      = 0;
  932.  
  933.     //Open file
  934.     pFile = fopen(filename, "rb");
  935.     if (pFile==NULL){
  936.         fputs ("File open error\n",stderr);
  937.         return (-1);
  938.     }
  939.  
  940.     //Get file size
  941.     fseek (pFile , 0 , SEEK_END);
  942.     lSize = ftell(pFile);
  943.     rewind (pFile);
  944.  
  945.     //Copy the file into the buffer
  946.     double* buffer_real = new double[lSize/2];
  947.     double* buffer_imag = new double[lSize/2];
  948.     result = fread(buffer_real,1,lSize/2,pFile);
  949.     if (result!=lSize/2){
  950.         fputs ("File read error\n",stderr);
  951.         return (-1);
  952.     }
  953.     result = fread(buffer_imag,1,lSize/2,pFile);
  954.     if (result!=lSize/2){
  955.         fputs ("File read error\n",stderr);
  956.         return (-1);
  957.     }
  958.     for(int l = 0; l < L; l++)
  959.         for(int m = 0; m < M; m++){
  960.             data[l][m][REAL] = buffer_real[i];
  961.             data[l][m][IMAG] = buffer_imag[i];
  962.             i++;
  963.         }
  964.  
  965.     //Close file
  966.     fclose(pFile);
  967.  
  968.     return 0;
  969. }
  970.  
  971.  
  972. int read_from_file_3D_array_complex(fftw_complex*** data, char* filename, int L, int M, int N){
  973.  
  974.     FILE*               pFile  = NULL;
  975.     unsigned long int   lSize  = 0;
  976.     size_t              result = 0;
  977.     long int            i      = 0;
  978.  
  979.     //Open file
  980.     pFile = fopen(filename, "rb");
  981.     if (pFile==NULL){
  982.         fputs ("File open error\n",stderr);
  983.         return (-1);
  984.     }
  985.  
  986.     //Get file size
  987.     fseek (pFile , 0 , SEEK_END);
  988.     lSize = ftell(pFile);
  989.     rewind (pFile);
  990.  
  991.     //Copy the file into the buffer
  992.     double* buffer_real = new double[lSize/2];
  993.     double* buffer_imag = new double[lSize/2];
  994.     result = fread(buffer_real,1,lSize/2,pFile);
  995.     if (result!=lSize/2){
  996.         fputs ("File read error\n",stderr);
  997.         return (-1);
  998.     }
  999.     result = fread(buffer_imag,1,lSize/2,pFile);
  1000.     if (result!=lSize/2){
  1001.         fputs ("File read error\n",stderr);
  1002.         return (-1);
  1003.     }
  1004.     for(int l = 0; l < L; l++)
  1005.         for(int m = 0; m < M; m++)
  1006.             for(int n = 0; n < N; n++){
  1007.                 data[l][m][n][REAL] = buffer_real[i];
  1008.                 data[l][m][n][IMAG] = buffer_imag[i];
  1009.                 i++;
  1010.             }
  1011.  
  1012.     //Close file
  1013.     fclose(pFile);
  1014.  
  1015.     return 0;
  1016. }
  1017.  
  1018.  
  1019.  
  1020. int write_to_file_3D_array_complex(fftw_complex*** data, char* filename, int L, int M, int N){
  1021.  
  1022.     FILE*       pFile  = NULL;
  1023.  
  1024.     //Open file
  1025.     pFile = fopen(filename, "wb");
  1026.     if (pFile==NULL){
  1027.         fputs ("File open error\n",stderr);
  1028.         return (-1);
  1029.     }
  1030.  
  1031.     //Write buffer to file
  1032.     for(int l = 0; l < L; l++)
  1033.         for(int m = 0; m < M; m++)
  1034.             for(int n = 0; n < N; n++){
  1035.                 fwrite(&data[l][m][n][REAL] , sizeof(double), 1, pFile);
  1036.                 fwrite(&data[l][m][n][IMAG] , sizeof(double), 1, pFile);
  1037.             }
  1038.  
  1039.     //Close file
  1040.     fclose(pFile);
  1041.  
  1042.     return 0;
  1043. }
  1044.  
  1045.  
  1046. int write_to_file_2D_array_complex(fftw_complex** data, char* filename, int L, int M){
  1047.  
  1048.     FILE*       pFile  = NULL;
  1049.  
  1050.     //Open file
  1051.     pFile = fopen(filename, "wb");
  1052.     if (pFile==NULL){
  1053.         fputs ("File open error\n",stderr);
  1054.         return (-1);
  1055.     }
  1056.  
  1057.     //Write buffer to file
  1058.     for(int l = 0; l < L; l++)
  1059.         for(int m = 0; m < M; m++){
  1060.             fwrite(&data[l][m][REAL] , sizeof(double), 1, pFile);
  1061.             fwrite(&data[l][m][IMAG] , sizeof(double), 1, pFile);
  1062.         }
  1063.  
  1064.     //Close file
  1065.     fclose(pFile);
  1066.  
  1067.     return 0;
  1068. }
  1069.  
  1070.  
  1071. int main(){
  1072.  
  1073.     int                 Nmod    = 1222;
  1074.     int                 Nsweeps = 235;//586;
  1075.     int                 Nch     = 4;
  1076.     int                 Nrange  = 0;
  1077.     fftw_complex***     RawData = NULL;
  1078.     fftw_complex***     Complex_Image = NULL;
  1079.     input_params        InputParams;
  1080.  
  1081.     InputParams.f0              = 13.6e9;
  1082.     InputParams.chirp_rate      = 2.0040e12;
  1083.     InputParams.Fs              = 2.5e6;
  1084.     InputParams.Tmod            = 4.8900e-4;
  1085.     InputParams.va              = -0.1164;
  1086.     InputParams.Nfft_range      = 2048;
  1087.     InputParams.azimuth_ax_size = 250;
  1088.     InputParams.xa              = linspace(-0.6822, 0.6822, Nsweeps);
  1089.     InputParams.xa_offset       = create_real_vector(Nch);
  1090.     InputParams.range_ax        = linspace(1.9924, 40.0055, InputParams.Nfft_range );
  1091.     InputParams.azimuth_ax      = linspace(-25,25,InputParams.azimuth_ax_size);
  1092.     InputParams.Win_range       = create_real_vector(Nmod);
  1093.     InputParams.Win_azimuth     = create_real_vector(InputParams.azimuth_ax_size);
  1094.     InputParams.rg_min_index    = 84;
  1095.     InputParams.rg_max_index    = 1325;
  1096.  
  1097.  
  1098.     //Define range window
  1099.     for(int i = 0; i < Nmod; i++)
  1100.         InputParams.Win_range[i] = 1;
  1101.  
  1102.     //Define azimuth window
  1103.     for(int j = 0; j < Nsweeps; j++)
  1104.         InputParams.Win_azimuth[j] = 1;
  1105.  
  1106.     //Initialize RaData
  1107.     RawData = zeros_3D_complex(Nmod, Nsweeps, Nch);
  1108. //  for(int ch = 0; ch < Nch; ch++)
  1109. //      for(int i = 0; i < Nmod; i++)
  1110. //          for(int j = 0; j < Nsweeps; j++){
  1111. //              RawData[i][j][ch][REAL] = i;
  1112. //              RawData[i][j][ch][IMAG] = j;
  1113. //          }
  1114.  
  1115.     //Read RawData from binary file
  1116.     read_from_file_3D_array_complex(RawData, "/home/slischi/Desktop/RawData_test.bin", Nmod, Nsweeps, Nch);
  1117. //  print_3D_array_complex("RawData", RawData, Nmod, Nsweeps, Nch);
  1118.  
  1119. //  //Write RawData to file
  1120. //  write_to_file_3D_array_complex(RawData, "/home/slischi/Desktop/RawData_readback_test.bin", Nmod, Nsweeps, Nch);
  1121.  
  1122.  
  1123.  
  1124.     //Frequency Scaling Algorithm
  1125.     clock_t start = clock();
  1126.     Complex_Image = Frequency_Scaling_Algorithm(InputParams, RawData, Nmod, Nsweeps, Nch, &Nrange);
  1127.     clock_t end = clock();
  1128.     printf("Tempo di esecuzione =  %f secondi \n", ((double)(end - start)) / CLOCKS_PER_SEC);
  1129.  
  1130.  
  1131.     //Write Complex_Image to file
  1132.     write_to_file_3D_array_complex(Complex_Image, "/home/slischi/Desktop/Image_test.bin", Nrange, InputParams.azimuth_ax_size, Nch);
  1133.  
  1134.     cout << "Done\n";
  1135.  
  1136.     return 0;
  1137. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!