Skip to content

Commit edc2c21

Browse files
ha7ilmha7ilm
committed
FastDDC code added, and now it builds!
1 parent cb1b6ac commit edc2c21

File tree

4 files changed

+254
-1
lines changed

4 files changed

+254
-1
lines changed

‎csdr.c

+101-1
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,8 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
4848
#include "ima_adpcm.h"
4949
#include <sched.h>
5050
#include <math.h>
51+
#include <strings.h>
52+
#include "fastddc.h"
5153

5254
char usage[]=
5355
"csdr - a simple commandline tool for Software Defined Radio receiver DSP.\n\n"
@@ -1269,7 +1271,6 @@ int main(int argc, char *argv[])
12691271
float high_cut;
12701272
float transition_bw;
12711273
window_t window = WINDOW_DEFAULT;
1272-
char window_string[256]; //TODO: nice buffer overflow opportunity
12731274

12741275
int fd;
12751276
if(fd=init_fifo(argc,argv))
@@ -1646,6 +1647,105 @@ int main(int argc, char *argv[])
16461647
}
16471648
}
16481649

1650+
if( !strcmp(argv[1],"fastddc_fwd_cc") ) //<decimation> [transition_bw [window]]
1651+
{
1652+
int decimation;
1653+
if(argc<=2) return badsyntax("need required parameter (decimation)");
1654+
sscanf(argv[2],"%d",&decimation);
1655+
1656+
float transition_bw = 0.05;
1657+
if(argc>=3) sscanf(argv[3],"%g",&transition_bw);
1658+
1659+
window_t window = WINDOW_DEFAULT;
1660+
if(argc>=4) window=firdes_get_window_from_string(argv[5]);
1661+
else fprintf(stderr,"fastddc_fwd_cc: window = %s\n",firdes_get_string_from_window(window));
1662+
1663+
fastddc_t ddc;
1664+
if(fastddc_init(&ddc, transition_bw, decimation, 0)) { badsyntax("error in fastddc_init()"); return 1; }
1665+
fastddc_print(&ddc);
1666+
1667+
if(!initialize_buffers()) return -2;
1668+
sendbufsize(ddc.fft_size);
1669+
1670+
//make FFT plan
1671+
complexf* input = (complexf*)fft_malloc(sizeof(complexf)*ddc.fft_size);
1672+
complexf* windowed = (complexf*)fft_malloc(sizeof(complexf)*ddc.fft_size);
1673+
complexf* output = (complexf*)fft_malloc(sizeof(complexf)*ddc.fft_size);
1674+
1675+
for(int i=0;i<ddc.fft_size;i++) iof(input,i)=qof(input,i)=0; //null the input buffer
1676+
1677+
int benchmark = 1;
1678+
if(benchmark) fprintf(stderr,"fastddc_fwd_cc: benchmarking FFT...");
1679+
FFT_PLAN_T* plan=make_fft_c2c(ddc.fft_size, windowed, output, 1, benchmark);
1680+
if(benchmark) fprintf(stderr," done\n");
1681+
1682+
for(;;)
1683+
{
1684+
FEOF_CHECK;
1685+
//overlapped FFT
1686+
for(int i=0;i<ddc.overlap_length;i++) input[i]=input[i+ddc.input_size];
1687+
fread(input+ddc.overlap_length, sizeof(complexf), ddc.input_size, stdin);
1688+
apply_window_c(input,windowed,ddc.fft_size,window);
1689+
fft_execute(plan);
1690+
fwrite(output, sizeof(complexf), ddc.fft_size, stdout);
1691+
TRY_YIELD;
1692+
}
1693+
}
1694+
1695+
if( !strcmp(argv[1],"fastddc_apply_cc") ) //<decimation> <shift_rate> [transition_bw [window]]
1696+
{
1697+
int decimation;
1698+
if(argc<=2) return badsyntax("need required parameter (decimation)");
1699+
sscanf(argv[2],"%d",&decimation);
1700+
1701+
float shift_rate;
1702+
if(argc>=3) sscanf(argv[3],"%g",&shift_rate);
1703+
1704+
float transition_bw = 0.05;
1705+
if(argc>=4) sscanf(argv[4],"%g",&transition_bw);
1706+
1707+
window_t window = WINDOW_DEFAULT;
1708+
if(argc>=5) window=firdes_get_window_from_string(argv[5]);
1709+
else fprintf(stderr,"fastddc_apply_cc: window = %s\n",firdes_get_string_from_window(window));
1710+
1711+
fastddc_t ddc;
1712+
if(fastddc_init(&ddc, transition_bw, decimation, shift_rate)) { badsyntax("error in fastddc_init()"); return 1; }
1713+
fastddc_print(&ddc);
1714+
1715+
if(!initialize_buffers()) return -2;
1716+
sendbufsize(ddc.output_size);
1717+
1718+
//prepare making the filter and doing FFT on it
1719+
complexf* taps=(complexf*)calloc(sizeof(complexf),ddc.fft_size); //initialize to zero
1720+
complexf* taps_fft=(complexf*)malloc(sizeof(complexf)*ddc.fft_size);
1721+
FFT_PLAN_T* plan_taps = make_fft_c2c(ddc.fft_size, taps, taps_fft, 1, 0); //forward, don't benchmark (we need this only once)
1722+
1723+
//make the filter
1724+
float filter_half_bw = 0.5/decimation;
1725+
firdes_bandpass_c(taps, ddc.taps_real_length, shift_rate-filter_half_bw, shift_rate+filter_half_bw, window);
1726+
fft_execute(plan_taps);
1727+
1728+
//make FFT plan
1729+
complexf* input = (complexf*)fft_malloc(sizeof(complexf)*ddc.fft_size);
1730+
complexf* output = (complexf*)fft_malloc(sizeof(complexf)*ddc.output_size);
1731+
1732+
int benchmark = 1;
1733+
if(benchmark) fprintf(stderr,"fastddc_apply_cc: benchmarking FFT...");
1734+
FFT_PLAN_T* plan_inverse = make_fft_c2c(ddc.fft_size, input, output, 0, 1); //inverse, do benchmark
1735+
if(benchmark) fprintf(stderr," done\n");
1736+
1737+
decimating_shift_addition_status_t shift_stat;
1738+
bzero(&shift_stat, sizeof(shift_stat));
1739+
for(;;)
1740+
{
1741+
FEOF_CHECK;
1742+
fread(input, sizeof(complexf), ddc.fft_size, stdin);
1743+
shift_stat = fastddc_apply_cc(input, output, &ddc, plan_inverse, taps_fft, shift_stat);
1744+
fwrite(output, sizeof(complexf), ddc.output_size, stdout);
1745+
TRY_YIELD;
1746+
}
1747+
}
1748+
16491749
if(!strcmp(argv[1],"none"))
16501750
{
16511751
return 0;

‎fastddc.c

+125
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
/*
2+
This software is part of libcsdr, a set of simple DSP routines for
3+
Software Defined Radio.
4+
5+
Copyright (c) 2014, Andras Retzler <[email protected]>
6+
All rights reserved.
7+
8+
Redistribution and use in source and binary forms, with or without
9+
modification, are permitted provided that the following conditions are met:
10+
* Redistributions of source code must retain the above copyright
11+
notice, this list of conditions and the following disclaimer.
12+
* Redistributions in binary form must reproduce the above copyright
13+
notice, this list of conditions and the following disclaimer in the
14+
documentation and/or other materials provided with the distribution.
15+
* Neither the name of the copyright holder nor the
16+
names of its contributors may be used to endorse or promote products
17+
derived from this software without specific prior written permission.
18+
19+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
20+
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21+
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22+
DISCLAIMED. IN NO EVENT SHALL ANDRAS RETZLER BE LIABLE FOR ANY
23+
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24+
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25+
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
26+
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27+
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28+
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29+
*/
30+
31+
#include "fastddc.h"
32+
33+
//DDC implementation based on:
34+
//http://www.3db-labs.com/01598092_MultibandFilterbank.pdf
35+
36+
inline int is_integer(float a) { return floorf(a) == a; }
37+
38+
int fastddc_init(fastddc_t* ddc, float transition_bw, int decimation, float shift_rate)
39+
{
40+
ddc->pre_decimation = 1; //this will be done in the frequency domain
41+
ddc->post_decimation = decimation; //this will be done in the time domain
42+
while( is_integer((float)ddc->post_decimation/2) && ddc->post_decimation/2 != 1)
43+
{
44+
ddc->post_decimation/=2;
45+
ddc->pre_decimation*=2;
46+
}
47+
ddc->taps_real_length = firdes_filter_len(transition_bw); //the number of non-zero taps
48+
ddc->taps_length = ceil(ddc->taps_real_length/(float)ddc->pre_decimation) * ddc->pre_decimation; //the number of taps must be a multiple of the decimation factor
49+
ddc->fft_size = next_pow2(ddc->taps_length * 4); //it is a good rule of thumb for performance (based on the article), but we should do benchmarks
50+
while (ddc->fft_size<ddc->pre_decimation) ddc->fft_size*=2; //fft_size should be a multiple of pre_decimation
51+
ddc->overlap_length = ddc->taps_length - 1;
52+
ddc->input_size = ddc->fft_size - ddc->overlap_length;
53+
ddc->fft_inv_size = ddc->fft_size / ddc->pre_decimation;
54+
55+
//Shift operation in the frequency domain: we can shift by a multiple of v.
56+
ddc->v = ddc->fft_size/ddc->overlap_length; //+-1 ? (or maybe ceil() this?) //TODO: why?
57+
int middlebin=ddc->fft_size / 2;
58+
ddc->startbin = middlebin + middlebin * shift_rate * 2;
59+
ddc->startbin = ddc->v * round( ddc->startbin / (float)ddc->v );
60+
ddc->offsetbin = ddc->startbin - middlebin;
61+
ddc->post_shift = ((float)ddc->offsetbin/ddc->fft_size) - shift_rate;
62+
ddc->pre_shift = ddc->offsetbin * ddc->v;
63+
64+
//Overlap is scraped, not added
65+
ddc->scrape=ddc->overlap_length/ddc->pre_decimation;
66+
ddc->output_size=ddc->fft_inv_size-ddc->scrape;
67+
68+
return ddc->fft_size<=2; //returns true on error
69+
}
70+
71+
72+
void fastddc_print(fastddc_t* ddc)
73+
{
74+
fprintf(stderr,
75+
"fastddc_print_sizes(): (fft_size = %d) = (taps_length = %d) + (input_size = %d) - 1\n"
76+
"\t(overlap_length = %d) = taps_length - 1, taps_real_length = %d\n"
77+
"\tdecimation = (pre_decimation = %d) * (post_decimation = %d), fft_inv_size = %d\n"
78+
"\tstartbin = %d, offsetbin = %d, v = %d, pre_shift = %g, post_shift = %g\n"
79+
,
80+
ddc->fft_size, ddc->taps_length, ddc->input_size,
81+
ddc->overlap_length, ddc->taps_real_length,
82+
ddc->pre_decimation, ddc->post_decimation, ddc->fft_inv_size,
83+
ddc->startbin, ddc->offsetbin, ddc->v, ddc->pre_shift, ddc->post_shift );
84+
}
85+
86+
decimating_shift_addition_status_t fastddc_apply_cc(complexf* input, complexf* output, fastddc_t* ddc, FFT_PLAN_T* plan_inverse, complexf* taps_fft, decimating_shift_addition_status_t shift_stat)
87+
{
88+
//implements DDC by using the overlap & scrape method
89+
//TODO: +/-1s on overlap_size et al
90+
//input shoud have ddc->fft_size number of elements
91+
92+
complexf* inv_input = plan_inverse->input;
93+
complexf* inv_output = plan_inverse->output;
94+
95+
//Initialize buffers for inverse FFT to zero
96+
for(int i=0;i<plan_inverse->size;i++)
97+
{
98+
iof(inv_input,i)=0;
99+
qof(inv_input,i)=0;
100+
}
101+
102+
//Alias & shift & filter at once
103+
// * no, we won't break this algorithm to parts that are easier to understand: now we go for speed
104+
for(int i=0;i<ddc->fft_size;i++)
105+
{
106+
int output_index = (ddc->startbin+i)%plan_inverse->size;
107+
int tap_index = (ddc->fft_size+i-ddc->offsetbin)%ddc->fft_size;
108+
cmultadd(inv_input+output_index, input+i, taps_fft+tap_index); //cmultadd(output, input1, input2): complex output += complex input1 * complex input 2
109+
}
110+
111+
fft_execute(plan_inverse);
112+
113+
//Normalize data
114+
for(int i=0;i<plan_inverse->size;i++) //@apply_ddc_fft_cc: normalize by size
115+
{
116+
iof(inv_output,i)/=plan_inverse->size;
117+
qof(inv_output,i)/=plan_inverse->size;
118+
}
119+
120+
//Overlap is scraped, not added
121+
//Shift correction
122+
shift_addition_data_t dsadata=decimating_shift_addition_init(ddc->post_shift, ddc->post_decimation); //this could be optimized (passed as parameter), but we would not win too much at all
123+
shift_stat=decimating_shift_addition_cc(plan_inverse->output+ddc->scrape, output, ddc->output_size, dsadata, ddc->post_decimation, shift_stat);
124+
return shift_stat;
125+
}

‎fastddc.h

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
#include <math.h>
2+
#include "libcsdr.h"
3+
#include "libcsdr_gpl.h"
4+
5+
typedef struct fastddc_s
6+
{
7+
int pre_decimation;
8+
int post_decimation;
9+
int taps_length;
10+
int taps_real_length;
11+
int overlap_length; //it is taps_length - 1
12+
int fft_size;
13+
int fft_inv_size;
14+
int input_size;
15+
int output_size;
16+
float pre_shift;
17+
int startbin; //for pre_shift
18+
int v; //step for pre_shift
19+
int offsetbin;
20+
float post_shift;
21+
int output_scrape;
22+
int scrape;
23+
} fastddc_t;
24+
25+
int fastddc_init(fastddc_t* ddc, float transition_bw, int decimation, float shift_rate);
26+
decimating_shift_addition_status_t fastddc_apply_cc(complexf* input, complexf* output, fastddc_t* ddc, FFT_PLAN_T* plan_inverse, complexf* taps_fft, decimating_shift_addition_status_t shift_stat);
27+
void fastddc_print(fastddc_t* ddc);

‎libcsdr_wrapper.c

+1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
#include "libcsdr.c"
22
#include "libcsdr_gpl.c"
33
#include "ima_adpcm.c"
4+
#include "fastddc.c"
45
//this wrapper helps parsevect.py to generate better output

0 commit comments

Comments
 (0)