Added dct2.c with 1-D and 2-D DCT and inverse DCT functions, based off of the existin...
authorphilayres <philayres@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 25 Jul 2017 15:48:13 +0000 (15:48 +0000)
committerphilayres <philayres@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 25 Jul 2017 15:48:13 +0000 (15:48 +0000)
Provided simple test cases for DCT functions.
Took out the dct2 and idct2 stubs from c2wideband.c
Split out some reusable test helper functions for reuse (in t_helper.c/h).
Formatted c2wideband.c with correct tab spacing

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

codec2-dev/src/CMakeLists.txt
codec2-dev/src/c2wideband.c
codec2-dev/src/dct2.c [new file with mode: 0644]
codec2-dev/src/dct2.h [new file with mode: 0644]
codec2-dev/unittest/CMakeLists.txt
codec2-dev/unittest/t_helpers.c [new file with mode: 0644]
codec2-dev/unittest/t_helpers.h [new file with mode: 0644]
codec2-dev/unittest/tc2wideband.c
codec2-dev/unittest/tdct2.c [new file with mode: 0644]

index 7d758d37233eb8579a4dbeb832b1437e80c6642d..66a6fe1595ebab5ba0255268e5292cef84f3b0f7 100644 (file)
@@ -230,6 +230,7 @@ set(CODEC2_SRCS
     varicode.c
     modem_stats.c
     c2wideband.c
+    dct2.c
 )
 
 set(CODEC2_PUBLIC_HEADERS
index c41a71581e1e58841bb988dd00e185cd7dd1afd2..684b7e488b4eb4408cba94bcfea04fad882fb1fa 100644 (file)
@@ -49,6 +49,7 @@
 #include "machdep.h"
 #include "codec2.h"
 #include "newamp1.h"
+#include "dct2.h"
 
 #include "c2wideband.h"
 #include "c2wideband_map.h"
@@ -67,519 +68,503 @@ void setup_map(WIDEBAND_MAP * wb_map, int Nt, int K);
 //TODO: move all this to the standard codec2.c functions
 #ifdef IGNORE
 
-int main(int argc, char *argv[]) {
-  struct CODEC2 *codec2;
-  FILE *fin;
-  FILE *fout;
-  short *buf;
-  unsigned char *bits;
-  int nsam, nbit, i, r;
-  int errno = -1;
-  for (i = 0; i < 10; i++) {
-    r = codec2_rand();
-    printf("[%d] r = %d\n", i, r);
-  }
-
-  if (argc != 3) {
-    printf("usage: %s InputRawSpeechFile OutputRawSpeechFile\n", argv[0]);
-    exit(1);
-  }
-
-  if ((fin = fopen(argv[1], "rb")) == NULL) {
-    fprintf(stderr, "Error opening input speech file: %s: %s.\n",
-            argv[1], strerror(errno));
-    exit(1);
-  }
-
-  if ((fout = fopen(argv[2], "wb")) == NULL) {
-    fprintf(stderr, "Error opening output speech file: %s: %s.\n",
-            argv[2], strerror(errno));
-    exit(1);
-  }
+int main(int argc, char *argv[])
+{
+    struct CODEC2 *codec2;
+    FILE *fin;
+    FILE *fout;
+    short *buf;
+    unsigned char *bits;
+    int nsam, nbit, i, r;
+    int errno = -1;
+    for (i = 0; i < 10; i++) {
+        r = codec2_rand();
+        printf("[%d] r = %d\n", i, r);
+    }
+
+    if (argc != 3) {
+        printf("usage: %s InputRawSpeechFile OutputRawSpeechFile\n", argv[0]);
+        exit(1);
+    }
+
+    if ((fin = fopen(argv[1], "rb")) == NULL) {
+        fprintf(stderr, "Error opening input speech file: %s: %s.\n",
+                argv[1], strerror(errno));
+        exit(1);
+    }
+
+    if ((fout = fopen(argv[2], "wb")) == NULL) {
+        fprintf(stderr, "Error opening output speech file: %s: %s.\n",
+                argv[2], strerror(errno));
+        exit(1);
+    }
 
 #ifdef DUMP
-  dump_on("c2demo");
+    dump_on("c2demo");
 #endif
 
-  /* Note only one set of Codec 2 states is required for an encoder
-     and decoder pair. */
+    /* Note only one set of Codec 2 states is required for an encoder
+       and decoder pair. */
 
-  codec2 = codec2_create(CODEC2_MODE_WB);
-  nsam = C2WB_SPERF; //codec2_samples_per_frame(codec2);
-  buf = (short*) malloc(nsam * sizeof (short));
-  nbit = C2WB_BPERF; //codec2_bits_per_frame(codec2);
-  bits = (unsigned char*) malloc(nbit * sizeof (char));
+    codec2 = codec2_create(CODEC2_MODE_WB);
+    nsam = C2WB_SPERF; //codec2_samples_per_frame(codec2);
+    buf = (short*) malloc(nsam * sizeof (short));
+    nbit = C2WB_BPERF; //codec2_bits_per_frame(codec2);
+    bits = (unsigned char*) malloc(nbit * sizeof (char));
 
-  while (fread(buf, sizeof (short), nsam, fin) == (size_t) nsam) {
-    codec2_encode_wb(codec2, bits, buf);
-    codec2_decode_wb(codec2, buf, bits);
-    fwrite(buf, sizeof (short), nsam, fout);
-  }
+    while (fread(buf, sizeof (short), nsam, fin) == (size_t) nsam) {
+        codec2_encode_wb(codec2, bits, buf);
+        codec2_decode_wb(codec2, buf, bits);
+        fwrite(buf, sizeof (short), nsam, fout);
+    }
 
-  free(buf);
-  free(bits);
-  codec2_destroy(codec2);
+    free(buf);
+    free(bits);
+    codec2_destroy(codec2);
 
-  fclose(fin);
-  fclose(fout);
+    fclose(fin);
+    fclose(fout);
 
-  return 0;
+    return 0;
 }
 #endif
 
-void codec2_encode_wb(struct CODEC2 *c2, unsigned char * bits, short speech[]) {
-  return;
+void codec2_encode_wb(struct CODEC2 *c2, unsigned char * bits, short speech[])
+{
+    return;
 }
 
-void codec2_decode_wb(struct CODEC2 *c2, short speech[], const unsigned char * bits) {
-  return;
+void codec2_decode_wb(struct CODEC2 *c2, short speech[], const unsigned char * bits)
+{
+    return;
 }
 
-void calculate_Am_freqs_kHz(float Wo, int L, float Am_freqs_kHz[L]) {
-  //(1:L)*Wo*4/pi;
-  int i;
-  for (i = 0; i < L; i++) {
-    Am_freqs_kHz[i] = i * Wo * 4 / PI;
-  }
+void calculate_Am_freqs_kHz(float Wo, int L, float Am_freqs_kHz[L])
+{
+    //(1:L)*Wo*4/pi;
+    int i;
+    for (i = 0; i < L; i++) {
+        Am_freqs_kHz[i] = i * Wo * 4 / PI;
+    }
 }
 
 
 // Adapted from newamp.m
 
-void resample_const_rate_f_mel(C2CONST *c2const, MODEL * model, float K, float* rate_K_surface, float* rate_K_sample_freqs_kHz) {
-  mel_sample_freqs_kHz(rate_K_sample_freqs_kHz, K, 100, 0.95 * c2const->Fs / 2);
-  resample_const_rate_f(c2const, model, rate_K_surface, rate_K_sample_freqs_kHz, K);
+void resample_const_rate_f_mel(C2CONST *c2const, MODEL * model, float K, float* rate_K_surface, float* rate_K_sample_freqs_kHz)
+{
+    mel_sample_freqs_kHz(rate_K_sample_freqs_kHz, K, 100, 0.95 * c2const->Fs / 2);
+    resample_const_rate_f(c2const, model, rate_K_surface, rate_K_sample_freqs_kHz, K);
 }
 
 
 // Updates rate_K_vec with the corrected values
 // function [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs)
 
-void correct_rate_K_vec(MODEL *model, float rate_K_vec[], float rate_K_sample_freqs_kHz[], float Am_freqs_kHz[], float orig_AmdB[], int K, float Wo, int L, int Fs, float rate_K_vec_corrected[]) {
+void correct_rate_K_vec(MODEL *model, float rate_K_vec[], float rate_K_sample_freqs_kHz[], float Am_freqs_kHz[], float orig_AmdB[], int K, float Wo, int L, int Fs, float rate_K_vec_corrected[])
+{
 
-  /*
-      % aliasing correction --------------------------------------
+    /*
+        % aliasing correction --------------------------------------
 
-      % The mel sample rate decreases as frequency increases. Look for
-      % any regions above 1000Hz where we have missed definition of a
-      % spectral peak (formant) due to aliasing.  Adjust the rate K
-      % sample levels to restore peaks.  Theory is that correct
-      % definition of a formant is less important than the frequency of
-      % the formant.  As long as we define a formant in that general
-      % frequency area it will sound OK.
-   */
+        % The mel sample rate decreases as frequency increases. Look for
+        % any regions above 1000Hz where we have missed definition of a
+        % spectral peak (formant) due to aliasing.  Adjust the rate K
+        % sample levels to restore peaks.  Theory is that correct
+        % definition of a formant is less important than the frequency of
+        % the formant.  As long as we define a formant in that general
+        % frequency area it will sound OK.
+     */
 
 
-  // this is passed in
-  //Am_freqs_kHz = (1:L)*Wo*4/pi;    
+    // this is passed in
+    //Am_freqs_kHz = (1:L)*Wo*4/pi;    
 
-  //% Lets see where we have made an error
-  float error[MAX_AMP + 1];
-  //float orig_error[MAX_AMP+1];
-  float AmdB_[MAX_AMP + 1];
+    //% Lets see where we have made an error
+    float error[MAX_AMP + 1];
+    //float orig_error[MAX_AMP+1];
+    float AmdB_[MAX_AMP + 1];
 
-  float nasty_error_freq;
-  float rate_K_prev_sample_kHz;
-  float rate_K_next_sample_kHz;
+    float nasty_error_freq;
+    float rate_K_prev_sample_kHz;
+    float rate_K_next_sample_kHz;
 
-  int Ncorrections = 3; //% maximum number of rate K samples to correct
-  int error_thresh = 3; //% only worry about errors larger than thresh    
+    int Ncorrections = 3; //% maximum number of rate K samples to correct
+    int error_thresh = 3; //% only worry about errors larger than thresh    
 
-  int closest_k = 0;
-  int i;
+    int closest_k = 0;
+    int i;
 
-  // regenerate the AmdB values from the updated model
-  for (i = 0; i < MAX_AMP + 1; i++) {
-    AmdB_[i] = 20 * log10(model->A[i]);
-  }
-
-  // calculate error between original AmdB and the new values
-  for (i = 0; i < MAX_AMP + 1; i++) {
-    int a = orig_AmdB[i] - AmdB_[i];
-    error[i] = a;
-    //  orig_error[i] = a;
-  }
-
-
-
-  //% first 1000Hz is densely sampled so ignore
-  int start_m = floor(L * 1000 / (Fs / 2));
-  for (i = 0; i < start_m; i++) {
-    error[i] = 0;
-  }
-  //start_m = floor(L*1000/(Fs/2));
-  //error(1:start_m) = 0;  % first 1000Hz is densly sampled so ignore
-
-  // ignore these
-  //float nasty_error_m_log = []; 
-  //float nasty_error_log = [];
+    // regenerate the AmdB values from the updated model
+    for (i = 0; i < MAX_AMP + 1; i++) {
+        AmdB_[i] = 20 * log10(model->A[i]);
+    }
 
+    // calculate error between original AmdB and the new values
+    for (i = 0; i < MAX_AMP + 1; i++) {
+        int a = orig_AmdB[i] - AmdB_[i];
+        error[i] = a;
+        //  orig_error[i] = a;
+    }
 
-  // could probably memcpy this, but I need to think about it...
-  // TODO copy??
-  for (i = 0; i < K; i++) {
-    //rate_K_vec_corrected = rate_K_vec
-    rate_K_vec_corrected[i] = rate_K_vec[i];
-  }
 
-  for (i = 0; i < Ncorrections; i++) {
 
-    // [mx mx_m] = max(error);
-    int mx_m = 0;
-    int mx = 0;
-    int m;
-    for (m = 0; m < start_m; m++) {
-      if (error[m] > mx) {
-        mx_m = m;
-        mx = error[m];
-      }
+    //% first 1000Hz is densely sampled so ignore
+    int start_m = floor(L * 1000 / (Fs / 2));
+    for (i = 0; i < start_m; i++) {
+        error[i] = 0;
     }
+    //start_m = floor(L*1000/(Fs/2));
+    //error(1:start_m) = 0;  % first 1000Hz is densly sampled so ignore
 
-    if (mx > error_thresh) {
-
-      // ignore these
-      //nasty_error_log = [nasty_error_log mx];
-      //nasty_error_m_log = [nasty_error_m_log mx_m];
-
-      //% find closest rate K sample to nasty error
+    // ignore these
+    //float nasty_error_m_log = []; 
+    //float nasty_error_log = [];
 
-      nasty_error_freq = (float) mx_m * Wo * Fs / (2 * PI * 1000);
 
-      //[tmp closest_k] = min(abs(rate_K_sample_freqs_kHz - nasty_error_freq));
-      closest_k = -1;
-      float cka = K;
-      int k;
-      for (k = 0; k < K; k++) {
-        float a = abs(rate_K_sample_freqs_kHz[k] - nasty_error_freq);
+    // could probably memcpy this, but I need to think about it...
+    // TODO copy??
+    for (i = 0; i < K; i++) {
+        //rate_K_vec_corrected = rate_K_vec
+        rate_K_vec_corrected[i] = rate_K_vec[i];
+    }
 
-        if (closest_k == -1 || a < cka) {
-          closest_k = k;
-          cka = a;
-        }
-      }
-
-      // this sets the value from the original AmdB value set for the model prior to recalculation
-      // TODO - check this is correct
-      rate_K_vec_corrected[closest_k] = orig_AmdB[mx_m];
-
-      //% zero out error in this region and look for another large error region
-
-      // here I'm assuming that the required result is a zero based index
-      // rather than an Octave 1 based index
-      // closest_k is already 0 based in this code
-      // TODO - check this
-      k = (0 > closest_k - 1 ? 0 : closest_k - 1);
-      //k = max(1, closest_k-1); 
-
-      rate_K_prev_sample_kHz = rate_K_sample_freqs_kHz[k];
-
-
-      // again, assuming the required result is 0 rather than 1 based
-      // TODO - check 
-      k = (K - 1 < closest_k + 1 ? K : closest_k + 1);
-      //k = min(K, closest_k+1); 
-
-      rate_K_next_sample_kHz = rate_K_sample_freqs_kHz[k];
-
-      int st = -1;
-      int st_m = 0;
-      int en = -1;
-      int en_m = 0;
-      for (m = 0; m < C2WB_K - 1; m++) {
-        //[tmp st_m] = min(abs(Am_freqs_kHz - rate_K_prev_sample_kHz));
-        //[tmp en_m] = min(abs(Am_freqs_kHz - rate_K_next_sample_kHz));
-        int pa = abs(Am_freqs_kHz[m] - rate_K_prev_sample_kHz);
-        int na = abs(abs(Am_freqs_kHz[m] - rate_K_next_sample_kHz));
-        if (st == -1 || pa < st) {
-          st = pa;
-          st_m = m;
+    for (i = 0; i < Ncorrections; i++) {
+
+        // [mx mx_m] = max(error);
+        int mx_m = 0;
+        int mx = 0;
+        int m;
+        for (m = 0; m < start_m; m++) {
+            if (error[m] > mx) {
+                mx_m = m;
+                mx = error[m];
+            }
         }
-        if (en == -1 || na < en) {
-          en = na;
-          en_m = m;
+
+        if (mx > error_thresh) {
+
+            // ignore these
+            //nasty_error_log = [nasty_error_log mx];
+            //nasty_error_m_log = [nasty_error_m_log mx_m];
+
+            //% find closest rate K sample to nasty error
+
+            nasty_error_freq = (float) mx_m * Wo * Fs / (2 * PI * 1000);
+
+            //[tmp closest_k] = min(abs(rate_K_sample_freqs_kHz - nasty_error_freq));
+            closest_k = -1;
+            float cka = K;
+            int k;
+            for (k = 0; k < K; k++) {
+                float a = abs(rate_K_sample_freqs_kHz[k] - nasty_error_freq);
+
+                if (closest_k == -1 || a < cka) {
+                    closest_k = k;
+                    cka = a;
+                }
+            }
+
+            // this sets the value from the original AmdB value set for the model prior to recalculation
+            // TODO - check this is correct
+            rate_K_vec_corrected[closest_k] = orig_AmdB[mx_m];
+
+            //% zero out error in this region and look for another large error region
+
+            // here I'm assuming that the required result is a zero based index
+            // rather than an Octave 1 based index
+            // closest_k is already 0 based in this code
+            // TODO - check this
+            k = (0 > closest_k - 1 ? 0 : closest_k - 1);
+            //k = max(1, closest_k-1); 
+
+            rate_K_prev_sample_kHz = rate_K_sample_freqs_kHz[k];
+
+
+            // again, assuming the required result is 0 rather than 1 based
+            // TODO - check 
+            k = (K - 1 < closest_k + 1 ? K : closest_k + 1);
+            //k = min(K, closest_k+1); 
+
+            rate_K_next_sample_kHz = rate_K_sample_freqs_kHz[k];
+
+            int st = -1;
+            int st_m = 0;
+            int en = -1;
+            int en_m = 0;
+            for (m = 0; m < C2WB_K - 1; m++) {
+                //[tmp st_m] = min(abs(Am_freqs_kHz - rate_K_prev_sample_kHz));
+                //[tmp en_m] = min(abs(Am_freqs_kHz - rate_K_next_sample_kHz));
+                int pa = abs(Am_freqs_kHz[m] - rate_K_prev_sample_kHz);
+                int na = abs(abs(Am_freqs_kHz[m] - rate_K_next_sample_kHz));
+                if (st == -1 || pa < st) {
+                    st = pa;
+                    st_m = m;
+                }
+                if (en == -1 || na < en) {
+                    en = na;
+                    en_m = m;
+                }
+            }
+            if (closest_k == K)
+                en_m = L;
+
+            for (m = st_m; m < en_m; m++) {
+                error[m] = 0;
+            }
+            //error(st_m:en_m) = 0;
         }
-      }
-      if (closest_k == K)
-        en_m = L;
-
-      for (m = st_m; m < en_m; m++) {
-        error[m] = 0;
-      }
-      //error(st_m:en_m) = 0;
     }
-  }
 
 
 
 }
 
-void diff_de(int rows, int cols, float D[rows][cols], float E[rows][cols], float diff_de[rows][cols]) {
+void diff_de(int rows, int cols, float D[rows][cols], float E[rows][cols], float diff_de[rows][cols])
+{
 
-  int col, row;
-  for (col = 0; col < cols; col++) {
-    for (row = 0; row < rows; row++) {
-      diff_de[row][col] = D[row][col] - E[row][col];
+    int col, row;
+    for (col = 0; col < cols; col++) {
+        for (row = 0; row < rows; row++) {
+            diff_de[row][col] = D[row][col] - E[row][col];
+        }
     }
-  }
 
 }
 
-float mean(float data[], int n) {
-  float sum = 0.0;
-  int i;
-  for (i = 0; i < n; ++i) {
-    sum += data[i];
-  }
-
-  return sum / n;
-}
-
-void array_col_to_row(int rows, int cols, float data[rows][cols], int col, float res[]) {
-  int row;
-  for (row = 0; row < rows; row++) {
-    res[row] = data[row][col];
-  }
-}
+float mean(float data[], int n)
+{
+    float sum = 0.0;
+    int i;
+    for (i = 0; i < n; ++i) {
+        sum += data[i];
+    }
 
-float std(float data[], int rows) {
-  float standardDeviation = 0.0;
-  int i;
-  for (i = 0; i < rows; ++i) {
-    standardDeviation += pow(data[i] - mean(data, rows), 2);
-  }
-  return sqrt(standardDeviation / rows);
+    return sum / n;
 }
 
-void std_on_cols(int rows, int cols, float data[rows][cols], float res[]) {
-  float row_from_col[cols];
-  int col;
-  for (col = 0; col < cols; col++) {
-    array_col_to_row(rows, cols, data, col, row_from_col);
-    res[col] = std(row_from_col, rows);
-  }
+void array_col_to_row(int rows, int cols, float data[rows][cols], int col, float res[])
+{
+    int row;
+    for (row = 0; row < rows; row++) {
+        res[row] = data[row][col];
+    }
 }
 
-float mean_std_diff_de(int rows, int cols, float D[rows][cols], float E[rows][cols]) {
-  float matrix_diff_de[rows][cols];
-  float std_array[cols];
-  diff_de(rows, cols, D, E, matrix_diff_de);
-  std_on_cols(rows, cols, matrix_diff_de, std_array);
-  return mean(std_array, cols);
+float std(float data[], int rows)
+{
+    float standardDeviation = 0.0;
+    int i;
+    for (i = 0; i < rows; ++i) {
+        standardDeviation += pow(data[i] - mean(data, rows), 2);
+    }
+    return sqrt(standardDeviation / rows);
 }
 
-//TODO : everything
-
-void dct2(int rows, int cols, float matrix[], float res[rows][cols]) {
-  int r,c;
-  for (r = 0; r < rows; r++) {
-    for (c = 0; c < cols; c++) {
-      res[r][c] = (float) rand() / RAND_MAX;
+void std_on_cols(int rows, int cols, float data[rows][cols], float res[])
+{
+    float row_from_col[cols];
+    int col;
+    for (col = 0; col < cols; col++) {
+        array_col_to_row(rows, cols, data, col, row_from_col);
+        res[col] = std(row_from_col, rows);
     }
-  }
-  /*
-     
-    if m == 1
-      y = dct (x.', n).';
-    elseif n == 1
-      y = dct (x, m);
-    else
-      y = dct (dct (x, m).', n).';
-    endif
-   */
+}
 
+float mean_std_diff_de(int rows, int cols, float D[rows][cols], float E[rows][cols])
+{
+    float matrix_diff_de[rows][cols];
+    float std_array[cols];
+    diff_de(rows, cols, D, E, matrix_diff_de);
+    std_on_cols(rows, cols, matrix_diff_de, std_array);
+    return mean(std_array, cols);
 }
 
-//TODO : everything
 
-void idct2(int rows, int cols, float inrks[rows][cols], float rate_K_surface_block[]) {
-  int K = C2WB_K;
-  int k;
-  for (k = 0; k < K; k++) {
-    rate_K_surface_block[k] = (float) rand() / RAND_MAX;
-  }
-}
 //% Encode/decoder a 160ms block of model parameters
 //% TODO: (i) quantisation of DCT coeffs (ii) break into separate encoder and decoder functions
 //
 //function [model_block_ dct2_sd qn rate_K_surface_block rate_K_surface_block_] = wideband_enc_dec(model_block, rmap, cmap)
 
 void wideband_enc_dec(C2CONST *c2const, int n_block_frames, MODEL model_block[n_block_frames], WIDEBAND_MAP * wb_map,
-        MODEL model_block_[n_block_frames], float dct2_sd[n_block_frames], int * p_qn, float rate_K_surface_block[n_block_frames][C2WB_K], float rate_K_surface_block_[n_block_frames][C2WB_K]) {
-
-  
-  //printf("starting wideband_enc_dec with model_block[0]->L = %d\n", model_block[0].L);
-  //    c2wideband_const;
-  //
-  //    sim_quant = 1; % used to simulate quantisation, set to 1,2,4, etc
-  //    dist_dB   = 2; % use enough coefficients to get this distortion ond DCT coeffs
-  int sim_quant = 1;
-  float dist_dB = 2;
+                      MODEL model_block_[n_block_frames], float dct2_sd[n_block_frames], int * p_qn, float rate_K_surface_block[n_block_frames][C2WB_K], float rate_K_surface_block_[n_block_frames][C2WB_K])
+{
 
-  //    int Fs = c2const->Fs; 
-  //    int L;   
-  //    float Wo;
 
-  int K = C2WB_K;
-  int Nt = C2WB_NT;
-
-  int rows = Nt;
-  int cols = K;
-  int dec = C2WB_DEC;
-  //printf("starting iteration\n");
-  // iterate through the frames in the block
-  int f;
-  for (f = 0; f < n_block_frames; f++) {
-    //printf("processing frame %d of %d\n", f, n_block_frames);
-
-    float rate_K_sample_freqs_kHz[C2WB_K];
-    //MODEL * model = &model_block[f];
-
-
-    //    % Resample variable rate L vectors to fixed length rate K.  We have
-    //    % left high end correction out for now, this is less of an issue
-    //    % with a higher K
-
-
-    //    [rate_K_surface_block rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model_block, K, Fs);
-    //printf("resample_const_rate_f_mel\n");
-    resample_const_rate_f_mel(c2const, &model_block[f], K, rate_K_surface_block[f], rate_K_sample_freqs_kHz);
-
-    //    % decimate down to 20ms time resolution, and DCT
-
-    //    D = dct2(rate_K_surface_block(1:dec:Nt*dec,:));
-
-    float D[rows][cols];
-    float E[rows][cols];
-    //printf("dct2\n");
-    dct2(rows, cols, rate_K_surface_block[f], D);
-
-    //    % So D is the 2D block of DCT coeffs at the encoder.  We want to
-    //    % create a quantised version at the "decoder" E.  This loop copies
-    //    % DCTs coeffs from D to E, until we get beneath a distortion
-    //    % threshold.
-    //
-    //    % This is essentially variable rate quantisation, but gives us
-    //    % some idea of the final bit rate.  In practice we will also need
-    //    % to constrain the total number of bits (ie bit rate), and
-    //    % quantise each coefficient.
+    //printf("starting wideband_enc_dec with model_block[0]->L = %d\n", model_block[0].L);
+    //    c2wideband_const;
     //
-    //    % Turns out than mean SD (across many blocks/frames) is about the
-    //    % same in the DCT domain as the rate K domain.  So we can just
-    //    % measure MSE between D and E to estimate mean SD of the rate K
-    //    % vectors after quantisation.
-
-    //    E = mapped = zeros(Nt,K);                
-    int r,c;
-    for (r = 0; r < rows; r++) {
-      memset(E[r], '\0', cols * sizeof (float));
-    }
+    //    sim_quant = 1; % used to simulate quantisation, set to 1,2,4, etc
+    //    dist_dB   = 2; % use enough coefficients to get this distortion ond DCT coeffs
+    int sim_quant = 1;
+    float dist_dB = 2;
+
+    //    int Fs = c2const->Fs; 
+    //    int L;   
+    //    float Wo;
+
+    int K = C2WB_K;
+    int Nt = C2WB_NT;
+
+    int rows = Nt;
+    int cols = K;
+    int dec = C2WB_DEC;
+    //printf("starting iteration\n");
+    // iterate through the frames in the block
+    int f;
+    for (f = 0; f < n_block_frames; f++) {
+        //printf("processing frame %d of %d\n", f, n_block_frames);
+
+        float rate_K_sample_freqs_kHz[C2WB_K];
+        //MODEL * model = &model_block[f];
+
+
+        //    % Resample variable rate L vectors to fixed length rate K.  We have
+        //    % left high end correction out for now, this is less of an issue
+        //    % with a higher K
+
+
+        //    [rate_K_surface_block rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model_block, K, Fs);
+        //printf("resample_const_rate_f_mel\n");
+        resample_const_rate_f_mel(c2const, &model_block[f], K, rate_K_surface_block[f], rate_K_sample_freqs_kHz);
+
+        //    % decimate down to 20ms time resolution, and DCT
+
+        //    D = dct2(rate_K_surface_block(1:dec:Nt*dec,:));
+
+        float D[rows][cols];
+        float E[rows][cols];
+        //printf("dct2\n");
+        dct2(rows, cols, &rate_K_surface_block[f], D);
+
+        //    % So D is the 2D block of DCT coeffs at the encoder.  We want to
+        //    % create a quantised version at the "decoder" E.  This loop copies
+        //    % DCTs coeffs from D to E, until we get beneath a distortion
+        //    % threshold.
+        //
+        //    % This is essentially variable rate quantisation, but gives us
+        //    % some idea of the final bit rate.  In practice we will also need
+        //    % to constrain the total number of bits (ie bit rate), and
+        //    % quantise each coefficient.
+        //
+        //    % Turns out than mean SD (across many blocks/frames) is about the
+        //    % same in the DCT domain as the rate K domain.  So we can just
+        //    % measure MSE between D and E to estimate mean SD of the rate K
+        //    % vectors after quantisation.
+
+        //    E = mapped = zeros(Nt,K);                
+        int r, c;
+        for (r = 0; r < rows; r++) {
+            memset(E[r], '\0', cols * sizeof (float));
+        }
 
-    //    qn = 0;
-    *p_qn = 0;
+        //    qn = 0;
+        *p_qn = 0;
 
-    //    adct2_sd = mean(std(D-E));
+        //    adct2_sd = mean(std(D-E));
 
-    float adct2_sd;
+        float adct2_sd;
 
 
-    adct2_sd = mean_std_diff_de(rows, cols, D, E);
+        adct2_sd = mean_std_diff_de(rows, cols, D, E);
 
-    //printf("while adct2_sd > dist_dB (dist_dB=%f)\n", dist_dB);
-    //    while adct2_sd > dist_dB
-    //TODO : validate the addition of  *p_qn limit check is correct
-    while (adct2_sd > dist_dB && *p_qn < C2WB_NT) {
+        //printf("while adct2_sd > dist_dB (dist_dB=%f)\n", dist_dB);
+        //    while adct2_sd > dist_dB
+        //TODO : validate the addition of  *p_qn limit check is correct
+        while (adct2_sd > dist_dB && *p_qn < C2WB_NT) {
 
-      //      qn++;                        
-      // --- this has moved to the end to cope with 0 base
-      assert(*p_qn < C2WB_NT);
-      int rmapqn = wb_map->rmap[*p_qn];
-      int cmapqn = wb_map->cmap[*p_qn];
-      assert(rmapqn < rows);
-      assert(cmapqn < cols);
-      //      E(rmap(qn), cmap(qn)) = sim_quant*round(D(rmap(qn), cmap(qn))/sim_quant);
-      E[rmapqn][cmapqn] = sim_quant * round(D[rmapqn][cmapqn] / sim_quant);
-      //      adct2_sd = mean(std(D-E));
-      adct2_sd = mean_std_diff_de(rows, cols, D, E);
+            //      qn++;                        
+            // --- this has moved to the end to cope with 0 base
+            assert(*p_qn < C2WB_NT);
+            int rmapqn = wb_map->rmap[*p_qn];
+            int cmapqn = wb_map->cmap[*p_qn];
+            assert(rmapqn < rows);
+            assert(cmapqn < cols);
+            //      E(rmap(qn), cmap(qn)) = sim_quant*round(D(rmap(qn), cmap(qn))/sim_quant);
+            E[rmapqn][cmapqn] = sim_quant * round(D[rmapqn][cmapqn] / sim_quant);
+            //      adct2_sd = mean(std(D-E));
+            adct2_sd = mean_std_diff_de(rows, cols, D, E);
 
-      //printf("qn %d %f\n", *p_qn, adct2_sd);
+            //printf("qn %d %f\n", *p_qn, adct2_sd);
 
-      (*p_qn)++;
-      //      %printf("qn %d %f\n", qn, adct2_sd);
-      //    end
-    }
+            (*p_qn)++;
+            //      %printf("qn %d %f\n", qn, adct2_sd);
+            //    end
+        }
 
-    //    % note neat trick to interpolate to 10ms frames despite dec to 20ms, this means
-    //    % we don't need a separate decode side interpolator.
+        //    % note neat trick to interpolate to 10ms frames despite dec to 20ms, this means
+        //    % we don't need a separate decode side interpolator.
 
-    //    dct2_sd = mean(std(D-E));
+        //    dct2_sd = mean(std(D-E));
 
 
-    dct2_sd[f] = mean_std_diff_de(rows, cols, D, E);
+        dct2_sd[f] = mean_std_diff_de(rows, cols, D, E);
 
 
-    //    rate_K_surface_block_ = idct2([sqrt(dec)*E; zeros(Nt*(dec-1), K)]);
-    float inrks[rows * dec][K];
+        //    rate_K_surface_block_ = idct2([sqrt(dec)*E; zeros(Nt*(dec-1), K)]);
+        float inrks[rows * dec][K];
 
-    //printf("setting inrks\n");
-    for (r = 0; r < rows; r++) {
-      for (c = 0; c < cols; c++) {
-        inrks[r][c] = sqrt(dec) * E[r][c];
-      }
-    }
-    for (r = 0; r < Nt * (dec - 1); r++) {
-      for (c = 0; c < K; c++) {
-        inrks[r + rows][c] = 0;
-      }
-    }
+        //printf("setting inrks\n");
+        for (r = 0; r < rows; r++) {
+            for (c = 0; c < cols; c++) {
+                inrks[r][c] = sqrt(dec) * E[r][c];
+            }
+        }
+        for (r = 0; r < Nt * (dec - 1); r++) {
+            for (c = 0; c < K; c++) {
+                inrks[r + rows][c] = 0;
+            }
+        }
 
 
-    //TODO ???
-    //= [sqrt(dec)*E; zeros(Nt*(dec-1), K)];
-    //printf("idct2\n");
-    idct2(rows, cols, inrks, rate_K_surface_block_[f]);
+        //TODO ???
+        //= [sqrt(dec)*E; zeros(Nt*(dec-1), K)];
+        //printf("idct2\n");
+        idct2(rows, cols, inrks, &rate_K_surface_block_[f]);
 
-    //    model_block_ = resample_rate_L(model_block, rate_K_surface_block_, rate_K_sample_freqs_kHz, Fs);        
-    //printf("resample_rate_L\n");
-    resample_rate_L(c2const, &model_block[f], rate_K_surface_block_[f], rate_K_sample_freqs_kHz, K);
-    //endfunction 
-  }
+        //    model_block_ = resample_rate_L(model_block, rate_K_surface_block_, rate_K_sample_freqs_kHz, Fs);        
+        //printf("resample_rate_L\n");
+        resample_rate_L(c2const, &model_block[f], rate_K_surface_block_[f], rate_K_sample_freqs_kHz, K);
+        //endfunction 
+    }
 }
 
-void setup_map(WIDEBAND_MAP * wb_map, int Nt, int K) {
+void setup_map(WIDEBAND_MAP * wb_map, int Nt, int K)
+{
 
-  /*
-    % map that defines order we read out and quantise DCT coeffs
-    % TODO: for C port we need an Octave function to write Map to a C
-    % include file
+    /*
+      % map that defines order we read out and quantise DCT coeffs
+      % TODO: for C port we need an Octave function to write Map to a C
+      % include file
 
-    map = load("c2wideband_map");
+      map = load("c2wideband_map");
 
-    % create arrays to reverse map quantiser_num to r,c Luts
+      % create arrays to reverse map quantiser_num to r,c Luts
 
-    rmap = cmap = zeros(1,Nt*K);
-    for r=1:Nt
-      for c=1:K
-         quantiser_num = map(r,c);
-         rmap(quantiser_num) = r;
-         cmap(quantiser_num) = c;
+      rmap = cmap = zeros(1,Nt*K);
+      for r=1:Nt
+        for c=1:K
+           quantiser_num = map(r,c);
+           rmap(quantiser_num) = r;
+           cmap(quantiser_num) = c;
+        end
       end
-    end
 
-   */
+     */
 
-  //printf("setup_map(wb_map, Nt, K)");
+    //printf("setup_map(wb_map, Nt, K)");
     int quantiser_num, r, c;
 
-  memset(wb_map->rmap, '\0', Nt * K * sizeof *wb_map->rmap);
-  memset(wb_map->cmap, '\0', Nt * K * sizeof *wb_map->cmap);
+    memset(wb_map->rmap, '\0', Nt * K * sizeof *wb_map->rmap);
+    memset(wb_map->cmap, '\0', Nt * K * sizeof *wb_map->cmap);
 
-  for (r = 0; r < Nt; r++) {
-    for (c = 0; c < K; c++) {
-      quantiser_num = c2wideband_map[r][c];
-      wb_map->rmap[quantiser_num] = r;
-      wb_map->cmap[quantiser_num] = c;
+    for (r = 0; r < Nt; r++) {
+        for (c = 0; c < K; c++) {
+            quantiser_num = c2wideband_map[r][c];
+            wb_map->rmap[quantiser_num] = r;
+            wb_map->cmap[quantiser_num] = c;
+        }
     }
-  }
 
 }
 
@@ -592,122 +577,124 @@ function [model_ rate_K_surface] = experiment_rate_K_dct2(model, plots=1)
 
 //TODO - review and produce a unit test for this
 // this is really just a draft at this point
-void experiment_rate_K_dct2(C2CONST *c2const, MODEL model_frames[], const int total_frames) {
 
-  printf("experiment_rate_K_dct2 with frames: %d\n", total_frames);
-  //  newamp;
-  //  c2wideband_const;
-  //  [frames nc] = size(model);
+void experiment_rate_K_dct2(C2CONST *c2const, MODEL model_frames[], const int total_frames)
+{
+
+    printf("experiment_rate_K_dct2 with frames: %d\n", total_frames);
+    //  newamp;
+    //  c2wideband_const;
+    //  [frames nc] = size(model);
+
+    const int n_block_frames = C2WB_NT * C2WB_DEC;
+
+    const int K = C2WB_K;
+    const int Nt = C2WB_NT;
+    const int dec = C2WB_DEC;
+    const int Tf = C2WB_TF;
+    WIDEBAND_MAP wb_map;
+
+
+    //  % break into blocks of (Nt time samples) x (K freq samples)
 
-  const int n_block_frames = C2WB_NT * C2WB_DEC;
+    // Nblocks = floor(frames/(Nt*dec));
 
-  const int K = C2WB_K;
-  const int Nt = C2WB_NT;
-  const int dec = C2WB_DEC;    
-  const int Tf = C2WB_TF;
-  WIDEBAND_MAP wb_map;
+    const int Nblocks = floor(total_frames / n_block_frames);
 
+    // number of frames that can be processed (if the final block is not a full set of frames)
+    int frames;
 
-  //  % break into blocks of (Nt time samples) x (K freq samples)
+    frames = Nblocks * n_block_frames;
 
-  // Nblocks = floor(frames/(Nt*dec));
-  
-  const int Nblocks = floor(total_frames/n_block_frames);
+    printf("total_frames: %d processable frames: %d Nblocks: %d\n", total_frames, frames, Nblocks);
 
-  // number of frames that can be processed (if the final block is not a full set of frames)
-  int frames;
-  
-  frames = Nblocks * n_block_frames;
-  
-  printf("total_frames: %d processable frames: %d Nblocks: %d\n", total_frames, frames, Nblocks);
-  
 
-  setup_map(&wb_map, Nt, K);
+    setup_map(&wb_map, Nt, K);
 
 
 
-  //% per-block processing ----------------------------------------------------
+    //% per-block processing ----------------------------------------------------
 
-  //  % init a bunch of output variables
+    //  % init a bunch of output variables
 
-  //  rate_K_surface_ = zeros(Nblocks*Nt*dec, K);  
-  // sumnz = zeros(1,Nblocks);
-  // dct2_sd = zeros(1,Nblocks);
-  // model_ = [];
+    //  rate_K_surface_ = zeros(Nblocks*Nt*dec, K);  
+    // sumnz = zeros(1,Nblocks);
+    // dct2_sd = zeros(1,Nblocks);
+    // model_ = [];
 
 
-  float rate_K_surface_block[total_frames][K]; // rate K vecs for each frame, form a surface that makes pretty graphs   
-  float rate_K_surface_block_[total_frames][K];
+    float rate_K_surface_block[total_frames][K]; // rate K vecs for each frame, form a surface that makes pretty graphs   
+    float rate_K_surface_block_[total_frames][K];
 
 #ifdef C2WB_PLOT    
-//  float rate_K_surface[total_frames][K];
-//  float rate_K_surface_[total_frames][K];
-//  MODEL model_[total_frames];
+    //  float rate_K_surface[total_frames][K];
+    //  float rate_K_surface_[total_frames][K];
+    //  MODEL model_[total_frames];
 #endif
 
-  float sumnz[Nblocks];
-  float dct2_sd[Nblocks];
-  int qn;
+    float sumnz[Nblocks];
+    float dct2_sd[Nblocks];
+    int qn;
 
-  MODEL model_block_[n_block_frames];
+    MODEL model_block_[n_block_frames];
 
-  //for n=1:Nblocks
-  // Step through the model in blocks of 16 frames
-  int n_block = 0;
-  int f;
-  for (f = 0; f < frames; f += n_block_frames) {
-    MODEL * model_block = &model_frames[f];
+    //for n=1:Nblocks
+    // Step through the model in blocks of 16 frames
+    int n_block = 0;
+    int f;
+    for (f = 0; f < frames; f += n_block_frames) {
+        MODEL * model_block = &model_frames[f];
 
-  
-    //st = (n-1)*dec*Nt+1; en = st + dec*Nt - 1;
-    //printf("st: %d en: %d\n", st, en);
-    // effectively handled by the iterations through the block in wideband_enc_dec
 
-    //[model_block_ adct2_sd qn rate_K_surface_block rate_K_surface_block_] = wideband_enc_dec(model(st:en,:), rmap, cmap);
-//    void wideband_enc_dec(C2CONST *c2const, int n_block_frames, MODEL model_block[n_block_frames], WIDEBAND_MAP * wb_map,
-//        MODEL model_block_[n_block_frames], float dct2_sd[n_block_frames], int * p_qn, float rate_K_surface_block[n_block_frames][C2WB_K], float rate_K_surface_block_[n_block_frames][C2WB_K]) {
+        //st = (n-1)*dec*Nt+1; en = st + dec*Nt - 1;
+        //printf("st: %d en: %d\n", st, en);
+        // effectively handled by the iterations through the block in wideband_enc_dec
 
-    wideband_enc_dec(c2const, n_block_frames, model_block, &wb_map,
-            model_block_, &dct2_sd[n_block], &qn, &rate_K_surface_block[f], &rate_K_surface_block_[f]);
+        //[model_block_ adct2_sd qn rate_K_surface_block rate_K_surface_block_] = wideband_enc_dec(model(st:en,:), rmap, cmap);
+        //    void wideband_enc_dec(C2CONST *c2const, int n_block_frames, MODEL model_block[n_block_frames], WIDEBAND_MAP * wb_map,
+        //        MODEL model_block_[n_block_frames], float dct2_sd[n_block_frames], int * p_qn, float rate_K_surface_block[n_block_frames][C2WB_K], float rate_K_surface_block_[n_block_frames][C2WB_K]) {
+
+        wideband_enc_dec(c2const, n_block_frames, model_block, &wb_map,
+                         model_block_, &dct2_sd[n_block], &qn, &rate_K_surface_block[f], &rate_K_surface_block_[f]);
 
 #ifdef C2WB_PLOT
-    //model_ = [model_; model_block_];
-    //% log these for plotting/development
-
-    //rate_K_surface(st:en,:) = rate_K_surface_block;
-    //rate_K_surface_(st:en,:) = rate_K_surface_block_;
-    
-//    for (int p = 0; p < n_block_frames; p++) {
-      
-//      model_[f+p] = model_block_[p];
-      
-//      for (int k = 0; k < K; k++) {
-//        rate_K_surface[f+p][k] = rate_K_surface_block[n_block + p][k];
-//        rate_K_surface_[f+p][k] = rate_K_surface_block_[n_block + p][k];
-//      }
-//    }
-    //dct2_sd(n) = adct2_sd;            
-    sumnz[n_block] = (float) qn;
-    // already handled in call
+        //model_ = [model_; model_block_];
+        //% log these for plotting/development
+
+        //rate_K_surface(st:en,:) = rate_K_surface_block;
+        //rate_K_surface_(st:en,:) = rate_K_surface_block_;
+
+        //    for (int p = 0; p < n_block_frames; p++) {
+
+        //      model_[f+p] = model_block_[p];
+
+        //      for (int k = 0; k < K; k++) {
+        //        rate_K_surface[f+p][k] = rate_K_surface_block[n_block + p][k];
+        //        rate_K_surface_[f+p][k] = rate_K_surface_block_[n_block + p][k];
+        //      }
+        //    }
+        //dct2_sd(n) = adct2_sd;            
+        sumnz[n_block] = (float) qn;
+        // already handled in call
 #endif
-    n_block++;
-    //end
-  }
+        n_block++;
+        //end
+    }
 
 
 #ifdef C2WB_PLOT
-  printf("average dct spectral distortion: %3.2f dB\n", mean(dct2_sd, Nblocks));
-  printf("mean number of coeffs/DCT: %3.2f/%d\n", mean(sumnz, Nblocks), Nt * K);
-  printf("coeffs/second: %3.2f\n", mean(sumnz, Nblocks) / (Nt * Tf * dec));
-  printf("bits/s: %3.2f\n", 2.9 * mean(sumnz, Nblocks) / (Nt * Tf * dec));
-
-  //TODO
-  //dist = std((rate_K_surface_(1:dec:Nblocks*Nt*dec,:) - rate_K_surface(1:dec:Nblocks*Nt*dec,:))');
-  //float dist = 
-  //if plots
-  //figure(1); clf; plot(dist); title('Rate K SD');
-  //printf("Rate K spectral distortion mean: %3.2f dB var: %3.2f\n", mean(dist), var(dist));
-  //end
+    printf("average dct spectral distortion: %3.2f dB\n", mean(dct2_sd, Nblocks));
+    printf("mean number of coeffs/DCT: %3.2f/%d\n", mean(sumnz, Nblocks), Nt * K);
+    printf("coeffs/second: %3.2f\n", mean(sumnz, Nblocks) / (Nt * Tf * dec));
+    printf("bits/s: %3.2f\n", 2.9 * mean(sumnz, Nblocks) / (Nt * Tf * dec));
+
+    //TODO
+    //dist = std((rate_K_surface_(1:dec:Nblocks*Nt*dec,:) - rate_K_surface(1:dec:Nblocks*Nt*dec,:))');
+    //float dist = 
+    //if plots
+    //figure(1); clf; plot(dist); title('Rate K SD');
+    //printf("Rate K spectral distortion mean: %3.2f dB var: %3.2f\n", mean(dist), var(dist));
+    //end
 #endif
 
 
diff --git a/codec2-dev/src/dct2.c b/codec2-dev/src/dct2.c
new file mode 100644 (file)
index 0000000..641bba4
--- /dev/null
@@ -0,0 +1,165 @@
+#include "dct2.h"
+
+/*
+ * 1-dimensional DCT, taking and returning only the real components of the result.
+ * The implementation relies on the real FFT used elsewhere within Codec 2, ensuring
+ * that any optimisations of the complex code are reused, rather than reimplemented.
+ */
+void dct(const int N, float y[N], float res[N])
+{
+
+    int i;
+    float y2[2 * N]; // input to FFT will be double the size of the input
+    COMP c[N + 1];
+    COMP phi[N];
+
+
+    codec2_fftr_cfg cfg = codec2_fftr_alloc(2 * N, 0, NULL, NULL);
+
+    for (i = 0; i < N; i++) {
+        y2[i] = y[i];
+        y2[N + i] = y[N - 1 - i];
+    }
+
+    codec2_fftr(cfg, y2, c);
+
+    for (i = 0; i < N; i++) {
+        float p;
+        if (i == 0)
+            p = 0;
+        else
+            p = -1 * PI * i / (2 * N);
+        phi[i] = comp_exp_j(p);
+        res[i] = cmult(phi[i], c[i]).real;
+    }
+
+}
+
+/* 
+ * 2-dimensional DCT, using and returning only the real components.
+ * Built on the 1-D DCT.
+ */
+
+void dct2(const int M, const int N, float y[M][N], float res[M][N])
+{
+    float a[M][N];
+
+    int i;
+    int j;
+    for (i = 0; i < M; i++) {
+        dct(N, y[i], a[i]);
+    }
+
+    float a_col[M];
+    float b_col[M];
+
+    for (j = 0; j < N; j++) {
+        for (i = 0; i < M; i++) {
+            a_col[i] = a[i][j];
+        }
+
+        dct(M, a_col, b_col);
+
+        for (i = 0; i < M; i++) {
+            res[i][j] = b_col[i];
+        }
+    }
+}
+
+//def idct(a):
+//    N = len(a)
+//    c = empty(N+1,complex)
+//
+//    phi = exp(1j*pi*arange(N)/(2*N))
+//    c[:N] = phi*a
+//    c[N] = 0.0
+//    return irfft(c)[:N]
+
+/* 
+ * Inverse DCT 1-dimensional.
+ */
+
+void idct(const int N, float a[N], float res[N])
+{
+    
+    int nfft = 2 * N ;
+    
+    // if nfft, the size of the inverse FFT output is odd, extend it by one
+    if(nfft & 1)
+        nfft++;
+    
+    int i;
+    COMP ac;
+    COMP c[N + 1];
+    COMP phi[N];
+    float res_fft[nfft];
+
+    codec2_fftr_cfg cfg = codec2_fftr_alloc(nfft, 1, NULL, NULL);
+
+    assert(cfg);
+    
+    for (i = 0; i < N; i++) {
+        float p;
+        p = PI * i / (2 * N);
+        phi[i] = comp_exp_j(p);
+        ac.real = a[i];
+        ac.imag = 0;
+        c[i] = cmult(phi[i], ac);
+    }
+
+    c[N].real = 0;
+    c[N].imag = 0;
+    
+    codec2_fftri(cfg, c, res_fft);
+    
+    // Scale the result
+    for(i=0;i<N; i++){
+        res[i] = res_fft[i] / nfft;
+    }
+
+
+}
+
+//def idct2(b):
+//    M = b.shape[0]
+//    N = b.shape[1]
+//    a = empty([M,N],float)
+//    y = empty([M,N],float)
+//
+//    for i in range(M):
+//        a[i,:] = idct(b[i,:])
+//    for j in range(N):
+//        y[:,j] = idct(a[:,j])
+//
+//    return y
+
+/*
+ *  Inverse DCT 2-dimensional
+ */
+void idct2(int M, int N, float y[M][N], float res[M][N])
+{
+
+    float a[M][N];
+
+    int i;
+    int j;
+    for (i = 0; i < M; i++) {
+        idct(N, y[i], a[i]);
+    }
+
+    float a_col[M];
+    float b_col[M];
+
+    for (j = 0; j < N; j++) {
+        for (i = 0; i < M; i++) {
+            a_col[i] = a[i][j];
+        }
+
+        idct(M, a_col, b_col);
+
+        for (i = 0; i < M; i++) {
+            res[i][j] = b_col[i];
+        }
+    }
+
+}
diff --git a/codec2-dev/src/dct2.h b/codec2-dev/src/dct2.h
new file mode 100644 (file)
index 0000000..2ecc90e
--- /dev/null
@@ -0,0 +1,21 @@
+/* 
+ * File:   dct2.h
+ * Author: Phil Ayres
+ *
+ * Created on 21 July 2017, 14:24
+ */
+
+#ifndef DCT2_H
+#define        DCT2_H
+
+#include "codec2_fft.h"
+#include "comp.h"
+#include "comp_prim.h"
+
+void dct(const int N, float y[], float res[]);
+void dct2(const int M, const int N, float y[M][N], float res[M][N]);
+void idct(const int N, float a[N], float res[N]);
+void idct2(int M, int N, float y[M][N], float res[M][N]);
+
+#endif /* DCT2_H */
+
index a4bd8d53f6e39e1aa92f5c586b923a3088a22e49..45fa5fb72f1214a291e260fa1c683b3922a90cf3 100644 (file)
@@ -96,5 +96,9 @@ target_link_libraries(tsrc samplerate)
 add_executable(tlininterp tlininterp.c)
 add_executable(tdec tdec.c)
 
-add_executable(tc2wideband tc2wideband.c ../src/c2wideband.c ../src/octave.c)
+add_executable(tc2wideband tc2wideband.c t_helpers.c ../src/c2wideband.c ../src/octave.c)
 target_link_libraries(tc2wideband codec2)
+
+add_executable(tdct2 tdct2.c t_helpers.c ../src/dct2.c)
+target_link_libraries(tdct2 codec2)
+
diff --git a/codec2-dev/unittest/t_helpers.c b/codec2-dev/unittest/t_helpers.c
new file mode 100644 (file)
index 0000000..4f20743
--- /dev/null
@@ -0,0 +1,38 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "t_helpers.h"
+
+void test(char * tfn)
+{
+    fn = tfn;
+    printf("========================================\n");
+    printf("test function: %s\n", fn);
+    printf("========================================\n");
+}
+
+void test_failed()
+{
+    printf("Failed to calculate %s.\n", fn);
+    exit(1);
+}
+
+void test_failed_s(char * expected, char * res)
+{
+
+    printf("Failed to calculate %s.\n", fn);
+
+    printf("expected: %s\ngot: %s\n", expected, res);
+    exit(1);
+}
+
+void test_failed_f(float expected, float res)
+{
+
+    printf("Failed to calculate %s.\n", fn);
+    printf("expected: %f\ngot: %f\n", expected, res);
+    exit(1);
+}
+
diff --git a/codec2-dev/unittest/t_helpers.h b/codec2-dev/unittest/t_helpers.h
new file mode 100644 (file)
index 0000000..87665ee
--- /dev/null
@@ -0,0 +1,20 @@
+/* 
+ * File:   t_helpers.h
+ * Author: phil
+ *
+ * Created on 21 July 2017, 14:20
+ */
+
+#ifndef T_HELPERS_H
+#define        T_HELPERS_H
+
+void test(char * tfn);
+void test_failed();
+void test_failed_s(char * expected, char * res);
+void test_failed_f(float expected, float res);
+
+char *fn;
+
+
+#endif /* T_HELPERS_H */
+
index 7b41cd039f3d1bf160466fb489c7c716ea4a30f9..4ec01e74b7d6bc96a257e620c98e722d4d9fb3a1 100644 (file)
@@ -5,6 +5,7 @@
 #include <math.h>
 #include <unistd.h>
 
+#include "t_helpers.h"
 #include "c2wideband.h"
 #include "codec2.h"
 #include "defines.h"
@@ -29,219 +30,193 @@ void test_wideband_enc_dec();
 void setup_map(WIDEBAND_MAP * wb_map, int Nt, int K);
 void test_with_real_data(int argc, char *argv[]);
 
-char *fn;
+int main(int argc, char *argv[])
+{
 
-void test(char * tfn) {
-  fn = tfn;
-  printf("test function: %s\n", fn);
-}
-
-void test_failed() {
-  printf("Failed to calculate %s.\n", fn);
-  exit(1);
-}
+    printf("Testing file c2wideband\n");
 
-void test_failed_s(char * expected, char * res) {
+    test("mean");
+    int n = 5;
+    float in[] = {1.0, 2.0, 4.0, 6.0, 10.0};
+    float res_f;
+    float expect_f;
 
-  printf("Failed to calculate %s.\n", fn);
-
-  printf("expected: %s\ngot: %s\n", expected, res);
-  exit(1);
-}
+    res_f = mean(in, n);
+    expect_f = 23.0 / 5;
+    if (res_f != expect_f) {
+        test_failed_f(expect_f, res_f);
+    }
 
-void test_failed_f(float expected, float res) {
+    test("std");
+    res_f = std(in, n);
+    expect_f = 3.2000000000000002;
+    if (res_f != expect_f) {
+        test_failed_f(expect_f, res_f);
+    }
 
-  printf("Failed to calculate %s.\n", fn);
-  printf("expected: %f\ngot: %f\n", expected, res);
-  exit(1);
-}
 
-int main(int argc, char *argv[]) {
+    test("diff_de");
+    int rows = 5;
+    int cols = 3;
+    float D[rows][cols];
+    float E[rows][cols];
+    float res_diff_de[rows][cols];
+    float expect_diff_de[rows][cols];
+    int r, c;
+
+    for (r = 0; r < rows; r++) {
+        for (c = 0; c < cols; c++) {
+            float d = rand();
+            float e = rand();
+            D[r][c] = d;
+            E[r][c] = e;
+            expect_diff_de[r][c] = d - e;
+        }
+    }
 
-  printf("Testing file c2wideband\n");
+    diff_de(rows, cols, D, E, res_diff_de);
 
-  test("mean");
-  int n = 5;
-  float in[] = {1.0, 2.0, 4.0, 6.0, 10.0};
-  float res_f;
-  float expect_f;
+    if (memcmp(res_diff_de, expect_diff_de, rows * cols * sizeof (float)) != 0) {
+        test_failed();
+    }
 
-  res_f = mean(in, n);
-  expect_f = 23.0 / 5;
-  if (res_f != expect_f) {
-    test_failed_f(expect_f, res_f);
-  }
 
-  test("std");
-  res_f = std(in, n);
-  expect_f = 3.2000000000000002;
-  if (res_f != expect_f) {
-    test_failed_f(expect_f, res_f);
-  }
+    test("array_col_to_row");
+    float data[rows][cols];
+    float res_data[rows];
+    float expect_data[cols][rows];
+    for (r = 0; r < rows; r++) {
+        for (c = 0; c < cols; c++) {
+            float d = rand();
+            data[r][c] = d;
+            expect_data[c][r] = d;
+        }
+    }
+    for (c = 0; c < cols; c++) {
 
+        array_col_to_row(rows, cols, data, c, res_data);
+        if (memcmp(res_data, expect_data[c], cols * sizeof (float)) != 0) {
+            test_failed();
+        }
+    }
 
-  test("diff_de");
-  int rows = 5;
-  int cols = 3;
-  float D[rows][cols];
-  float E[rows][cols];
-  float res_diff_de[rows][cols];
-  float expect_diff_de[rows][cols];
-  int r,c;
 
-  for (r = 0; r < rows; r++) {
-    for (c = 0; c < cols; c++) {
-      float d = rand();
-      float e = rand();
-      D[r][c] = d;
-      E[r][c] = e;
-      expect_diff_de[r][c] = d - e;
-    }
-  }
+    test("std_on_cols");
 
-  diff_de(rows, cols, D, E, res_diff_de);
+    cols = 4;
 
-  if (memcmp(res_diff_de, expect_diff_de, rows * cols * sizeof (float)) != 0) {
-    test_failed();
-  }
+    float data_std[5][4] = {
+        {1.0, 1.0, 4.0, 6.0},
+        {1.0, 2.0, 5.5, 6.0},
+        {1.0, 4.0, 6.0, 6.0},
+        {1.0, 6.0, 3.3, 6.0},
+        {1.0, 10.0, -0.2, 6.0}
+    };
 
+    float res_data_std[cols];
+    float expect_data_std[] = {0.0, 3.2000000000000002, 2.1903424389807178, 0.0};
 
-  test("array_col_to_row");
-  float data[rows][cols];
-  float res_data[rows];
-  float expect_data[cols][rows];
-  for (r = 0; r < rows; r++) {
-    for (c = 0; c < cols; c++) {
-      float d = rand();
-      data[r][c] = d;
-      expect_data[c][r] = d;
+    std_on_cols(rows, cols, data_std, res_data_std);
+    if (memcmp(res_data_std, expect_data_std, cols * sizeof (float)) != 0) {
+        test_failed(fn);
     }
-  }
-  for (c = 0; c < cols; c++) {
 
-    array_col_to_row(rows, cols, data, c, res_data);
-    if (memcmp(res_data, expect_data[c], cols * sizeof (float)) != 0) {
-      test_failed();
-    }
-  }
+    test("mean_std_diff_de");
 
+    cols = 4;
 
-  test("std_on_cols");
+    float data_std_d[5][4] = {
+        {2.0, 1.0, 5.0, 6.0},
+        {1.0, 2.0, 10.5, 6.0},
+        {4.0, 4.0, 12.0, 6.0},
+        {1.0, 6.0, 4.6, 6.0},
+        {6.5, 10.0, 0.0, 6.0}
+    };
 
-  cols = 4;
+    float data_std_e[5][4] = {
+        {1.0, 0.0, 1.0, 0.0},
+        {0.0, 0.0, 5.0, 0.0},
+        {3.0, 0.0, 6.0, 0.0},
+        {0.0, 0.0, 1.3, 0.0},
+        {5.5, 0.0, 0.2, 0.0}
+    };
 
-  float data_std[5][4] = {
-    {1.0, 1.0, 4.0, 6.0},
-    {1.0, 2.0, 5.5, 6.0},
-    {1.0, 4.0, 6.0, 6.0},
-    {1.0, 6.0, 3.3, 6.0},
-    {1.0, 10.0, -0.2, 6.0}
-  };
+    float expected_msd = (3.2000000000000002 + 2.1903424389807178) / cols;
 
-  float res_data_std[cols];
-  float expect_data_std[] = {0.0, 3.2000000000000002, 2.1903424389807178, 0.0};
+    float res_msd;
+    res_msd = mean_std_diff_de(rows, cols, data_std_d, data_std_e);
+    if (abs(res_msd - expected_msd) > 0.00001) {
+        test_failed_f(res_msd, expected_msd);
+    }
 
-  std_on_cols(rows, cols, data_std, res_data_std);
-  if (memcmp(res_data_std, expect_data_std, cols * sizeof (float)) != 0) {
-    test_failed(fn);
-  }
+    test_wideband_enc_dec();
 
-  test("mean_std_diff_de");
+    test_with_real_data(argc, argv);
 
-  cols = 4;
+    return 1;
+}
 
-  float data_std_d[5][4] = {
-    {2.0, 1.0, 5.0, 6.0},
-    {1.0, 2.0, 10.5, 6.0},
-    {4.0, 4.0, 12.0, 6.0},
-    {1.0, 6.0, 4.6, 6.0},
-    {6.5, 10.0, 0.0, 6.0}
-  };
+void test_wideband_enc_dec()
+{
 
-  float data_std_e[5][4] = {
-    {1.0, 0.0, 1.0, 0.0},
-    {0.0, 0.0, 5.0, 0.0},
-    {3.0, 0.0, 6.0, 0.0},
-    {0.0, 0.0, 1.3, 0.0},
-    {5.5, 0.0, 0.2, 0.0}
-  };
+    test("wideband_enc_dec");
 
-  float expected_msd = (3.2000000000000002 + 2.1903424389807178) / cols;
+    short *buf;
+    unsigned char *bits;
+    //int            nsam, nbit, i, r;
+    int nbit;
+    int nsam;
+    int n_block_frames = C2WB_NT * C2WB_DEC;
+    int K = C2WB_K;
+    int Nt = C2WB_NT;
+    int Fs = C2WB_FS;
+    int dec = C2WB_DEC;
+    int Nblocks = n_block_frames;
+    int rate_K_surface_size = Nblocks * Nt * dec;
 
-  float res_msd;
-  res_msd = mean_std_diff_de(rows, cols, data_std_d, data_std_e);
-  if (abs(res_msd - expected_msd) > 0.00001) {
-    test_failed_f(res_msd, expected_msd);
-  }
+    struct CODEC2 * codec2 = codec2_create(CODEC2_MODE_WB);
 
-  test_wideband_enc_dec();
+    C2CONST c2const = c2const_create(Fs);
+    nsam = C2WB_SPERF; //codec2_samples_per_frame(codec2);
+    nbit = C2WB_BPERF; //codec2_bits_per_frame(codec2);
 
-  test_with_real_data(argc, argv);
+    buf = (short*) malloc(nsam * sizeof (short));
+    bits = (unsigned char*) malloc(nbit * sizeof (char));
 
-  return 1;
-}
+    WIDEBAND_MAP wb_map;
+    setup_map(&wb_map, Nt, K);
 
-void test_wideband_enc_dec() {
-
-  test("wideband_enc_dec");
-
-  short *buf;
-  unsigned char *bits;
-  //int            nsam, nbit, i, r;
-  int nbit;
-  int nsam;
-  int n_block_frames = C2WB_NT * C2WB_DEC;
-  int K = C2WB_K;
-  int Nt = C2WB_NT;
-  int Fs = C2WB_FS;
-  int dec = C2WB_DEC;
-  int Nblocks = n_block_frames;
-  int rate_K_surface_size = Nblocks * Nt * dec;
-
-  struct CODEC2 * codec2 = codec2_create(CODEC2_MODE_WB);
-
-  C2CONST c2const = c2const_create(Fs);
-  nsam = C2WB_SPERF; //codec2_samples_per_frame(codec2);
-  nbit = C2WB_BPERF; //codec2_bits_per_frame(codec2);
-
-  buf = (short*) malloc(nsam * sizeof (short));
-  bits = (unsigned char*) malloc(nbit * sizeof (char));
-
-  WIDEBAND_MAP wb_map;
-  setup_map(&wb_map, Nt, K);
-
-  MODEL model_block[n_block_frames];
-  MODEL model_block_[n_block_frames];
-  float dct2_sd[Nblocks];
-  int qn;
-  float rate_K_surface_block[rate_K_surface_size][K]; // rate K vecs for each frame, form a surface that makes pretty graphs   
-  float rate_K_surface_block_[rate_K_surface_size][K];
-  int n,i;
-
-  for (n = 0; n < n_block_frames; n++) {
-    model_block[n].L = 10;
-    model_block[n].Wo = (float) rand() / RAND_MAX;
-    for (i = 0; i < MAX_AMP + 1; i++) {
-      model_block[n].phi[i] = (float) rand() / RAND_MAX;
-      model_block[n].A[i] = (float) rand() / RAND_MAX;
+    MODEL model_block[n_block_frames];
+    MODEL model_block_[n_block_frames];
+    float dct2_sd[Nblocks];
+    int qn;
+    float rate_K_surface_block[rate_K_surface_size][K]; // rate K vecs for each frame, form a surface that makes pretty graphs   
+    float rate_K_surface_block_[rate_K_surface_size][K];
+    int n, i;
+
+    for (n = 0; n < n_block_frames; n++) {
+        model_block[n].L = 10;
+        model_block[n].Wo = (float) rand() / RAND_MAX;
+        for (i = 0; i < MAX_AMP + 1; i++) {
+            model_block[n].phi[i] = (float) rand() / RAND_MAX;
+            model_block[n].A[i] = (float) rand() / RAND_MAX;
+        }
+        model_block[n].voiced = round((float) rand() / RAND_MAX);
     }
-    model_block[n].voiced = round((float) rand() / RAND_MAX);
-  }
 
-  n = 0;
+    n = 0;
 
-  printf("setup complete. now calling the function\n");
+    printf("setup complete. now calling the function\n");
 
-  wideband_enc_dec(&c2const, n_block_frames, model_block, &wb_map,
-          model_block_, dct2_sd, &qn, &rate_K_surface_block[n], &rate_K_surface_block_[n]);
+    wideband_enc_dec(&c2const, n_block_frames, model_block, &wb_map,
+                     model_block_, dct2_sd, &qn, &rate_K_surface_block[n], &rate_K_surface_block_[n]);
 
-  codec2_destroy(codec2);
-  free(bits);
-  free(buf);
-  printf("made it to the end\n");
+    codec2_destroy(codec2);
+    free(bits);
+    free(buf);
+    printf("made it to the end\n");
 }
 
-
 /*---------------------------------------------------------------------------*\
 
   FUNCTION....: test_with_real_data()
@@ -257,53 +232,54 @@ void test_wideband_enc_dec() {
 
 \*---------------------------------------------------------------------------*/
 
-void test_with_real_data(int argc, char *argv[]) {
+void test_with_real_data(int argc, char *argv[])
+{
     int Nt = C2WB_NT;
     int Fs = C2WB_FS;
     int K = C2WB_K;
 
     C2CONST c2const = c2const_create(Fs);
-    int   n_samp = c2const.n_samp;
-    int   m_pitch = c2const.m_pitch;
-    short buf[n_samp];         /* input/output buffer                   */
-    float Sn[m_pitch];         /* float input speech samples            */
-    COMP  Sw[FFT_ENC];         /* DFT of Sn[]                           */
+    int n_samp = c2const.n_samp;
+    int m_pitch = c2const.m_pitch;
+    short buf[n_samp]; /* input/output buffer                   */
+    float Sn[m_pitch]; /* float input speech samples            */
+    COMP Sw[FFT_ENC]; /* DFT of Sn[]                           */
     codec2_fft_cfg fft_fwd_cfg; /* fwd FFT states                        */
-    float w[m_pitch];          /* time domain hamming window            */
-    COMP  W[FFT_ENC];          /* DFT of w[]                            */
+    float w[m_pitch]; /* time domain hamming window            */
+    COMP W[FFT_ENC]; /* DFT of w[]                            */
     MODEL model[FRAMES];
     void *nlp_states;
     float pitch, prev_f0;
-    int   i,m,f;
-    
+    int i, m, f;
+
     if (argc != 2) {
         printf("usage: ./tnewamp1 RawFile\n");
         exit(1);
     }
     nlp_states = nlp_create(&c2const);
-    prev_f0 = 1.0/P_MAX_S;
-    fft_fwd_cfg = codec2_fft_alloc(FFT_ENC, 0, NULL, NULL); 
-    make_analysis_window(&c2const,fft_fwd_cfg, w, W);
+    prev_f0 = 1.0 / P_MAX_S;
+    fft_fwd_cfg = codec2_fft_alloc(FFT_ENC, 0, NULL, NULL);
+    make_analysis_window(&c2const, fft_fwd_cfg, w, W);
 
-    for(i=0; i<m_pitch; i++) {
-       Sn[i] = 1.0;
+    for (i = 0; i < m_pitch; i++) {
+        Sn[i] = 1.0;
     }
 
     WIDEBAND_MAP wb_map;
     setup_map(&wb_map, Nt, K);
     int n_block_frames = C2WB_NT * C2WB_DEC;
 
-    float model_octave[FRAMES][MAX_AMP+2];    // model params in matrix format, useful for C <-> Octave  
-    float rate_K_surface[FRAMES][K];          // rate K vecs for each frame, form a surface that makes pretty graphs
+    float model_octave[FRAMES][MAX_AMP + 2]; // model params in matrix format, useful for C <-> Octave  
+    float rate_K_surface[FRAMES][K]; // rate K vecs for each frame, form a surface that makes pretty graphs
     float rate_K_surface_[FRAMES][K];
     MODEL model_block_[n_block_frames];
 
-    int Nblocks = FRAMES/n_block_frames;
+    int Nblocks = FRAMES / n_block_frames;
     float dct2_sd[Nblocks];
     int qn;
 
-    for(f=0; f<FRAMES; f++) {
-        for(m=0; m<MAX_AMP+2; m++) {
+    for (f = 0; f < FRAMES; f++) {
+        for (m = 0; m < MAX_AMP + 2; m++) {
             model_octave[f][m] = 0.0;
         }
     }
@@ -314,49 +290,49 @@ void test_with_real_data(int argc, char *argv[]) {
         exit(1);
     }
 
-    for(f=0; f<FRAMES; f++) {
-        assert(fread(buf,sizeof(short),n_samp,fin) == n_samp);
+    for (f = 0; f < FRAMES; f++) {
+        assert(fread(buf, sizeof (short), n_samp, fin) == n_samp);
 
         /* shift buffer of input samples, and insert new samples */
 
-       for(i=0; i<m_pitch-n_samp; i++) {
-           Sn[i] = Sn[i+n_samp];
-       }
-       for(i=0; i<n_samp; i++) {
-           Sn[i+m_pitch-n_samp] = buf[i];
+        for (i = 0; i < m_pitch - n_samp; i++) {
+            Sn[i] = Sn[i + n_samp];
+        }
+        for (i = 0; i < n_samp; i++) {
+            Sn[i + m_pitch - n_samp] = buf[i];
         }
 
-       /* Estimate Sinusoidal Model Parameters ----------------------*/
+        /* Estimate Sinusoidal Model Parameters ----------------------*/
 
-       nlp(nlp_states, Sn, n_samp, &pitch, Sw, W, &prev_f0);
-       model[f].Wo = TWO_PI/pitch;
+        nlp(nlp_states, Sn, n_samp, &pitch, Sw, W, &prev_f0);
+        model[f].Wo = TWO_PI / pitch;
 
-       dft_speech(&c2const, fft_fwd_cfg, Sw, Sn, w);
-       two_stage_pitch_refinement(&c2const, &model[f], Sw);
-       estimate_amplitudes(&model[f], Sw, W, 1);
+        dft_speech(&c2const, fft_fwd_cfg, Sw, Sn, w);
+        two_stage_pitch_refinement(&c2const, &model[f], Sw);
+        estimate_amplitudes(&model[f], Sw, W, 1);
         est_voicing_mbe(&c2const, &model[f], Sw, W);
 
-        fprintf(stderr,"f: %d Wo: %4.3f L: %d v: %d\n", f, model[f].Wo, model[f].L, model[f].voiced);
+        fprintf(stderr, "f: %d Wo: %4.3f L: %d v: %d\n", f, model[f].Wo, model[f].L, model[f].voiced);
 
         /* log some vectors for sinusoidal model */
+
         model_octave[f][0] = model[f].Wo;
         model_octave[f][1] = model[f].L;
 
-        for(m=1; m<=model[f].L; m++) {
-            model_octave[f][m+1] = model[f].A[m];
-        }        
+        for (m = 1; m <= model[f].L; m++) {
+            model_octave[f][m + 1] = model[f].A[m];
+        }
 
         /* once we have collected a block of samples process ----------*/
 
-        if (((f+1) % n_block_frames) == 0) {
+        if (((f + 1) % n_block_frames) == 0) {
 
             /* wideband processing ----------------------------------------*/
 
-            int block_st = f-n_block_frames+1;
+            int block_st = f - n_block_frames + 1;
             wideband_enc_dec(&c2const, n_block_frames, &model[block_st], &wb_map,
-                             model_block_, dct2_sd, &qn, 
-                             &rate_K_surface[block_st], 
+                             model_block_, dct2_sd, &qn,
+                             &rate_K_surface[block_st],
                              &rate_K_surface_[block_st]);
 
             fprintf(stderr, "  Performed a wideband_enc_dec() call, qn: %d\n", qn);
@@ -369,11 +345,11 @@ void test_with_real_data(int argc, char *argv[]) {
 
     /* save vectors in Octave format */
 
-    FILE *fout = fopen("tc2wideband_out.txt","wt");
+    FILE *fout = fopen("tc2wideband_out.txt", "wt");
     assert(fout != NULL);
     fprintf(fout, "# Created by tc2wideband.c\n");
-    octave_save_float(fout, "rate_K_surface_c", (float*)rate_K_surface, FRAMES, K, K);
-    octave_save_float(fout, "model_c", (float*)model_octave, FRAMES, MAX_AMP+2, MAX_AMP+2);
+    octave_save_float(fout, "rate_K_surface_c", (float*) rate_K_surface, FRAMES, K, K);
+    octave_save_float(fout, "model_c", (float*) model_octave, FRAMES, MAX_AMP + 2, MAX_AMP + 2);
     fclose(fout);
 
     printf("Done! Now run\n  octave:1> tc2wideband(\"../path/to/tc2wideband_out.txt\")\n");
diff --git a/codec2-dev/unittest/tdct2.c b/codec2-dev/unittest/tdct2.c
new file mode 100644 (file)
index 0000000..d369026
--- /dev/null
@@ -0,0 +1,99 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "t_helpers.h"
+#include "dct2.h"
+
+int main()
+{
+    int i;
+    int j;
+
+    test("dct");
+
+    int N = 4;
+
+    float y[] = {1, 2, 3, 4};
+    float expect_dct[] = {20, -6.30864406, 0, -0.44834153};
+    float res_dct[N];
+
+    dct(N, y, res_dct);
+
+
+    for (i = 0; i < N; i++) {
+        //printf("Result = %f : %f\n", res_dct[i].real, res_dct[i].imag);
+        //printf("Expected = %f\n", expect_dct[i]);
+        if (abs(expect_dct[i] - res_dct[i]) > 0.00001)
+            test_failed_f(expect_dct[i], res_dct[i]);
+    }
+
+    test("idct");
+    float *expect_idct = y;
+    float res_idct[N];
+
+    idct(N, res_dct, res_idct);
+
+
+    for (i = 0; i < N; i++) {
+        //printf("Result = %f : %f\n", res_dct[i].real, res_dct[i].imag);
+        //printf("Expected = %f\n", expect_dct[i]);
+        if (abs(expect_idct[i] - res_idct[i]) > 0.00001)
+            test_failed_f(expect_idct[i], res_idct[i]);
+    }
+
+    test("dct2");
+
+    int M = 3;
+    float expect_dct2[3][4] = {
+        { 180, -26.76530997, 8.48528137, 1.90215201},
+        { 3.46410162, -9.60123774, -7.34846923, -3.97696289},
+        { -66, 5.5432772, 4.24264069, 2.29610059}
+    };
+
+    float y2[3][4] = {
+        {1, 2, 3, 4},
+        {5, 6, 7, 8},
+        {3, 1, 2, 3}
+    };
+
+    
+    float res_dct2[M][N];
+    dct2(M, N, y2, res_dct2);
+
+
+
+    for (i = 0; i < M; i++) {
+        for (j = 0; j < N; j++) {
+            //printf("Result = %f : %f\n", res_dct[i].real, res_dct[i].imag);
+            //printf("Expected = %f\n", expect_dct[i]);
+            if (abs(expect_dct2[i][j] - res_dct2[i][j]) > 0.00001)
+                test_failed_f(expect_dct2[i][j], res_dct2[i][j]);
+        }
+    }
+    
+    test("idct2");
+
+
+    float expect_idct2[3][4] = {
+        {1, 2, 3, 4},
+        {5, 6, 7, 8},
+        {3, 1, 2, 3}
+    };
+
+    
+    float res_idct2[M][N];
+    idct2(M, N, res_dct2, res_idct2);
+
+    for (i = 0; i < M; i++) {
+        for (j = 0; j < N; j++) {
+            //printf("Result = %f : %f\n", res_dct[i].real, res_dct[i].imag);
+            //printf("Expected = %f\n", expect_dct[i]);
+            if (abs(expect_idct2[i][j] - res_idct2[i][j]) > 0.00001)
+                test_failed_f(expect_idct2[i][j], res_idct2[i][j]);
+        }
+    }
+}
\ No newline at end of file