Merge pull request #28 from suttonpd/master

Cleanup and bugfix of dft module
master
Ismael Gómez-Miguelez 11 years ago
commit 8f08380d46

@ -29,42 +29,42 @@
#ifndef DFT_H_
#define DFT_H_
#include <fftw3.h>
#include <stdbool.h>
#include "liblte/config.h"
/* dft is a frontend to the fftw3 library. It facilitates the computation of
* complex or real DFT, power spectral density, normalization, etc.
* It also supports the creation of multiple FFT plans for different FFT sizes
* or options, selecting a different one at runtime.
/* Generic DFT module.
* Supports one-dimensional complex and real transforms. Options are set
* using the dft_plan_set_x functions.
*
* Options (default is false):
* mirror - Rearranges negative and positive frequency bins. Swaps after
* transform for FORWARD, swaps before transform for BACKWARD.
* db - Provides output in dB (10*log10(x)).
* norm - Normalizes output (by sqrt(len) for complex, len for real).
* dc - Handles insertion and removal of null DC carrier internally.
*/
typedef enum {
COMPLEX_2_COMPLEX, REAL_2_REAL, COMPLEX_2_REAL
COMPLEX, REAL
}dft_mode_t;
typedef enum {
FORWARD, BACKWARD
}dft_dir_t;
#define DFT_MIRROR_PRE 1
#define DFT_PSD 2
#define DFT_OUT_DB 4
#define DFT_MIRROR_POS 8
#define DFT_NORMALIZE 16
#define DFT_DC_OFFSET 32
typedef struct LIBLTE_API {
int size;
int sign;
void *in;
void *out;
void *p;
int options;
dft_dir_t dir;
dft_mode_t mode;
int size; // DFT length
void *in; // Input buffer
void *out; // Output buffer
void *p; // DFT plan
bool forward; // Forward transform?
bool mirror; // Shift negative and positive frequencies?
bool db; // Provide output in dB?
bool norm; // Normalize output?
bool dc; // Handle insertion/removal of null DC carrier internally?
dft_dir_t dir; // Forward/Backward
dft_mode_t mode; // Complex/Real
}dft_plan_t;
typedef _Complex float dft_c_t;
@ -72,32 +72,24 @@ typedef float dft_r_t;
/* Create DFT plans */
LIBLTE_API int dft_plan(dft_plan_t *plan, const int dft_points,
dft_mode_t mode, dft_dir_t dir);
LIBLTE_API int dft_plan_c2c(dft_plan_t *plan, const int dft_points, dft_dir_t dir);
LIBLTE_API int dft_plan_r2r(dft_plan_t *plan, const int dft_points, dft_dir_t dir);
LIBLTE_API int dft_plan_c2r(dft_plan_t *plan, const int dft_points, dft_dir_t dir);
LIBLTE_API int dft_plan(dft_plan_t *plan, const int dft_points, dft_dir_t dir,
dft_mode_t type);
LIBLTE_API int dft_plan_c(dft_plan_t *plan, const int dft_points, dft_dir_t dir);
LIBLTE_API int dft_plan_r(dft_plan_t *plan, const int dft_points, dft_dir_t dir);
LIBLTE_API void dft_plan_free(dft_plan_t *plan);
/* Set options */
/* Create a vector of DFT plans */
LIBLTE_API int dft_plan_vector(dft_plan_t *plans, const int *dft_points,
dft_mode_t *modes, dft_dir_t *dirs, int nof_plans);
LIBLTE_API int dft_plan_multi_c2c(dft_plan_t *plans, const int *dft_points,
dft_dir_t dir, int nof_plans);
LIBLTE_API int dft_plan_multi_c2r(dft_plan_t *plans, const int *dft_points,
dft_dir_t dir, int nof_plans);
LIBLTE_API int dft_plan_multi_r2r(dft_plan_t *plans, const int *dft_points,
dft_dir_t dir, int nof_plans);
LIBLTE_API void dft_plan_free_vector(dft_plan_t *plans, int nof_plans);
LIBLTE_API void dft_plan_set_mirror(dft_plan_t *plan, bool val);
LIBLTE_API void dft_plan_set_db(dft_plan_t *plan, bool val);
LIBLTE_API void dft_plan_set_norm(dft_plan_t *plan, bool val);
LIBLTE_API void dft_plan_set_dc(dft_plan_t *plan, bool val);
/* Compute DFT */
LIBLTE_API void dft_run(dft_plan_t *plan, void *in, void *out);
LIBLTE_API void dft_run_c2c(dft_plan_t *plan, dft_c_t *in, dft_c_t *out);
LIBLTE_API void dft_run_r2r(dft_plan_t *plan, dft_r_t *in, dft_r_t *out);
LIBLTE_API void dft_run_c2r(dft_plan_t *plan, dft_c_t *in, dft_r_t *out);
LIBLTE_API void dft_run_c(dft_plan_t *plan, dft_c_t *in, dft_c_t *out);
LIBLTE_API void dft_run_r(dft_plan_t *plan, dft_r_t *in, dft_r_t *out);
#endif // DFT_H_

@ -42,7 +42,7 @@ int lte_fft_init_(lte_fft_t *q, lte_cp_t cp_type, int nof_prb, dft_dir_t dir) {
fprintf(stderr, "Error: Invalid nof_prb=%d\n", nof_prb);
return -1;
}
if (dft_plan_c2c(&q->fft_plan, symbol_sz, dir)) {
if (dft_plan_c(&q->fft_plan, symbol_sz, dir)) {
fprintf(stderr, "Error: Creating DFT plan\n");
return -1;
}
@ -52,12 +52,10 @@ int lte_fft_init_(lte_fft_t *q, lte_cp_t cp_type, int nof_prb, dft_dir_t dir) {
return -1;
}
q->fft_plan.options = DFT_NORMALIZE;
if (dir==FORWARD) {
q->fft_plan.options |= DFT_DC_OFFSET | DFT_MIRROR_POS;
} else {
q->fft_plan.options |= DFT_DC_OFFSET | DFT_MIRROR_PRE;
}
dft_plan_set_mirror(&q->fft_plan, true);
dft_plan_set_norm(&q->fft_plan, true);
dft_plan_set_dc(&q->fft_plan, true);
q->symbol_sz = symbol_sz;
q->nof_symbols = CP_NSYMB(cp_type);
q->cp_type = cp_type;
@ -109,7 +107,7 @@ void lte_fft_run(lte_fft_t *q, cf_t *input, cf_t *output) {
int i;
for (i=0;i<q->nof_symbols;i++) {
input += CP_ISNORM(q->cp_type)?CP_NORM(i, q->symbol_sz):CP_EXT(q->symbol_sz);
dft_run_c2c(&q->fft_plan, input, q->tmp);
dft_run_c(&q->fft_plan, input, q->tmp);
memcpy(output, &q->tmp[q->nof_guards], q->nof_re * sizeof(cf_t));
input += q->symbol_sz;
output += q->nof_re;
@ -124,7 +122,7 @@ void lte_ifft_run(lte_fft_t *q, cf_t *input, cf_t *output) {
for (i=0;i<q->nof_symbols;i++) {
cp_len = CP_ISNORM(q->cp_type)?CP_NORM(i, q->symbol_sz):CP_EXT(q->symbol_sz);
memcpy(&q->tmp[q->nof_guards], input, q->nof_re * sizeof(cf_t));
dft_run_c2c(&q->fft_plan, q->tmp, &output[cp_len]);
dft_run_c(&q->fft_plan, q->tmp, &output[cp_len]);
input += q->nof_re;
/* add CP */
memcpy(output, &output[q->symbol_sz], cp_len * sizeof(cf_t));

@ -136,6 +136,5 @@ int main(int argc, char **argv) {
n_prb++;
}
fftwf_cleanup();
exit(0);
}

@ -212,7 +212,6 @@ int main(int argc, char **argv) {
n = pbch_decode(&pbch, fft_buffer, ce, 1, &mib);
base_free();
fftwf_cleanup();
if (n < 0) {
fprintf(stderr, "Error decoding PBCH\n");

@ -233,7 +233,6 @@ int main(int argc, char **argv) {
printf("cfi: %d, distance: %d\n", cfi, distance);
base_free();
fftwf_cleanup();
if (n < 0) {
fprintf(stderr, "Error decoding PCFICH\n");

@ -304,6 +304,5 @@ int main(int argc, char **argv) {
} while (nof_frames <= max_frames);
base_free();
fftwf_cleanup();
exit(ret);
}

@ -329,6 +329,5 @@ int main(int argc, char **argv) {
goout:
base_free();
fftwf_cleanup();
exit(ret);
}

@ -272,7 +272,6 @@ int main(int argc, char **argv) {
}
base_free();
fftwf_cleanup();
if (n < 0) {
fprintf(stderr, "Error decoding phich\n");

@ -69,7 +69,7 @@ void sss_synch_m0m1(sss_synch_t *q, cf_t *input, int *m0, float *m0_value,
int i;
dft_run_c2c(&q->dftp_input, input, input_fft);
dft_run_c(&q->dftp_input, input, input_fft);
for (i = 0; i < N_SSS; i++) {
y[0][i] = input_fft[SSS_POS_SYMBOL + 2 * i];

@ -183,12 +183,13 @@ int pss_synch_set_N_id_2(pss_synch_t *q, int N_id_2) {
memset(q->pss_signal_freq, 0, PSS_LEN_FREQ * sizeof(cf_t));
memcpy(&pss_signal_pad[33], pss_signal_time, PSS_LEN * sizeof(cf_t));
if (dft_plan(&plan, PSS_LEN_FREQ - 1, COMPLEX_2_COMPLEX, BACKWARD)) {
if (dft_plan(&plan, PSS_LEN_FREQ - 1, BACKWARD, COMPLEX)) {
return -1;
}
plan.options = DFT_MIRROR_PRE | DFT_DC_OFFSET;
dft_plan_set_mirror(&plan, true);
dft_plan_set_dc(&plan, true);
dft_run_c2c(&plan, pss_signal_pad, q->pss_signal_freq);
dft_run_c(&plan, pss_signal_pad, q->pss_signal_freq);
vec_sc_prod_cfc(q->pss_signal_freq, (float) 1 / (PSS_LEN_FREQ - 1),
pss_signal_pad, PSS_LEN_FREQ);

@ -29,6 +29,7 @@
#include <strings.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "liblte/phy/sync/sss.h"
@ -42,11 +43,12 @@ void generate_N_id_1_table(int table[30][30]);
int sss_synch_init(sss_synch_t *q) {
bzero(q, sizeof(sss_synch_t));
if (dft_plan(&q->dftp_input, SSS_DFT_LEN, COMPLEX_2_COMPLEX, FORWARD)) {
if (dft_plan(&q->dftp_input, SSS_DFT_LEN, FORWARD, COMPLEX)) {
return -1;
}
generate_N_id_1_table(q->N_id_1_table);
q->dftp_input.options = DFT_MIRROR_POS | DFT_DC_OFFSET;
dft_plan_set_mirror(&q->dftp_input, true);
dft_plan_set_dc(&q->dftp_input, true);
return 0;
}

@ -158,8 +158,6 @@ int main(int argc, char **argv) {
sync_free(&sync);
lte_ifft_free(&ifft);
fftwf_cleanup();
printf("Ok\n");
exit(0);
}

@ -44,13 +44,13 @@ int conv_fft_cc_init(conv_fft_cc_t *state, int input_len, int filter_len) {
if (!state->input_fft || !state->filter_fft || !state->output_fft) {
return -1;
}
if (dft_plan(&state->input_plan,state->output_len,COMPLEX_2_COMPLEX,FORWARD)) {
if (dft_plan(&state->input_plan,state->output_len,FORWARD,COMPLEX)) {
return -2;
}
if (dft_plan(&state->filter_plan,state->output_len,COMPLEX_2_COMPLEX,FORWARD)) {
if (dft_plan(&state->filter_plan,state->output_len,FORWARD,COMPLEX)) {
return -3;
}
if (dft_plan(&state->output_plan,state->output_len,COMPLEX_2_COMPLEX,BACKWARD)) {
if (dft_plan(&state->output_plan,state->output_len,BACKWARD,COMPLEX)) {
return -4;
}
return 0;
@ -73,12 +73,12 @@ void conv_fft_cc_free(conv_fft_cc_t *state) {
int conv_fft_cc_run(conv_fft_cc_t *state, _Complex float *input, _Complex float *filter, _Complex float *output) {
dft_run_c2c(&state->input_plan, input, state->input_fft);
dft_run_c2c(&state->filter_plan, filter, state->filter_fft);
dft_run_c(&state->input_plan, input, state->input_fft);
dft_run_c(&state->filter_plan, filter, state->filter_fft);
vec_prod_ccc(state->input_fft,state->filter_fft,state->output_fft,state->output_len);
dft_run_c2c(&state->output_plan, state->output_fft, output);
dft_run_c(&state->output_plan, state->output_fft, output);
return state->output_len;

@ -33,73 +33,15 @@
#include "liblte/phy/utils/dft.h"
#define div(a,b) ((a-1)/b+1)
#define dft_ceil(a,b) ((a-1)/b+1)
#define dft_floor(a,b) (a/b)
int dft_plan_vector(dft_plan_t *plans, const int *dft_points,
dft_mode_t *modes, dft_dir_t *dirs, int nof_plans) {
int i;
for (i=0;i<nof_plans;i++) {
if (dft_plan(&plans[i], dft_points[i],modes[i],dirs[i])) {
return -1;
}
}
return 0;
}
int dft_plan_multi_c2c(dft_plan_t *plans, const int *dft_points,
dft_dir_t dir, int nof_plans) {
int i;
for (i=0;i<nof_plans;i++) {
if (dft_plan(&plans[i],dft_points[i],COMPLEX_2_COMPLEX,dir)) {
return -1;
}
}
return 0;
}
int dft_plan_multi_c2r(dft_plan_t *plans, const int *dft_points,
dft_dir_t dir, int nof_plans) {
int i;
for (i=0;i<nof_plans;i++) {
if (dft_plan(&plans[i], dft_points[i],COMPLEX_2_REAL,dir)) {
return -1;
}
}
return 0;
}
int dft_plan_multi_r2r(dft_plan_t *plans, const int *dft_points,
dft_dir_t dir, int nof_plans) {
int i;
for (i=0;i<nof_plans;i++) {
if (dft_plan(&plans[i], dft_points[i],REAL_2_REAL,dir)) {
return -1;
}
}
return 0;
}
int dft_plan(dft_plan_t *plan, const int dft_points,
dft_mode_t mode, dft_dir_t dir) {
switch(mode) {
case COMPLEX_2_COMPLEX:
if (dft_plan_c2c(plan,dft_points,dir)) {
return -1;
}
break;
case REAL_2_REAL:
if (dft_plan_r2r(plan,dft_points,dir)) {
return -1;
}
break;
case COMPLEX_2_REAL:
if (dft_plan_c2r(plan,dft_points,dir)) {
return -1;
}
break;
int dft_plan(dft_plan_t *plan, const int dft_points, dft_dir_t dir,
dft_mode_t mode) {
if(mode == COMPLEX){
return dft_plan_c(plan,dft_points,dir);
} else {
return dft_plan_r(plan,dft_points,dir);
}
return 0;
}
@ -109,166 +51,134 @@ static void allocate(dft_plan_t *plan, int size_in, int size_out, int len) {
plan->out = fftwf_malloc(size_out*len);
}
int dft_plan_c2c(dft_plan_t *plan, const int dft_points, dft_dir_t dir) {
int sign;
sign = (dir == FORWARD) ? FFTW_FORWARD : FFTW_BACKWARD;
int dft_plan_c(dft_plan_t *plan, const int dft_points, dft_dir_t dir) {
allocate(plan,sizeof(fftwf_complex),sizeof(fftwf_complex), dft_points);
int sign = (dir == FORWARD) ? FFTW_FORWARD : FFTW_BACKWARD;
plan->p = fftwf_plan_dft_1d(dft_points, plan->in, plan->out, sign, 0U);
if (!plan->p) {
return -1;
}
plan->size = dft_points;
plan->mode = COMPLEX_2_COMPLEX;
plan->mode = COMPLEX;
plan->dir = dir;
plan->forward = (dir==FORWARD)?true:false;
plan->mirror = false;
plan->db = false;
plan->norm = false;
return 0;
}
int dft_plan_r2r(dft_plan_t *plan, const int dft_points, dft_dir_t dir) {
int sign;
sign = (dir == FORWARD) ? FFTW_R2HC : FFTW_HC2R;
int dft_plan_r(dft_plan_t *plan, const int dft_points, dft_dir_t dir) {
allocate(plan,sizeof(float),sizeof(float), dft_points);
int sign = (dir == FORWARD) ? FFTW_R2HC : FFTW_HC2R;
plan->p = fftwf_plan_r2r_1d(dft_points, plan->in, plan->out, sign, 0U);
if (!plan->p) {
return -1;
}
plan->size = dft_points;
plan->mode = REAL_2_REAL;
plan->mode = REAL;
plan->dir = dir;
plan->forward = (dir==FORWARD)?true:false;
plan->mirror = false;
plan->db = false;
plan->norm = false;
return 0;
}
int dft_plan_c2r(dft_plan_t *plan, const int dft_points, dft_dir_t dir) {
if (dft_plan_c2c(plan, dft_points, dir)) {
return -1;
}
plan->mode = COMPLEX_2_REAL;
return 0;
void dft_plan_set_mirror(dft_plan_t *plan, bool val){
plan->mirror = val;
}
void dft_plan_set_db(dft_plan_t *plan, bool val){
plan->db = val;
}
void dft_plan_set_norm(dft_plan_t *plan, bool val){
plan->norm = val;
}
void dft_plan_set_dc(dft_plan_t *plan, bool val){
plan->dc = val;
}
static void copy(char *dst, char *src, int size_d, int len, int mirror, int dc_offset) {
int offset=dc_offset?1:0;
int hlen;
if (mirror == DFT_MIRROR_PRE) {
hlen = div(len,2);
static void copy_pre(char *dst, char *src, int size_d, int len,
bool forward, bool mirror, bool dc) {
int offset = dc?1:0;
if(mirror && !forward){
int hlen = dft_floor(len,2);
memset(dst,0,size_d*offset);
memcpy(&dst[offset*size_d], &src[size_d*hlen], size_d*(hlen-offset));
memcpy(&dst[hlen*size_d], src, size_d*(len - hlen));
} else if (mirror == DFT_MIRROR_POS) {
hlen = div(len,2);
memcpy(dst, &src[size_d*hlen], size_d*hlen);
memcpy(&dst[hlen*size_d], &src[size_d*offset], size_d*(len - hlen));
memcpy(&dst[size_d*offset], &src[size_d*hlen], size_d*(len-hlen-offset));
memcpy(&dst[(len-hlen)*size_d], src, size_d*hlen);
} else {
memcpy(dst,src,size_d*len);
}
}
static void copy_post(char *dst, char *src, int size_d, int len,
bool forward, bool mirror, bool dc) {
int offset = dc?1:0;
if(mirror && forward){
int hlen = dft_ceil(len,2);
memcpy(dst, &src[size_d*hlen], size_d*(len-hlen));
memcpy(&dst[(len-hlen)*size_d], &src[size_d*offset], size_d*(hlen-offset));
} else {
memcpy(dst,src,size_d*len);
}
}
void dft_run(dft_plan_t *plan, void *in, void *out) {
switch(plan->mode) {
case COMPLEX_2_COMPLEX:
dft_run_c2c(plan,in,out);
break;
case REAL_2_REAL:
dft_run_r2r(plan,in,out);
break;
case COMPLEX_2_REAL:
dft_run_c2r(plan,in,out);
break;
if(plan->mode == COMPLEX) {
dft_run_c(plan,in,out);
} else {
dft_run_r(plan,in,out);
}
}
void dft_run_c2c(dft_plan_t *plan, dft_c_t *in, dft_c_t *out) {
void dft_run_c(dft_plan_t *plan, dft_c_t *in, dft_c_t *out) {
float norm;
int i;
fftwf_complex *f_out = plan->out;
copy((char*) plan->in,(char*) in,sizeof(dft_c_t),plan->size,plan->options & DFT_MIRROR_PRE,
plan->options & DFT_DC_OFFSET);
copy_pre((char*)plan->in, (char*)in, sizeof(dft_c_t), plan->size,
plan->forward, plan->mirror, plan->dc);
fftwf_execute(plan->p);
if (plan->options & DFT_NORMALIZE) {
if (plan->norm) {
/**FIXME: Use VOLK */
norm = sqrtf(plan->size);
for (i=0;i<plan->size;i++) {
f_out[i] /= norm;
}
}
if (plan->options & DFT_OUT_DB) {
if (plan->db) {
for (i=0;i<plan->size;i++) {
f_out[i] = 10*log10(f_out[i]);
}
}
copy((char*) out,(char*) plan->out,sizeof(dft_c_t),plan->size,plan->options & DFT_MIRROR_POS,
plan->options & DFT_DC_OFFSET);
copy_post((char*)out, (char*)plan->out, sizeof(dft_c_t), plan->size,
plan->forward, plan->mirror, plan->dc);
}
void dft_run_r2r(dft_plan_t *plan, dft_r_t *in, dft_r_t *out) {
void dft_run_r(dft_plan_t *plan, dft_r_t *in, dft_r_t *out) {
float norm;
int i;
int len = plan->size;
float *f_out = plan->out;
copy((char*) plan->in,(char*) in,sizeof(dft_r_t),plan->size,plan->options & DFT_MIRROR_PRE,
plan->options & DFT_DC_OFFSET);
memcpy(plan->in,in,sizeof(dft_r_t)*plan->size);
fftwf_execute(plan->p);
if (plan->options & DFT_NORMALIZE) {
if (plan->norm) {
norm = plan->size;
for (i=0;i<len;i++) {
f_out[i] /= norm;
}
}
if (plan->options & DFT_PSD) {
for (i=0;i<(len+1)/2-1;i++) {
out[i] = sqrtf(f_out[i]*f_out[i]+f_out[len-i-1]*f_out[len-i-1]);
}
}
if (plan->options & DFT_OUT_DB) {
if (plan->db) {
for (i=0;i<len;i++) {
out[i] = 10*log10(out[i]);
}
}
}
void dft_run_c2r(dft_plan_t *plan, dft_c_t *in, dft_r_t *out) {
int i;
float norm;
float *f_out = plan->out;
copy((char*) plan->in,(char*) in,sizeof(dft_r_t),plan->size,plan->options & DFT_MIRROR_PRE,
plan->options & DFT_DC_OFFSET);
fftwf_execute(plan->p);
if (plan->options & DFT_NORMALIZE) {
norm = plan->size;
for (i=0;i<plan->size;i++) {
f_out[i] /= norm;
}
}
if (plan->options & DFT_PSD) {
for (i=0;i<plan->size;i++) {
out[i] = (__real__ f_out[i])*(__real__ f_out[i])+
(__imag__ f_out[i])*(__imag__ f_out[i]);
if (!(plan->options & DFT_OUT_DB)) {
out[i] = sqrtf(out[i]);
}
}
}
if (plan->options & DFT_OUT_DB) {
for (i=0;i<plan->size;i++) {
out[i] = 10*log10(out[i]);
f_out[i] = 10*log10(f_out[i]);
}
}
memcpy(out,plan->out,sizeof(dft_r_t)*plan->size);
}
void dft_plan_free(dft_plan_t *plan) {
if (!plan) return;
if (!plan->size) return;
@ -278,12 +188,5 @@ void dft_plan_free(dft_plan_t *plan) {
bzero(plan, sizeof(dft_plan_t));
}
void dft_plan_free_vector(dft_plan_t *plan, int nof_plans) {
int i;
for (i=0;i<nof_plans;i++) {
dft_plan_free(&plan[i]);
}
}

@ -0,0 +1,36 @@
#
# Copyright 2012-2013 The libLTE Developers. See the
# COPYRIGHT file at the top-level directory of this distribution.
#
# This file is part of the libLTE library.
#
# libLTE is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# libLTE 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 Lesser General Public License for more details.
#
# A copy of the GNU Lesser General Public License can be found in
# the LICENSE file in the top-level directory of this distribution
# and at http://www.gnu.org/licenses/.
#
########################################################################
# DFT TEST
########################################################################
ADD_EXECUTABLE(dft_test dft_test.c)
TARGET_LINK_LIBRARIES(dft_test lte_phy)
ADD_TEST(dft_test dft_test)
ADD_TEST(dft_reverse dft_test -b) # Backwards first
ADD_TEST(dft_mirror dft_test -m) # Mirror the freq bins
ADD_TEST(dft_norm dft_test -n) # Normalize each transform
ADD_TEST(dft_dc dft_test -b -d) # Backwards first & handle dc internally
ADD_TEST(dft_odd dft_test -N 255) # Odd-length
ADD_TEST(dft_odd_dc dft_test -N 255 -b -d) # Odd-length, backwards first, handle dc

@ -0,0 +1,137 @@
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <complex.h>
#include "liblte/phy/utils/dft.h"
typedef _Complex float cf_t;
int N = 256;
bool forward = true;
bool mirror = false;
bool norm = false;
bool dc = false;
void usage(char *prog) {
printf("Usage: %s\n", prog);
printf("\t-N Transform size [Default 256]\n");
printf("\t-b Backwards transform first [Default Forwards]\n");
printf("\t-m Mirror the transform freq bins [Default false]\n");
printf("\t-n Normalize the transform output [Default false]\n");
printf("\t-d Handle insertion/removal of null DC carrier internally [Default false]\n");
}
void parse_args(int argc, char **argv) {
int opt;
while ((opt = getopt(argc, argv, "Nbmnd")) != -1) {
switch (opt) {
case 'N':
N = atoi(argv[optind]);
break;
case 'b':
forward = false;
break;
case 'm':
mirror = true;
break;
case 'n':
norm = true;
break;
case 'd':
dc = true;
break;
default:
usage(argv[0]);
exit(-1);
}
}
}
void print(cf_t* in, int len)
{
for(int i=0;i<len;i++)
{
float re = crealf(in[i]);
float im = cimagf(in[i]);
printf("%f+%fi, ", re, im);
}
printf("\n\n");
}
int test_dft(cf_t* in){
int res = 0;
dft_plan_t plan;
if(forward){
dft_plan(&plan, N, FORWARD, COMPLEX);
} else {
dft_plan(&plan, N, BACKWARD, COMPLEX);
}
dft_plan_set_mirror(&plan, mirror);
dft_plan_set_norm(&plan, norm);
dft_plan_set_dc(&plan, dc);
cf_t* out1 = malloc(sizeof(cf_t)*N);
cf_t* out2 = malloc(sizeof(cf_t)*N);
bzero(out1, sizeof(cf_t)*N);
bzero(out2, sizeof(cf_t)*N);
print(in, N);
dft_run(&plan, in, out1);
print(out1, N);
dft_plan_t plan_rev;
if(!forward){
dft_plan(&plan_rev, N, FORWARD, COMPLEX);
} else {
dft_plan(&plan_rev, N, BACKWARD, COMPLEX);
}
dft_plan_set_mirror(&plan_rev, mirror);
dft_plan_set_norm(&plan_rev, norm);
dft_plan_set_dc(&plan_rev, dc);
dft_run(&plan_rev, out1, out2);
print(out2, N);
if(!norm){
cf_t n = N+0*I;
for(int i=0;i<N;i++)
out2[i] /= n;
}
for(int i=0;i<N;i++){
float diff = cabsf(in[i] - out2[i]);
if(diff > 0.01)
res = -1;
}
dft_plan_free(&plan);
dft_plan_free(&plan_rev);
free(out1);
free(out2);
return res;
}
int main(int argc, char **argv) {
parse_args(argc, argv);
cf_t* in = malloc(sizeof(cf_t)*N);
bzero(in, sizeof(cf_t)*N);
for(int i=1;i<N-1;i++)
{
float re = 100*(float)rand()/RAND_MAX;
float im = 100*(float)rand()/RAND_MAX;
in[i] = re + im*I;
}
if(test_dft(in) != 0)
return -1;
free(in);
printf("Done\n");
exit(0);
}
Loading…
Cancel
Save