int rows = Nt;
int cols = K;
int dec = C2WB_DEC;
+
+ // one time configuration of DCT & IDCT
+
+ codec2_dct_cfg dct_cfg_n = dct_config(cols);
+ codec2_dct_cfg dct_cfg_m = dct_config(rows);
+ codec2_dct_cfg idct_cfg_n = idct_config(cols);
+ codec2_dct_cfg idct_cfg_m = idct_config(rows);
+
//printf("starting iteration\n");
// iterate through the frames in the block
int f;
float D[rows][cols];
float E[rows][cols];
//printf("dct2\n");
- dct2(rows, cols, &rate_K_surface_block[f], D);
+ dct2(dct_cfg_m, dct_cfg_n, 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
//TODO ???
//= [sqrt(dec)*E; zeros(Nt*(dec-1), K)];
//printf("idct2\n");
- idct2(rows, cols, inrks, &rate_K_surface_block_[f]);
+ idct2(idct_cfg_m, idct_cfg_n, 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");
+/*---------------------------------------------------------------------------*\
+
+ FILE........: dct2.c
+ AUTHOR......: Phil Ayres
+ DATE CREATED: July 2017
+
+ * DCT functions based on existing Codec 2 FFT
+ *
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright David Rowe 2017
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License version 2.1, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+ */
+
#include "dct2.h"
+
+/*
+ * Call dct_config with the number of parameters to be supplied to a subsequent
+ * DCT or DCT2 function. Returns a configuration that must be
+ * passed as the initial argument to the subsequent calls.
+ */
+codec2_dct_cfg dct_config(int P)
+{
+ return codec2_fftr_alloc(P * 2, 0, NULL, NULL);
+}
+
+/*
+ * Call idct_config with the number of parameters to be supplied to a subsequent
+ * inverse DCT or inverse DCT2 function. Returns a configuration that must be
+ * passed as the initial argument to the subsequent calls.
+ */
+codec2_dct_cfg idct_config(int P)
+{
+ return codec2_fftr_alloc(P * 2, 1, NULL, NULL);
+}
+
+/*
+ * Free memory from previous dct_config() and idct_config() calls
+ */
+void dct_cfg_free(codec2_dct_cfg cfg)
+{
+ codec2_fft_free((codec2_fft_cfg)cfg);
+}
+
+
/*
* 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])
+void dct(codec2_dct_cfg cfg, const int N, float y[N], float res[N])
{
int i;
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];
* Built on the 1-D DCT.
*/
-void dct2(const int M, const int N, float y[M][N], float res[M][N])
+void dct2(codec2_dct_cfg cfg_m, codec2_dct_cfg cfg_n, 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]);
+ dct(cfg_n, N, y[i], a[i]);
}
- float a_col[M];
- float b_col[M];
+
+ float a_tran[N][M];
- for (j = 0; j < N; j++) {
- for (i = 0; i < M; i++) {
- a_col[i] = a[i][j];
+ float b_tran[N][M];
+// for (j = 0; j < N; j++) {
+// for (i = 0; i < M; i++) {
+// a_col[i] = a[i][j];
+// }
+//
+// dct(cfg, M, a_col, b_col);
+//
+// for (i = 0; i < M; i++) {
+// res[i][j] = b_col[i];
+// }
+// }
+
+// for i in range(M):
+ // a[i,:] = dct(y[i,:])
+
+ for(i=0; i<M; i++){
+ dct(cfg_n, N, y[i], a[i]);
+
+ for(j=0; j<N;j++){
+ a_tran[j][i] = a[i][j];
}
+
+// printf("y[i,:]= ");
+// for(j=0; j<N;j++){
+// printf("%f ", y[i][j]);
+// }
+// printf("\n");
+//
+// printf("a[i,:]= ");
+// for(j=0; j<N;j++){
+// printf("%f ", a[i][j]);
+// }
+// printf("\n");
+
+ }
+
+
+
+// for j in range(N):
+ // b[:,j] = dct(a[:,j])
+
- dct(M, a_col, b_col);
-
- for (i = 0; i < M; i++) {
- res[i][j] = b_col[i];
+ for(j=0; j<N;j++){
+
+ dct(cfg_m, M, a_tran[j], b_tran[j]);
+
+ for(i=0; i<M; i++){
+ res[i][j] = b_tran[j][i];
}
+
+// printf("a[:,j]= ");
+// for(i=0; i<M;i++){
+// printf("%f ", a_tran[j][i]);
+// }
+// printf("\n");
+//
+// printf("b[:,j]= ");
+// for(i=0; i<M;i++){
+// printf("%f ", b_tran[j][i]);
+// }
+// printf("\n");
+
}
+
+
}
-//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])
+void idct(codec2_dct_cfg cfg, 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 nfft = N * 2;
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);
+ //printf("phi: %f+%fj\n", phi[i].real, phi[i].imag);
ac.real = a[i];
ac.imag = 0;
c[i] = cmult(phi[i], ac);
+ //printf("c: %f+%fj\n", c[i].real, c[i].imag);
}
c[N].real = 0;
// Scale the result
for(i=0;i<N; i++){
res[i] = res_fft[i] / nfft;
+ //printf("res: %f\n", res[i]);
}
}
-//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])
+void idct2(codec2_dct_cfg cfg_m, codec2_dct_cfg cfg_n, 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]);
+ idct(cfg_n, N, y[i], a[i]);
}
float a_col[M];
a_col[i] = a[i][j];
}
- idct(M, a_col, b_col);
+ idct(cfg_m, M, a_col, b_col);
for (i = 0; i < M; i++) {
res[i][j] = b_col[i];
-/*
- * File: dct2.h
- * Author: Phil Ayres
- *
- * Created on 21 July 2017, 14:24
+/*---------------------------------------------------------------------------*\
+
+ FILE........: dct2.h
+ AUTHOR......: Phil Ayres
+ DATE CREATED: July 2017
+
+ * DCT functions based on existing Codec 2 FFT
+ *
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright David Rowe 2017
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License version 2.1, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>.
+
*/
#ifndef DCT2_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]);
+typedef codec2_fftr_cfg codec2_dct_cfg;
+
+void dct(codec2_dct_cfg cfg, const int N, float y[], float res[]);
+void dct2(codec2_dct_cfg cfg_m, codec2_dct_cfg cfg_n, const int M, const int N, float y[M][N], float res[M][N]);
+void idct(codec2_dct_cfg cfg, const int N, float a[N], float res[N]);
+void idct2(codec2_dct_cfg cfg_m, codec2_dct_cfg cfg_n, int M, int N, float y[M][N], float res[M][N]);
+codec2_dct_cfg dct_config(int P);
+codec2_dct_cfg idct_config(int P);
+void dct_cfg_free(codec2_dct_cfg cfg);
#endif /* DCT2_H */
-/*
- * File: t_helpers.h
- * Author: phil
- *
- * Created on 21 July 2017, 14:20
+/*---------------------------------------------------------------------------*\
+
+ FILE........: t_helpers.c
+ AUTHOR......: Phil Ayres
+ DATE CREATED: July 2017
+
+ * Simple helper functions for unit tests
+ *
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright David Rowe 2017
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License version 2.1, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>.
+
*/
#ifndef T_HELPERS_H
+/*---------------------------------------------------------------------------*\
+
+ FILE........: tdct2.c
+ AUTHOR......: Phil Ayres
+ DATE CREATED: July 2017
+
+ * Unit test for DCT & DCT2 functions
+ *
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright David Rowe 2017
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License version 2.1, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+ */
+
+
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
int N = 4;
+ codec2_dct_cfg dct_cfg = dct_config(N);
+
float y[] = {1, 2, 3, 4};
float expect_dct[] = {20, -6.30864406, 0, -0.44834153};
float res_dct[N];
- dct(N, y, res_dct);
+ dct(dct_cfg, N, y, res_dct);
for (i = 0; i < N; i++) {
test_failed_f(expect_dct[i], res_dct[i]);
}
+ dct_cfg_free(dct_cfg);
+
test("idct");
+
+ codec2_dct_cfg idct_cfg = idct_config(N);
+
float *expect_idct = y;
float res_idct[N];
- idct(N, res_dct, res_idct);
+ idct(idct_cfg, N, res_dct, res_idct);
for (i = 0; i < N; i++) {
test_failed_f(expect_idct[i], res_idct[i]);
}
+
+ dct_cfg_free(idct_cfg);
+
test("dct2");
int M = 3;
+
+ codec2_dct_cfg dct_cfg_n = dct_config(N);
+ codec2_dct_cfg dct_cfg_m = dct_config(M);
+
+
float expect_dct2[3][4] = {
{ 180, -26.76530997, 8.48528137, 1.90215201},
{ 3.46410162, -9.60123774, -7.34846923, -3.97696289},
{3, 1, 2, 3}
};
-
- float res_dct2[M][N];
- dct2(M, N, y2, res_dct2);
-
-
+ float res_dct2[M][N];
+ dct2(dct_cfg_m, dct_cfg_n, M, N, y2, res_dct2);
+
+ /*
+ printf("Result\n");
+ for (i = 0; i < M; i++) {
+ for (j = 0; j < N; j++) {
+ printf("%f ", res_dct2[i][j]);
+
+ //if (abs(expect_dct2[i][j] - res_dct2[i][j]) > 0.00001)
+ // test_failed_f(expect_dct2[i][j], res_dct2[i][j]);
+ }
+ printf("\n");
+ }
+ printf("Expected\n");
+ for (i = 0; i < M; i++) {
+ for (j = 0; j < N; j++) {
+ printf("%f ", expect_dct2[i][j]);
+
+ //if (abs(expect_dct2[i][j] - res_dct2[i][j]) > 0.00001)
+ // test_failed_f(expect_dct2[i][j], res_dct2[i][j]);
+ }
+ printf("\n");
+ }
+ */
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]);
}
+
}
-
+
+
+
+ dct_cfg_free(dct_cfg_m);
+ dct_cfg_free(dct_cfg_n);
test("idct2");
+ codec2_dct_cfg idct_cfg_n = idct_config(N);
+ codec2_dct_cfg idct_cfg_m = idct_config(M);
+
float expect_idct2[3][4] = {
{1, 2, 3, 4},
{3, 1, 2, 3}
};
-
+
float res_idct2[M][N];
- idct2(M, N, res_dct2, res_idct2);
+ idct2(idct_cfg_m, idct_cfg_n, 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]);
+ //printf("Result = %f \n", res_idct2[i][j]);
+ //printf("Expected = %f\n", expect_idct2[i][j]);
if (abs(expect_idct2[i][j] - res_idct2[i][j]) > 0.00001)
test_failed_f(expect_idct2[i][j], res_idct2[i][j]);
}
}
+
+ dct_cfg_free(idct_cfg_m);
+ dct_cfg_free(idct_cfg_n);
+
+
+ printf("OK!\n");
}
\ No newline at end of file
projects / freetel-svn-tracking.git / commitdiff