Skip to content

Commit

Permalink
clean up spectrum mode, enforce limits on bin counts
Browse files Browse the repository at this point in the history 
Still needs work
  • Loading branch information
@ka9q
ka9q committed Apr 30, 2024
1 parent 3e48a57 commit 897233c
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 54 deletions.
78 changes: 41 additions & 37 deletions filter.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,8 @@ struct filter_in *create_filter_input(struct filter_in *master,int const L,int c
// Set up output (slave) side of filter (possibly one of several sharing the same input master)
// These output filters should be deleted before their masters
// Segfault will occur if filter_in is deleted and execute_filter_output is executed
struct filter_out *create_filter_output(struct filter_out *slave,struct filter_in * master,complex float * const response,int olen, enum filtertype const out_type){
// Special case: for type == SPECTRUM, 'len' is the number of FFT bins, not the number of output time domain points (since there aren't any)
struct filter_out *create_filter_output(struct filter_out *slave,struct filter_in * master,complex float * const response,int len, enum filtertype const out_type){
assert(master != NULL);
if(master == NULL)
return NULL;
Expand All @@ -225,61 +226,63 @@ struct filter_out *create_filter_output(struct filter_out *slave,struct filter_i
if(slave == NULL)
return NULL;

assert(olen > 0);
if(olen > master->ilen)
olen = master->ilen; // Interpolation not yet supported - are callers prepared for this?
assert(len > 0);

// Share all but output fft bins, response, output and output type
slave->master = master;
slave->out_type = out_type;
slave->olen = olen;

float const overlap = (float)(master->ilen + master->impulse_length - 1) / master->ilen; // Total FFT time points / used time points
int const osize = round(olen * overlap); // Total number of time-domain FFT points including overlap

// N / L = Total FFT points / time domain points
float const overlap = (float)(master->ilen + master->impulse_length - 1) / master->ilen;
slave->response = response;
if(response != NULL)
slave->noise_gain = noise_gain(slave);
else
slave->noise_gain = NAN;
slave->noise_gain = (response == NULL) ? NAN : noise_gain(slave);

pthread_mutex_lock(&FFTW_planning_mutex);
fftwf_plan_with_nthreads(1); // IFFTs are always small, use only one internal thread
switch(slave->out_type){
default:
case COMPLEX:
case CROSS_CONJ:
slave->bins = osize; // Same as total number of time domain points
slave->fdomain = lmalloc(sizeof(complex float) * slave->bins);
slave->output_buffer.c = lmalloc(sizeof(complex float) * osize);
assert(slave->output_buffer.c != NULL);
slave->output_buffer.r = NULL; // catch erroneous references
slave->output.c = slave->output_buffer.c + osize - olen;
if((slave->rev_plan = fftwf_plan_dft_1d(osize,slave->fdomain,slave->output_buffer.c,FFTW_BACKWARD,FFTW_WISDOM_ONLY|FFTW_planning_level)) == NULL){
suggest(FFTW_planning_level,osize,FFTW_BACKWARD,COMPLEX);
slave->rev_plan = fftwf_plan_dft_1d(osize,slave->fdomain,slave->output_buffer.c,FFTW_BACKWARD,FFTW_MEASURE);
{
slave->olen = len;
slave->bins = round(len * overlap); // Total number of time-domain FFT points including overlap
slave->fdomain = lmalloc(sizeof(complex float) * slave->bins);
slave->output_buffer.c = lmalloc(sizeof(complex float) * slave->bins);
assert(slave->output_buffer.c != NULL);
slave->output_buffer.r = NULL; // catch erroneous references
slave->output.c = slave->output_buffer.c + slave->bins - len;
if((slave->rev_plan = fftwf_plan_dft_1d(slave->bins,slave->fdomain,slave->output_buffer.c,FFTW_BACKWARD,FFTW_WISDOM_ONLY|FFTW_planning_level)) == NULL){
suggest(FFTW_planning_level,slave->bins,FFTW_BACKWARD,COMPLEX);
slave->rev_plan = fftwf_plan_dft_1d(slave->bins,slave->fdomain,slave->output_buffer.c,FFTW_BACKWARD,FFTW_MEASURE);
}
}
if(fftwf_export_wisdom_to_filename(Wisdom_file) == 0)
fprintf(stdout,"fftwf_export_wisdom_to_filename(%s) failed\n",Wisdom_file);
break;
case SPECTRUM: // Like complex, but no IFFT or output time domain buffer
slave->bins = osize;
slave->fdomain = lmalloc(sizeof(complex float) * slave->bins); // User reads this directly
assert(slave->fdomain != NULL);
// Note: No time domain buffer; slave->output, etc, all NULL
// Also don't set up an IFFT
{
slave->olen = 0;
slave->bins = len;
slave->fdomain = lmalloc(sizeof(complex float) * slave->bins); // User reads this directly
assert(slave->fdomain != NULL);
// Note: No time domain buffer; slave->output, etc, all NULL
// Also don't set up an IFFT
}
break;
case REAL:
slave->bins = osize / 2 + 1;
slave->fdomain = lmalloc(sizeof(complex float) * slave->bins); // Not really needed for SPECTRUM?
assert(slave->fdomain != NULL);
slave->output_buffer.r = lmalloc(sizeof(float) * osize);
assert(slave->output_buffer.r != NULL);
slave->output_buffer.c = NULL;
slave->output.r = slave->output_buffer.r + osize - olen;
if((slave->rev_plan = fftwf_plan_dft_c2r_1d(osize,slave->fdomain,slave->output_buffer.r,FFTW_WISDOM_ONLY|FFTW_planning_level)) == NULL){
suggest(FFTW_planning_level,osize,FFTW_BACKWARD,REAL);
slave->rev_plan = fftwf_plan_dft_c2r_1d(osize,slave->fdomain,slave->output_buffer.r,FFTW_MEASURE);
{
slave->olen = len;
slave->bins = round(len * overlap) / 2 + 1;
slave->fdomain = lmalloc(sizeof(complex float) * slave->bins);
assert(slave->fdomain != NULL);
slave->output_buffer.r = lmalloc(sizeof(float) * slave->bins);
assert(slave->output_buffer.r != NULL);
slave->output_buffer.c = NULL;
slave->output.r = slave->output_buffer.r + slave->bins - len;
if((slave->rev_plan = fftwf_plan_dft_c2r_1d(slave->bins,slave->fdomain,slave->output_buffer.r,FFTW_WISDOM_ONLY|FFTW_planning_level)) == NULL){
suggest(FFTW_planning_level,slave->bins,FFTW_BACKWARD,REAL);
slave->rev_plan = fftwf_plan_dft_c2r_1d(slave->bins,slave->fdomain,slave->output_buffer.r,FFTW_MEASURE);
}
}
if(fftwf_export_wisdom_to_filename(Wisdom_file) == 0)
fprintf(stdout,"fftwf_export_wisdom_to_filename(%s) failed\n",Wisdom_file);
Expand Down Expand Up @@ -751,6 +754,7 @@ int window_filter(int const L,int const M,complex float * const response,float c
memcpy(buffer,response,N * sizeof(*buffer));
fftwf_execute(rev_filter_plan);
fftwf_destroy_plan(rev_filter_plan);
rev_filter_plan = NULL;
#ifdef FILTER_DEBUG
fprintf(stderr,"window_filter raw time domain\n");
for(int n=0; n < N; n++){
Expand Down Expand Up @@ -783,7 +787,7 @@ int window_filter(int const L,int const M,complex float * const response,float c
// Now back to frequency domain
fftwf_execute(fwd_filter_plan);
fftwf_destroy_plan(fwd_filter_plan);

fwd_filter_plan = NULL;
#ifdef FILTER_DEBUG
fprintf(stderr,"window_filter filter response amplitude\n");
for(int n=0;n<N;n++)
Expand Down
42 changes: 25 additions & 17 deletions spectrum.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@
// Spectrum analysis thread
void *demod_spectrum(void *arg){
assert(arg != NULL);
struct channel * const chan = arg;
struct channel * const chan = arg;

{
char name[100];
snprintf(name,sizeof(name),"spect %u",chan->output.rtp.ssrc);
Expand All @@ -29,7 +29,6 @@ void *demod_spectrum(void *arg){
FREE(chan->filter.energies);
FREE(chan->spectrum.bin_data);
delete_filter_output(&chan->filter.out);
pthread_mutex_unlock(&chan->status.lock);

if(chan->spectrum.bin_count <= 0)
chan->spectrum.bin_count = 64; // Force a reasonable number of bins
Expand All @@ -46,39 +45,48 @@ void *demod_spectrum(void *arg){
int const t = roundf(chan->spectrum.bin_bw / fft_bin_spacing);
int const binsperbin = (t == 0) ? 1 : t; // Force reasonable value
chan->spectrum.bin_bw = binsperbin * fft_bin_spacing; // Force to integer multiple of fft_bin_spacing
int const fft_bins = chan->spectrum.bin_count * binsperbin;
int fft_bins = chan->spectrum.bin_count * binsperbin;
if(fft_bins > Frontend.in.bins){
// Too many, limit to total available
fft_bins = Frontend.in.bins;
chan->spectrum.bin_count = fft_bins / binsperbin;
}

// Special filter without a response curve or IFFT
if(create_filter_output(&chan->filter.out,&Frontend.in,NULL,fft_bins,SPECTRUM) == NULL)
assert(0);

// Although we don't use filter_output, chan->filter.min_IF and max_IF still need to be set
// so radio.c:set_freq() will set the front end tuner properly
chan->filter.max_IF = (chan->spectrum.bin_count * chan->spectrum.bin_bw)/2;
chan->filter.min_IF = -chan->filter.max_IF;

// Special filter without a response curve or IFFT
// the output size arg to create_filter_output refers only to usable output time points. We want to access the frequency domain
// points directly, so we decrease to correct for the overlap factor
// I know all this can be simplified
int const olen = fft_bins * (float)Frontend.L / N;
create_filter_output(&chan->filter.out,&Frontend.in,NULL,olen,SPECTRUM);

// If it's already allocated (why should it be?) we don't know how big it is
if(chan->spectrum.bin_data != NULL)
FREE(chan->spectrum.bin_data);
chan->spectrum.bin_data = calloc(chan->spectrum.bin_count,sizeof(*chan->spectrum.bin_data));
// chan->spectrum.bin_data = calloc(chan->spectrum.bin_count,sizeof(*chan->spectrum.bin_data));
// experiment - make array largest possible to temp avoid memory corruption
chan->spectrum.bin_data = calloc(Frontend.in.bins,sizeof(*chan->spectrum.bin_data));

set_freq(chan,chan->tune.freq); // retune front end if needed to cover requested bandwidth
pthread_mutex_unlock(&chan->status.lock);

if(Verbose > 1)
fprintf(stdout,"spectrum %d: freq %'lf bin_bw %'f binsperbin %'d bin_count %'d\n",chan->output.rtp.ssrc,chan->tune.freq,chan->spectrum.bin_bw,binsperbin,chan->spectrum.bin_count);

// Still need to clean up code to force radio freq to be multiple of FFT bin spacing
while(downconvert(chan) == 0){
int binp = 0;
int binp = 0;
for(int i=0; i < chan->spectrum.bin_count; i++){ // For each noncoherent integration bin above center freq
double p = 0;
for(int j=0; j < binsperbin; j++){ // Add energy of each fft bin that's part of this user integration bin
for(int j=0; j < binsperbin; j++) // Add energy of each fft bin that's part of this user integration bin
p += cnrmf(chan->filter.out.fdomain[binp++]);
}

// Accumulate energy until next poll
chan->spectrum.bin_data[i] += p;
}
}
FREE(chan->spectrum.bin_data);
FREE(chan->status.command);
FREE(chan->filter.energies);
delete_filter_output(&chan->filter.out);
return NULL;
}

0 comments on commit 897233c

Please sign in to comment.