From f4da40ffa7cbf1a075a537ce7b107802ae66631c Mon Sep 17 00:00:00 2001 From: ha7ilm Date: Sat, 18 Feb 2017 15:49:04 +0100 Subject: [PATCH] fractional_decimator_ff gets new algorithm based on Lagrange interpolator formula --- .gitignore | 1 + csdr.c | 48 +- grc_tests/test_fractional_decimator.grc | 1532 ++++++++++++----------- libcsdr.c | 94 +- libcsdr.h | 25 +- 5 files changed, 957 insertions(+), 743 deletions(-) diff --git a/.gitignore b/.gitignore index 1d4be5f..71766fd 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ ddcd *.so tags dumpvect.*.vect +grc_tests/top_block.py diff --git a/csdr.c b/csdr.c index 5506ab8..6035ac9 100644 --- a/csdr.c +++ b/csdr.c @@ -97,6 +97,7 @@ char usage[]= " agc_ff [hang_time [reference [attack_rate [decay_rate [max_gain [attack_wait [filter_alpha]]]]]]]\n" " fastagc_ff [block_size [reference]]\n" " rational_resampler_ff [transition_bw [window]]\n" +" old_fractional_decimator_ff [transition_bw [window]]\n" " fractional_decimator_ff [transition_bw [window]]\n" " fft_cc [window [--octave] [--benchmark]]\n" " logpower_cf [add_db]\n" @@ -1343,14 +1344,57 @@ int main(int argc, char *argv[]) firdes_lowpass_f(taps, taps_length, 0.59*0.5/(rate-transition_bw), window); //0.6 const to compensate rolloff //for(int=0;i initialized to zero + fprintf(stderr,"fractional_decimator_ff: not using taps\n"); + fractional_decimator_ff_t d = fractional_decimator_ff_init(rate, 4, NULL, 0); for(;;) { FEOF_CHECK; if(d.input_processed==0) d.input_processed=the_bufsize; else memcpy(input_buffer, input_buffer+d.input_processed, sizeof(float)*(the_bufsize-d.input_processed)); fread(input_buffer+(the_bufsize-d.input_processed), sizeof(float), d.input_processed, stdin); - d = fractional_decimator_ff(input_buffer, output_buffer, the_bufsize, rate, taps, taps_length, d); + fractional_decimator_ff(input_buffer, output_buffer, the_bufsize, &d); + fwrite(output_buffer, sizeof(float), d.output_size, stdout); + TRY_YIELD; + } + } + + if(!strcmp(argv[1],"old_fractional_decimator_ff")) + { + //Process the params + if(argc<=2) return badsyntax("need required parameters (rate)"); + float rate; + sscanf(argv[2],"%g",&rate); + + float transition_bw=0.03; + if(argc>=4) sscanf(argv[3],"%g",&transition_bw); + + window_t window = WINDOW_DEFAULT; + if(argc>=5) + { + window = firdes_get_window_from_string(argv[4]); + } + else fprintf(stderr,"old_fractional_decimator_ff: window = %s\n",firdes_get_string_from_window(window)); + + if(!initialize_buffers()) return -2; + sendbufsize(the_bufsize / rate); + + if(rate==1) clone_(the_bufsize); //copy input to output in this special case (and stick in this function). + + //Generate filter taps + int taps_length = firdes_filter_len(transition_bw); + fprintf(stderr,"old_fractional_decimator_ff: taps_length = %d\n",taps_length); + float* taps = (float*)malloc(sizeof(float)*taps_length); + firdes_lowpass_f(taps, taps_length, 0.59*0.5/(rate-transition_bw), window); //0.6 const to compensate rolloff + //for(int=0;i initialized to zero + for(;;) + { + FEOF_CHECK; + if(d.input_processed==0) d.input_processed=the_bufsize; + else memcpy(input_buffer, input_buffer+d.input_processed, sizeof(float)*(the_bufsize-d.input_processed)); + fread(input_buffer+(the_bufsize-d.input_processed), sizeof(float), d.input_processed, stdin); + d = old_fractional_decimator_ff(input_buffer, output_buffer, the_bufsize, rate, taps, taps_length, d); fwrite(output_buffer, sizeof(float), d.output_size, stdout); TRY_YIELD; } diff --git a/grc_tests/test_fractional_decimator.grc b/grc_tests/test_fractional_decimator.grc index 3b587fe..c57bd3c 100644 --- a/grc_tests/test_fractional_decimator.grc +++ b/grc_tests/test_fractional_decimator.grc @@ -1,49 +1,49 @@ - - + + Tue Nov 25 18:16:05 2014 options - - id - top_block - - - _enabled - True - - - title - - author - - description - - window_size 1280, 1024 - - generate_options - wx_gui - category Custom - run_options - prompt + comment + - run + description + + + + _enabled True + + _coordinate + (10, 10) + + + _rotation + 0 + + + generate_options + wx_gui + + + id + top_block + max_nouts 0 @@ -53,36 +53,32 @@ - alias + run_options + prompt + + + run + True + + + thread_safe_setters - _coordinate - (10, 10) - - - _rotation - 0 + title + variable - id - decimation + comment + _enabled True - - value - 3.2 - - - alias - - _coordinate (376, 11) @@ -91,52 +87,84 @@ _rotation 0 - - - variable id - samp_rate_2 + decimation + + + value + 1.09352 + + + + variable_slider + + comment + + + + converver + float_converter + + + value + 4e3 _enabled True - - value - samp_rate/decimation - - - alias - - _coordinate - (272, 11) + (16, 291) _rotation 0 + + grid_pos + 2, 1, 1, 1 + + + id + input_freq + + + label + + + + max + samp_rate/2 + + + min + 0 + + + notebook + + + + num_steps + 100 + + + style + wx.SL_HORIZONTAL + variable - id - samp_rate + comment + _enabled True - - value - 240000 - - - alias - - _coordinate (176, 11) @@ -145,16 +173,91 @@ _rotation 0 + + id + samp_rate + + + value + 240000 + + + + variable + + comment + + + + _enabled + True + + + _coordinate + (272, 11) + + + _rotation + 0 + + + id + samp_rate_2 + + + value + samp_rate/decimation + analog_sig_source_x + + amp + 1 + + + alias + + + + comment + + + + affinity + + + + _enabled + True + + + freq + input_freq + + + _coordinate + (8, 83) + + + _rotation + 0 + id analog_sig_source_x_0 - _enabled - True + maxoutbuf + 0 + + + minoutbuf + 0 + + + offset + 0 type @@ -168,309 +271,24 @@ waveform analog.GR_COS_WAVE - - freq - input_freq - - - amp - 1 - - - offset - 0 - - - alias - - - - affinity - - - - minoutbuf - 0 - - - maxoutbuf - 0 - - - _coordinate - (8, 83) - - - _rotation - 0 - - - - wxgui_fftsink2 - - id - wxgui_fftsink2_1 - - - _enabled - True - - - type - float - - - title - Original signal - - - samp_rate - samp_rate - - - baseband_freq - 0 - - - y_per_div - 10 - - - y_divs - 10 - - - ref_level - 0 - - - ref_scale - 2.0 - - - fft_size - 1024 - - - fft_rate - 15 - - - peak_hold - False - - - average - False - - - avg_alpha - 0 - - - win - None - - - win_size - - - - grid_pos - - - - notebook - nb1,1 - - - freqvar - None - - - alias - - - - affinity - - - - _coordinate - (360, 355) - - - _rotation - 0 - - - - ha5kfu_execproc_xx - - id - ha5kfu_execproc_xx_0 - - - _enabled - True - - - type - ff - - - commandline - "csdr fractional_decimator_ff "+str(decimation) - - - alias - - - - affinity - - - - minoutbuf - 0 - - - maxoutbuf - 0 - - - _coordinate - (352, 115) - - - _rotation - 0 - - - - wxgui_scopesink2 - - id - wxgui_scopesink2_0_0 - - - _enabled - True - - - type - float - - - title - Original signal - - - samp_rate - samp_rate - - - v_scale - 0 - - - v_offset - 0 - - - t_scale - 0 - - - ac_couple - False - - - xy_mode - False - - - num_inputs - 1 - - - win_size - - - - grid_pos - - - - notebook - nb1,0 - - - trig_mode - wxgui.TRIG_MODE_AUTO - - - y_axis_label - Counts - - - alias - - - - affinity - - - - _coordinate - (360, 235) - - - _rotation - 0 - blocks_throttle - - id - blocks_throttle_0 - - - _enabled - True - - - type - float - - - samples_per_second - samp_rate - - - vlen - 1 - - - ignoretag - True - alias + + comment + + affinity - minoutbuf - 0 - - - maxoutbuf - 0 + _enabled + True _coordinate @@ -480,202 +298,72 @@ _rotation 0 - - - wxgui_fftsink2 id - wxgui_fftsink2_0 + blocks_throttle_0 - _enabled + ignoretag True + + maxoutbuf + 0 + + + minoutbuf + 0 + + + samples_per_second + samp_rate + type float - title - Resampled signal (csdr) - - - samp_rate - samp_rate_2 - - - baseband_freq - 0 - - - y_per_div - 10 - - - y_divs - 10 - - - ref_level - 0 - - - ref_scale - 2.0 - - - fft_size - 1024 - - - fft_rate - 15 - - - peak_hold - False - - - average - False - - - avg_alpha - 0 - - - win - None - - - win_size - - - - grid_pos - - - - notebook - nb0,1 - - - freqvar - None - - - alias - - - - affinity - - - - _coordinate - (640, 235) - - - _rotation - 0 - - - - wxgui_scopesink2 - - id - wxgui_scopesink2_0 - - - _enabled - True - - - type - float - - - title - Resampled signal (csdr) - - - samp_rate - samp_rate_2 - - - v_scale - 0 - - - v_offset - 0 - - - t_scale - 0 - - - ac_couple - False - - - xy_mode - False - - - num_inputs + vlen 1 - - win_size - - - - grid_pos - - - - notebook - nb0,0 - - - trig_mode - wxgui.TRIG_MODE_AUTO - - - y_axis_label - Counts - - - alias - - - - affinity - - - - _coordinate - (640, 115) - - - _rotation - 0 - fractional_resampler_xx - id - fractional_resampler_xx_0 + alias + + + + comment + + + + affinity + _enabled True - type - float + _coordinate + (360, 608) + + + _rotation + 0 + + + id + fractional_resampler_xx_0 + + + maxoutbuf + 0 + + + minoutbuf + 0 phase_shift @@ -685,147 +373,248 @@ resamp_ratio decimation + + type + float + + + + ha5kfu_execproc_xx alias + + commandline + "csdr fractional_decimator_ff "+str(decimation) + + + comment + + affinity - minoutbuf + _enabled + True + + + _coordinate + (352, 115) + + + _rotation 0 + + id + ha5kfu_execproc_xx_0 + maxoutbuf 0 - _coordinate - (360, 608) + minoutbuf + 0 - _rotation - 0 + type + ff - wxgui_scopesink2 + notebook - id - wxgui_scopesink2_0_1 + alias + + + + comment + _enabled True - type - float + _coordinate + (16, 683) - title - Resampled signal (GNU Radio) - - - samp_rate - samp_rate_2 - - - v_scale + _rotation 0 - - v_offset - 0 - - - t_scale - 0 - - - ac_couple - False - - - xy_mode - False - - - num_inputs - 1 - - - win_size - - grid_pos - + 1,2,1,1 + + + id + nb0 + + + labels + ['scope','fft'] notebook - nb2,0 + - trig_mode - wxgui.TRIG_MODE_AUTO + style + wx.NB_TOP + + + + notebook + + alias + - y_axis_label - Counts + comment + + + + _enabled + True + + + _coordinate + (16, 579) + + + _rotation + 0 + + + grid_pos + 1,1,1,1 + + + id + nb1 + + + labels + ['scope','fft'] + + + notebook + + + + style + wx.NB_TOP + + + + notebook + + alias + + + + comment + + + + _enabled + True + + + _coordinate + (16, 475) + + + _rotation + 0 + + + grid_pos + 2,2,1,1 + + + id + nb2 + + + labels + ['scope','fft'] + + + notebook + + + + style + wx.NB_TOP + + + + wxgui_fftsink2 + + avg_alpha + 0 + + + average + False + + + baseband_freq + 0 alias + + comment + + affinity + + _enabled + True + + + fft_size + 1024 + + + freqvar + None + _coordinate - (640, 459) + (640, 235) _rotation 0 - - - wxgui_fftsink2 + + grid_pos + + id - wxgui_fftsink2_0_0 + wxgui_fftsink2_0 - _enabled - True + notebook + nb0,1 - type - float - - - title - Resampled signal (GNU Radio) - - - samp_rate - samp_rate_2 - - - baseband_freq - 0 - - - y_per_div - 10 - - - y_divs - 10 + peak_hold + False ref_level @@ -835,54 +624,77 @@ ref_scale 2.0 - - fft_size - 1024 - fft_rate 15 - peak_hold - False + samp_rate + samp_rate_2 - average - False + title + Resampled signal (csdr) - avg_alpha - 0 - - - win - None + type + float win_size - grid_pos - - - - notebook - nb2,1 - - - freqvar + win None + + y_divs + 10 + + + y_per_div + 10 + + + + wxgui_fftsink2 + + avg_alpha + 0 + + + average + False + + + baseband_freq + 0 + alias + + comment + + affinity + + _enabled + True + + + fft_size + 1024 + + + freqvar + None + _coordinate (640, 579) @@ -891,192 +703,436 @@ _rotation 0 - - - notebook - - id - nb0 - - - _enabled - True - - - style - wx.NB_TOP - - - labels - ['scope','fft'] - grid_pos - 1,2,1,1 + + + + id + wxgui_fftsink2_0_0 notebook + nb2,1 + + + peak_hold + False + + + ref_level + 0 + + + ref_scale + 2.0 + + + fft_rate + 15 + + + samp_rate + samp_rate_2 + + + title + Resampled signal (GNU Radio) + + + type + float + + + win_size + + win + None + + + y_divs + 10 + + + y_per_div + 10 + + + + wxgui_fftsink2 + + avg_alpha + 0 + + + average + False + + + baseband_freq + 0 + alias - _coordinate - (16, 683) + comment + - _rotation - 0 - - - - notebook - - id - nb1 + affinity + _enabled True - style - wx.NB_TOP + fft_size + 1024 - labels - ['scope','fft'] + freqvar + None + + + _coordinate + (360, 355) + + + _rotation + 0 grid_pos - 1,1,1,1 + + + + id + wxgui_fftsink2_1 notebook + nb1,1 + + + peak_hold + False + + + ref_level + 0 + + + ref_scale + 2.0 + + + fft_rate + 15 + + + samp_rate + samp_rate + + + title + Original signal + + + type + float + + + win_size + + win + None + + + y_divs + 10 + + + y_per_div + 10 + + + + wxgui_scopesink2 + + ac_couple + False + alias - _coordinate - (16, 579) + comment + - _rotation - 0 - - - - notebook - - id - nb2 + affinity + _enabled True - style - wx.NB_TOP + _coordinate + (640, 115) - labels - ['scope','fft'] + _rotation + 0 grid_pos - 2,2,1,1 + + + + id + wxgui_scopesink2_0 notebook + nb0,0 + + + num_inputs + 1 + + + samp_rate + samp_rate_2 + + + t_scale + 0 + + + title + Resampled signal (csdr) + + + trig_mode + wxgui.TRIG_MODE_AUTO + + + type + float + + + v_offset + 0 + + + v_scale + 0 + + + win_size + + xy_mode + False + + + y_axis_label + Counts + + + + wxgui_scopesink2 + + ac_couple + False + alias - _coordinate - (16, 475) + comment + - _rotation - 0 - - - - variable_slider - - id - input_freq + affinity + _enabled True - label - + _coordinate + (360, 235) - value - 4e3 - - - min + _rotation 0 - - max - samp_rate/2 - - - num_steps - 100 - - - style - wx.SL_HORIZONTAL - - - converver - float_converter - grid_pos + + id + wxgui_scopesink2_0_0 + notebook + nb1,0 + + + num_inputs + 1 + + + samp_rate + samp_rate + + + t_scale + 0 + + + title + Original signal + + + trig_mode + wxgui.TRIG_MODE_AUTO + + + type + float + + + v_offset + 0 + + + v_scale + 0 + + + win_size + + xy_mode + False + + + y_axis_label + Counts + + + + wxgui_scopesink2 + + ac_couple + False + alias + + comment + + + + affinity + + + + _enabled + True + _coordinate - (16, 355) + (640, 459) _rotation 0 + + grid_pos + + + + id + wxgui_scopesink2_0_1 + + + notebook + nb2,0 + + + num_inputs + 1 + + + samp_rate + samp_rate_2 + + + t_scale + 0 + + + title + Resampled signal (GNU Radio) + + + trig_mode + wxgui.TRIG_MODE_AUTO + + + type + float + + + v_offset + 0 + + + v_scale + 0 + + + win_size + + + + xy_mode + False + + + y_axis_label + Counts + - ha5kfu_execproc_xx_0 - wxgui_scopesink2_0 + analog_sig_source_x_0 + blocks_throttle_0 0 0 - analog_sig_source_x_0 - blocks_throttle_0 + blocks_throttle_0 + fractional_resampler_xx_0 0 0 @@ -1086,18 +1142,6 @@ 0 0 - - blocks_throttle_0 - wxgui_scopesink2_0_0 - 0 - 0 - - - ha5kfu_execproc_xx_0 - wxgui_fftsink2_0 - 0 - 0 - blocks_throttle_0 wxgui_fftsink2_1 @@ -1106,7 +1150,13 @@ blocks_throttle_0 - fractional_resampler_xx_0 + wxgui_scopesink2_0_0 + 0 + 0 + + + fractional_resampler_xx_0 + wxgui_fftsink2_0_0 0 0 @@ -1117,8 +1167,14 @@ 0 - fractional_resampler_xx_0 - wxgui_fftsink2_0_0 + ha5kfu_execproc_xx_0 + wxgui_fftsink2_0 + 0 + 0 + + + ha5kfu_execproc_xx_0 + wxgui_scopesink2_0 0 0 diff --git a/libcsdr.c b/libcsdr.c index 972d216..6b8fc62 100644 --- a/libcsdr.c +++ b/libcsdr.c @@ -654,7 +654,7 @@ float inline fir_one_pass_ff(float* input, float* taps, int taps_length) return acc; } -fractional_decimator_ff_t fractional_decimator_ff(float* input, float* output, int input_size, float rate, float *taps, int taps_length, fractional_decimator_ff_t d) +old_fractional_decimator_ff_t old_fractional_decimator_ff(float* input, float* output, int input_size, float rate, float *taps, int taps_length, old_fractional_decimator_ff_t d) { if(rate<=1.0) return d; //sanity check, can't decimate <=1.0 //This routine can handle floating point decimation rates. @@ -687,6 +687,98 @@ fractional_decimator_ff_t fractional_decimator_ff(float* input, float* output, i return d; } +fractional_decimator_ff_t fractional_decimator_ff_init(float rate, int num_poly_points, float* taps, int taps_length) +{ + fractional_decimator_ff_t d; + d.num_poly_points = num_poly_points&~1; //num_poly_points needs to be even! + d.poly_precalc_denomiator = (float*)malloc(d.num_poly_points*sizeof(float)); + //x0..x3 + //-1,0,1,2 + //-(4/2)+1 + //x0..x5 + //-2,-1,0,1,2,3 + d.xifirst=-(num_poly_points/2)+1, d.xilast=num_poly_points/2; + int id = 0; //index in poly_precalc_denomiator + for(int xi=d.xifirst;xi<=d.xilast;xi++) + { + d.poly_precalc_denomiator[id]=1; + for(int xj=d.xifirst;xj<=d.xilast;xj++) + { + if(xi!=xj) d.poly_precalc_denomiator[id] *= (xi-xj); //poly_precalc_denomiator could be integer as well. But that would later add a necessary conversion. + } + id++; + } + d.where=-d.xifirst; + d.coeffs_buf=(float*)malloc(d.num_poly_points*sizeof(float)); + d.filtered_buf=(float*)malloc(d.num_poly_points*sizeof(float)); + //d.last_inputs_circbuf = (float)malloc(d.num_poly_points*sizeof(float)); + //d.last_inputs_startsat = 0; + //d.last_inputs_samplewhere = -1; + //for(int i=0;irate); + if(DEBUG_ASSERT) assert(d->rate > 1.0); + int oi=0; //output index + int index_high; +#define FD_INDEX_LOW (index_high-1) + //we optimize to calculate ceilf(where) only once every iteration, so we do it here: + for(;(index_high=ceilf(d->where))+d->taps_lengthwhere+=d->rate) //@fractional_decimator_ff + { + int sxifirst = FD_INDEX_LOW + d->xifirst; + int sxilast = FD_INDEX_LOW + d->xilast; + if(d->taps) for(int wi=0;winum_poly_points;wi++) d->filtered_buf[wi] = fir_one_pass_ff(input+FD_INDEX_LOW+wi, d->taps, d->taps_length); + else for(int wi=0;winum_poly_points;wi++) d->filtered_buf[wi] = *(input+FD_INDEX_LOW+wi); + int id=0; + float xwhere = d->where - FD_INDEX_LOW; + for(int xi=d->xifirst;xi<=d->xilast;xi++) + { + d->coeffs_buf[id]=1; + for(int xj=d->xifirst;xj<=d->xilast;xj++) + { + if(xi!=xj) d->coeffs_buf[id] *= (xwhere-xj); + } + id++; + } + float acc = 0; + for(int i=0;inum_poly_points;i++) + { + acc += (d->coeffs_buf[i]/d->poly_precalc_denomiator[i])*d->filtered_buf[i]; //(xnom/xden)*yn + } + output[oi++]=acc; + } + d->input_processed=FD_INDEX_LOW + d->xifirst; + d->where-=d->input_processed; + if(DEBUG_ASSERT) assert(d->where >= -d->xifirst); + d->output_size=oi; +} + +/* + int last_input_samplewhere_shouldbe = (index_high-1)+xifirst; + int last_input_offset = last_input_samplewhere_shouldbe - d->last_input_samplewhere; + if(last_input_offset < num_poly_points) + { + //if we can move the last_input circular buffer, we move, and add the new samples at the end + d->last_inputs_startsat += last_input_offset; + d->last_inputs_startsat %= num_poly_points; + int num_copied_samples = 0; + for(int i=0; ilast_inputs_circbuf[i]= + } + d->last_input_samplewhere = d->las + } +*/ void apply_fir_fft_cc(FFT_PLAN_T* plan, FFT_PLAN_T* plan_inverse, complexf* taps_fft, complexf* last_overlap, int overlap_size) { diff --git a/libcsdr.h b/libcsdr.h index 4f45a43..b1afa1c 100644 --- a/libcsdr.h +++ b/libcsdr.h @@ -146,12 +146,33 @@ void accumulate_power_cf(complexf* input, float* output, int size); void log_ff(float* input, float* output, int size, float add_db); typedef struct fractional_decimator_ff_s +{ + float where; + int input_processed; + int output_size; + int num_poly_points; //number of samples that the Lagrange interpolator will use + float* poly_precalc_denomiator; //while we don't precalculate coefficients here as in a Farrow structure, because it is a fractional interpolator, but we rather precaculate part of the interpolator expression + //float* last_inputs_circbuf; //circular buffer to store the last (num_poly_points) number of input samples. + //int last_inputs_startsat; //where the circular buffer starts now + //int last_inputs_samplewhere; + float* coeffs_buf; + float* filtered_buf; + int xifirst; + int xilast; + float rate; + float *taps; + int taps_length; +} fractional_decimator_ff_t; +fractional_decimator_ff_t fractional_decimator_ff_init(float rate, int num_poly_points, float* taps, int taps_length); +void fractional_decimator_ff(float* input, float* output, int input_size, fractional_decimator_ff_t* d); + +typedef struct old_fractional_decimator_ff_s { float remain; int input_processed; int output_size; -} fractional_decimator_ff_t; -fractional_decimator_ff_t fractional_decimator_ff(float* input, float* output, int input_size, float rate, float *taps, int taps_length, fractional_decimator_ff_t d); +} old_fractional_decimator_ff_t; +old_fractional_decimator_ff_t old_fractional_decimator_ff(float* input, float* output, int input_size, float rate, float *taps, int taps_length, old_fractional_decimator_ff_t d); typedef struct shift_table_data_s {