From: philayres Date: Tue, 25 Jul 2017 15:48:13 +0000 (+0000) Subject: Added dct2.c with 1-D and 2-D DCT and inverse DCT functions, based off of the existin... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=5befa61a062e653737fb807a686aa18c37340637;p=freetel-svn-tracking.git Added dct2.c with 1-D and 2-D DCT and inverse DCT functions, based off of the existing real number FFT functions. 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 --- diff --git a/codec2-dev/src/CMakeLists.txt b/codec2-dev/src/CMakeLists.txt index 7d758d37..66a6fe15 100644 --- a/codec2-dev/src/CMakeLists.txt +++ b/codec2-dev/src/CMakeLists.txt @@ -230,6 +230,7 @@ set(CODEC2_SRCS varicode.c modem_stats.c c2wideband.c + dct2.c ) set(CODEC2_PUBLIC_HEADERS diff --git a/codec2-dev/src/c2wideband.c b/codec2-dev/src/c2wideband.c index c41a7158..684b7e48 100644 --- a/codec2-dev/src/c2wideband.c +++ b/codec2-dev/src/c2wideband.c @@ -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 index 00000000..641bba44 --- /dev/null +++ b/codec2-dev/src/dct2.c @@ -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 +#include +#include +#include + +#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 index 00000000..87665ee8 --- /dev/null +++ b/codec2-dev/unittest/t_helpers.h @@ -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 */ + diff --git a/codec2-dev/unittest/tc2wideband.c b/codec2-dev/unittest/tc2wideband.c index 7b41cd03..4ec01e74 100644 --- a/codec2-dev/unittest/tc2wideband.c +++ b/codec2-dev/unittest/tc2wideband.c @@ -5,6 +5,7 @@ #include #include +#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 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 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 index 00000000..d3690269 --- /dev/null +++ b/codec2-dev/unittest/tdct2.c @@ -0,0 +1,99 @@ +#include +#include +#include +#include +#include +#include + +#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