vladimir/wdsp
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
cpu-optimization
wdsp/analyzer.c
T
Uladzimir Karpenka 12dc701604 analyzer: block the dispatcher on a semaphore, stop spawning a thread per frame 
Two pieces of pure overhead, neither of them doing any DSP.

The POSIX QueueUserWorkItem() shim spawned a thread and immediately joined
it. That is a synchronous call -- the parallelism the Windows thread pool
provides is absent here either way -- so it only bought a pthread_create()
and its stack mmap, ~15 us, for every FFT frame. Call the function directly.
(The Windows build is untouched and still gets its thread pool; a real pool
for POSIX would be a separate change, and only pays off for num_stitch > 1.)

The dispatcher thread ran `for (each ss, LO) {...} Sleep(1);`, and Sleep(1)
is usleep(1000), so it woke 1000 times a second to re-read the same flags.
That burned 0.63% of a core even with no samples arriving. Give it a
semaphore instead, signalled by the four Spectrum*() entry points when new
samples land, and by SetAnalyzer/DestroyAnalyzer after they raise
end_dispatcher so the blocking wait always has a way out. SetAnalyzer holds
SetAnalyzerSection while it waits for the dispatcher to quit, and the
dispatcher never takes that section, so signalling from under it is safe.

Measured on an Apple M1 Pro; 16384-point complex FFT, 1024-sample buffers,
~40 pixel frames/s, CPU of all threads via getrusage:

    idle, no samples at all      0.63% of a core  ->  0.00%
    under load                   3.54%            ->  2.65%

Output pixels are bit-identical. Stressed with 3 create/destroy cycles and
24 on-the-fly SetAnalyzer reconfigurations while a second thread fed samples:
no deadlock.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
2026-07-09 23:59:20 +03:00

1901 lines
53 KiB
C
Raw Permalink Blame History

/* analyzer.c
This file is part of a program that implements a Spectrum Analyzer
used in conjunction with software-defined-radio hardware.
Copyright (C) 2012, 2013, 2014, 2016, 2023, 2025 Warren Pratt, NR0V
Copyright (C) 2012 David McQuate, WA8YWQ - Kaiser window & Bessel function added.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
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 General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
The author can be reached by email at
warren@wpratt.com
*/
#include "comm.h"
DP pdisp[dMAX_DISPLAYS];
double bessi0(double x)
{
double ax,ans;
double y;
if ((ax=fabs(x)) < 3.75) {
y = x / 3.75, y = y * y;
ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y* (1.2067492
+ y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
} else {
y = 3.75 / ax;
ans = (exp(ax) / sqrt(ax)) * (0.39894228 + y * (0.1328592e-1
+ y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
+ y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
+ y * 0.392377e-2))))))));
}
return ans;
}
void new_window(int disp, int type, int size, double PiAlpha)
{
DP a = pdisp[disp];
int i;
double arg0, arg1, cgsum, igsum;
switch (type)
{
case 0: // rectangular window
{
a->inv_coherent_gain = 1.0;
igsum = (double)size;
for (i = 0; i < size; i++)
a->window[i] = a->inv_coherent_gain * 1.0;
break;
}
case 1: // blackman-harris window (4 term)
{
arg0 = 2.0 * PI / ((double)size - 1.0);
cgsum = 0.0;
igsum = 0.0;
for (i = 0; i < size; i++)
{
arg1 = arg0 * (double)i;
a->window[i] = 0.35875 - 0.48829 * cos(arg1) + 0.14128 * cos(2.0 * arg1) - 0.01168 * cos(3.0 * arg1);
cgsum += a->window[i];
igsum += a->window[i] * a->window[i];
}
a->inv_coherent_gain = (double)size / cgsum;
for (i = 0; i < size; i++)
a->window[i] *= a->inv_coherent_gain;
break;
}
case 2: // hann window
{
arg0 = 2.0 * PI / ((double)size - 1.0);
cgsum = 0.0;
igsum = 0.0;
for (i = 0; i < size; i++)
{
a->window[i] = 0.5 * (1.0 - cos((double)i * arg0));
cgsum += a->window[i];
igsum += a->window[i] * a->window[i];
}
a->inv_coherent_gain = (double)size / cgsum;
for (i = 0; i < size; i++)
a->window[i] *= a->inv_coherent_gain;
break;
}
case 3: // flat-top window
{
arg0 = 2.0 * PI / ((double)size - 1.0);
cgsum = 0.0;
igsum = 0.0;
for (i = 0; i < size; i++)
{
arg1 = arg0 * (double)i;
a->window[i] = 0.21557895 - 0.41663158 * cos(arg1) + 0.277263158 * cos(2.0 * arg1) - 0.083578947 * cos(3.0 * arg1) + 0.006947368 * cos (4.0 * arg1);
cgsum += a->window[i];
igsum += a->window[i] * a->window[i];
}
a->inv_coherent_gain = (double)size / cgsum;
for (i = 0; i < size; i++)
a->window[i] *= a->inv_coherent_gain;
break;
}
case 4: // hamming window
{
arg0 = 2.0 * PI / ((double)size - 1.0);
cgsum = 0.0;
igsum = 0.0;
for (i = 0; i < size; i++)
{
a->window[i] = (0.54 - 0.46 * cos((double)i * arg0));
cgsum += a->window[i];
igsum += a->window[i] * a->window[i];
}
a->inv_coherent_gain = (double)size / cgsum;
for (i = 0; i < size; i++)
a->window[i] *= a->inv_coherent_gain;
break;
}
case 5: // Kaiser window
{ arg0 = bessi0(PiAlpha);
arg1 = (double)(size - 1);
cgsum = 0.0;
igsum = 0.0;
for (i = 0; i < size; ++i)
{
a->window[i] = bessi0(PiAlpha * sqrt(1.0 - pow(2.0 * (double)i / arg1 - 1.0, 2))) / arg0;
cgsum += a->window[i];
igsum += a->window[i] * a->window[i];
}
a->inv_coherent_gain = (double)size / cgsum;
for (i = 0; i < size; i++)
a->window[i] *= a->inv_coherent_gain;
break;
}
case 6: // Blackman-Harris window (7-term)
default:
{
arg0 = 2.0 * PI / ((double)size - 1.0);
cgsum = 0.0;
igsum = 0.0;
for (i = 0; i < size; ++i)
{
arg1 = cos (arg0 * (double)i);
a->window[i] = + 6.3964424114390378e-02
+ arg1 * ( - 2.3993864599352804e-01
+ arg1 * ( + 3.5015956323820469e-01
+ arg1 * ( - 2.4774111897080783e-01
+ arg1 * ( + 8.5438256055858031e-02
+ arg1 * ( - 1.2320203369293225e-02
+ arg1 * ( + 4.3778825791773474e-04 ))))));
cgsum += a->window[i];
igsum += a->window[i] * a->window[i];
}
a->inv_coherent_gain = (double)size / cgsum;
for (i = 0; i < size; i++)
a->window[i] *= a->inv_coherent_gain;
break;
}
}
a->inherent_power_gain = igsum / (double)size;
a->inv_enb = 1.0 / (a->inherent_power_gain * a->inv_coherent_gain * a->inv_coherent_gain);
// print_window_gain ("windows.txt", type, a->inv_coherent_gain, a->inherent_power_gain);
}
// spur elimination, REAL input data
void eliminate(int disp, int ss, int LO)
{
DP a = pdisp[disp];
int i, k, begin, end, ilim;
double mag;
if (ss == a->begin_ss)
begin = a->fscL + a->clip;
else
begin = a->clip;
if (ss == a->end_ss)
end = a->out_size - 1 -a->clip - a->fscH;
else
end = a->out_size - 1 - a->clip;
ilim = a->out_size - 1;
if (a->flip[LO])
for (i = ilim - begin, k = 0; i > ilim - end; i--, k++)
{
mag = (a->fft_out[ss][LO])[i][0] * (a->fft_out[ss][LO])[i][0] + (a->fft_out[ss][LO])[i][1] * (a->fft_out[ss][LO])[i][1];
if ((a->spec_flag[ss] == 0) || (mag < (a->result[ss])[k]))
(a->result[ss])[k] = mag;
}
else
for (i = begin, k = 0; i < end; i++, k++)
{
mag = (a->fft_out[ss][LO])[i][0] * (a->fft_out[ss][LO])[i][0] + (a->fft_out[ss][LO])[i][1] * (a->fft_out[ss][LO])[i][1];
if ((a->spec_flag[ss] == 0) || (mag < (a->result[ss])[k]))
(a->result[ss])[k] = mag;
}
a->ss_bins[ss] = k;
}
// spur elimination, COMPLEX input data
void Celiminate(int disp, int ss, int LO)
{
DP a = pdisp[disp];
int i, k, begin0, end0, begin1, end1, ilim;
double mag;
if (ss == a->begin_ss)
{
begin0 = a->out_size / 2 + 1 + a->clip + a->fscL;
if (begin0 > a->out_size)
begin1 = begin0 - a->out_size;
else
begin1 = 0;
}
else
{
begin0 = a->out_size / 2 + 1 + a->clip;
begin1 = 0;
}
if (ss == a->end_ss)
{
end1 = a->out_size / 2 - a->clip - a->fscH;
if (end1 < 0)
end0 = a->out_size + end1;
else
end0 = a->out_size;
}
else
{
end0 = a->out_size;
end1 = a->out_size / 2 - a->clip;
}
ilim = a->out_size - 1;
if (a->flip[LO])
{
for (i = ilim - begin0, k = 0; i > ilim - end0; i--, k++)
{
mag = (a->fft_out[ss][LO])[i][0] * (a->fft_out[ss][LO])[i][0] + (a->fft_out[ss][LO])[i][1] * (a->fft_out[ss][LO])[i][1];
if ((a->spec_flag[ss] == 0) || (mag < (a->result[ss])[k]))
(a->result[ss])[k] = mag;
}
for (i = ilim - begin1; i > ilim - end1; i--, k++)
{
mag = (a->fft_out[ss][LO])[i][0] * (a->fft_out[ss][LO])[i][0] + (a->fft_out[ss][LO])[i][1] * (a->fft_out[ss][LO])[i][1];
if ((a->spec_flag[ss] == 0) || (mag < (a->result[ss])[k]))
(a->result[ss])[k] = mag;
}
}
else
{
for (i = begin0, k = 0; i < end0; i++, k++)
{
mag = (a->fft_out[ss][LO])[i][0] * (a->fft_out[ss][LO])[i][0] + (a->fft_out[ss][LO])[i][1] * (a->fft_out[ss][LO])[i][1];
if ((a->spec_flag[ss] == 0) || (mag < (a->result[ss])[k]))
(a->result[ss])[k] = mag;
}
for (i = begin1; i < end1; i++, k++)
{
mag = (a->fft_out[ss][LO])[i][0] * (a->fft_out[ss][LO])[i][0] + (a->fft_out[ss][LO])[i][1] * (a->fft_out[ss][LO])[i][1];
if ((a->spec_flag[ss] == 0) || (mag < (a->result[ss])[k]))
(a->result[ss])[k] = mag;
}
}
a->ss_bins[ss] = k;
}
void detector ( int det_type, // detector type
int m, // number of bins
int num_pixels, // number of output pixels
double pix_per_bin, // pixels per bin
double bin_per_pix, // bins per pixel
double* bins, // input buffer
double* pixels, // output buffer
double inv_enb, // inverse equivalent noise bandwidth
double fsclipL,
double fsclipH,
double det_offset
)
{
int i, imin, ilim;
int pix_count = 0;
int rose, fell, next_pix_count, bcount, last_pix_count;
double prev_maxi, mini, maxi, psum;
if (pix_per_bin <= 1.0)
{
if (fsclipL == floor(fsclipL)) imin = 0;
else imin = 1;
if (fsclipH == floor(fsclipH)) ilim = m;
else ilim = m - 1;
switch (det_type)
{
case 0: // positive peak
for (i = 0; i < num_pixels; i++)
pixels[i] = - 1.0e300;
for (i = imin; i < ilim; i++)
{
pix_count = (int)(det_offset + (double)i * pix_per_bin);
if (pix_count >= num_pixels) pix_count = num_pixels - 1;
if (bins[i] > pixels[pix_count])
pixels[pix_count] = bins[i];
}
break;
case 1: // rosenfell
rose = 0;
fell = 0;
mini = + 1.0e300;
maxi = - 1.0e300;
prev_maxi = - 1.0e300;
for (i = imin; i < ilim; i++) // for each FFT bin
{
// determine the pixel number that this FFT bin goes into
pix_count = (int)(det_offset + (double)i * pix_per_bin);
if (pix_count >= num_pixels) pix_count = num_pixels - 1;
// determine the pixel number for the NEXT FFT bin
next_pix_count = (int)((double)(i + 1) * pix_per_bin);
// update the minimum and maximum of the set of bins within the pixel
if (bins[i] < mini) mini = bins[i];
if (bins[i] > maxi) maxi = bins[i];
// if the next bin is also within the pixel && there is a next bin,
// compare its value with the current bin and update rose and fell
if (next_pix_count == pix_count && i < ilim - 1)
{
// NOTE: when next_pix_count != pix_count, rose and fell do not get updated;
// that's OK because we do NOT need to know if there's a rise or fall across bins
if (bins[i + 1] > bins[i]) rose = 1;
if (bins[i + 1] < bins[i]) fell = 1;
}
// if the next bin is NOT within the pixel || there is no next bin, finalize the pixel
// value and reset parameters
else
{
if (rose && fell)
if (pix_count & 1) // odd pixel
pixels[pix_count] = max (prev_maxi, maxi);
else // even pixel
pixels[pix_count] = mini;
else
pixels[pix_count] = maxi;
rose = 0;
fell = 0;
prev_maxi = maxi;
mini = + 1.0e300;
maxi = - 1.0e300;
}
}
break;
case 2: // average - adjusted for window's equivalent noise bandwidth
psum = 0.0;
bcount = 0;
for (i = imin; i < ilim; i++)
{
last_pix_count = pix_count;
pix_count = (int)(det_offset + (double)i * pix_per_bin);
if (pix_count >= num_pixels) pix_count = num_pixels - 1;
if (pix_count == last_pix_count)
{
psum += bins[i];
bcount++;
}
else
{
pixels[last_pix_count] = psum / (double)bcount * inv_enb;
psum = bins[i];
bcount = 1;
}
if (i == ilim - 1)
{
pixels[pix_count] = psum / (double)bcount * inv_enb;
}
}
break;
case 3: // sample - adjusted for window's equivalent noise bandwidth
bcount = 0;
for (i = imin; i < ilim; i++)
{
last_pix_count = pix_count;
pix_count = (int)(det_offset + (double)i * pix_per_bin);
if (pix_count >= num_pixels) pix_count = num_pixels - 1;
if (pix_count == last_pix_count)
{
bcount++;
}
else
{
pixels[last_pix_count] = bins[i - bcount / 2 - 1] * inv_enb;
bcount = 1;
}
if (i == ilim - 1)
{
pixels[pix_count] = bins[i - bcount / 2] * inv_enb;
}
}
break;
case 4: // rms
psum = 0.0;
bcount = 0;
for (i = imin; i < ilim; i++)
{
last_pix_count = pix_count;
pix_count = (int)(det_offset + (double)i * pix_per_bin);
if (pix_count >= num_pixels) pix_count = num_pixels - 1;
if (pix_count == last_pix_count)
{
psum += bins[i] * bins[i];
bcount++;
}
else
{
pixels[last_pix_count] = sqrt (psum / (double)bcount) * inv_enb;
psum = bins[i] * bins[i];
bcount = 1;
}
if (i == ilim - 1)
{
pixels[pix_count] = sqrt (psum / (double)bcount) * inv_enb;
}
}
break;
}
}
else
{
double frac;
double pix_pos = fsclipL - floor(fsclipL);
int ampl_comp = (det_type == 2) || (det_type == 3) || (det_type == 4);
for (i = 1; i < m; i++)
{
while (pix_pos < ((double)i + 1.0e-06) && pix_count < num_pixels)
{
frac = pix_pos - (double)(i - 1);
pixels[pix_count] = bins[i - 1] * (1.0 - frac) + bins[i] * frac;
if (ampl_comp) pixels[pix_count] *= inv_enb;
pix_count++;
pix_pos += bin_per_pix;
}
}
}
}
void avenger ( int av_mode, // averaging mode
int num_pixels, // number of pixels
int* avail_frames, // number of available frames for window averaging
int num_average, // number of frames to average within a window
int* av_in_idx, // in index for av_buff
int* av_out_idx, // out index for av_buff
double av_backmult, // multiplier for recursive averaging
double scale, // scale factor
double* t_pixels, // input buffer
double* av_sum, // history buffer for averaging
double** av_buff, // frame buffer for window averaging
double* cd, // correction factor buffer
int norm, // if TRUE, normalize to one Hz bandwidth
double norm_oneHz, // normalization factor to add
dOUTREAL* pixels // output buffer
)
{
int i;
double factor;
switch (av_mode)
{
case -1: // peak-hold
{
for (i = 0; i < num_pixels; i++)
{
if (t_pixels[i] > av_sum[i])
av_sum[i] = t_pixels[i];
pixels[i] = (dOUTREAL)(10.0 * mlog10(scale * cd[i] * av_sum[i] + 1.0e-60));
}
break;
}
case 0: // no averaging
default:
{
for (i = 0; i < num_pixels; i++)
pixels[i] = (dOUTREAL)(10.0 * mlog10(scale * cd[i] * t_pixels[i] + 1.0e-60));
break;
}
case 1: // weighted averaging of linear data
{
double onem_avb = 1.0 - av_backmult;
for (i = 0; i < num_pixels; i++)
{
av_sum[i] = av_backmult * av_sum[i] + onem_avb * t_pixels[i];
pixels[i] = (dOUTREAL)(10.0 * mlog10(scale * cd[i] * av_sum[i] + 1.0e-60));
}
break;
}
case 2: // window averaging of linear data
{
if (*avail_frames < num_average)
{
factor = scale / (double)++(*avail_frames);
for (i = 0; i < num_pixels; i++)
{
av_sum[i] += t_pixels[i];
av_buff[*av_in_idx][i] = t_pixels[i];
pixels[i] = (dOUTREAL)(10.0 * mlog10(cd[i] * av_sum[i] * factor + 1.0e-60));
}
}
else
{
factor = scale / (double)(*avail_frames);
for (i = 0; i < num_pixels; i++)
{
av_sum[i] += t_pixels[i] - (av_buff[*av_out_idx])[i];
av_buff[*av_in_idx][i] = t_pixels[i];
pixels[i] = (dOUTREAL)(10.0 * mlog10(cd[i] * av_sum[i] * factor + 1.0e-60));
}
if (++(*av_out_idx) == dMAX_AVERAGE)
*av_out_idx = 0;
}
if (++(*av_in_idx) == dMAX_AVERAGE)
*av_in_idx = 0;
break;
}
case 3: // weighted averaging of log data - looks nice, not accurate for time-varying signals
{
double onem_avb = 1.0 - av_backmult;
for (i = 0; i < num_pixels; i++)
{
av_sum[i] = av_backmult * av_sum[i] + onem_avb * (10.0 * mlog10(scale * cd[i] * t_pixels[i] + 1e-60));
pixels[i] = (dOUTREAL)av_sum[i];
}
break;
}
}
if (norm)
for (i = 0; i < num_pixels; i++)
pixels[i] += (dOUTREAL)norm_oneHz;
}
void stitch(int disp)
{
DP a = pdisp[disp];
int i, j, k, n, m;
double* ptr;
// stitch
m = 0;
ptr = a->pre_av_out;
for (n = a->begin_ss; n <= a->end_ss; n++)
{
memcpy(ptr, a->result[n], a->ss_bins[n] * sizeof(double));
ptr += a->ss_bins[n];
m += a->ss_bins[n];
}
for (i = 0; i < a->num_pixout; i++) // for each output
{
EnterCriticalSection(&a->ResampleSection);
// if a detection of the same 'det_type' has already been done, use that result
j = i - 1;
k = i;
while (j >= 0)
{
if (a->det_type[i] == a->det_type[j])
k = j;
j--;
}
if (k == i)
// detect
detector (a->det_type[i], m, a->num_pixels, a->pix_per_bin, a->bin_per_pix, a->pre_av_out,
a->t_pixels[i], a->inv_enb, a->fsclipL, a->fsclipH, a->det_offset);
else
memcpy (a->t_pixels[i], a->t_pixels[k], a->num_pixels * sizeof (double));
// average & convert to dBm
avenger (a->av_mode[i], a->num_pixels, &a->avail_frames[i], a->num_average[i], &a->av_in_idx[i], &a->av_out_idx[i],
a->av_backmult[i], a->scale, a->t_pixels[i], a->av_sum[i], a->av_buff[i], a->cd, a->normalize[i], a->norm_oneHz,
a->pixels[i][a->w_pix_buff[i]]);
LeaveCriticalSection(&a->ResampleSection);
EnterCriticalSection(&a->PB_ControlsSection[i]);
a->last_pix_buff[i] = a->w_pix_buff[i];
while ((a->w_pix_buff[i] = (a->w_pix_buff[i] + 1) % dNUM_PIXEL_BUFFS) == a->r_pix_buff[i]);
LeaveCriticalSection(&a->PB_ControlsSection[i]);
InterlockedBitTestAndSet(&(a->pb_ready[i][a->last_pix_buff[i]]), 0);
}
}
DWORD WINAPI spectra (void *pargs)
{
int i, j;
int disp = ((int)(uintptr_t)pargs) >> 12;
int ss = (((int)(uintptr_t)pargs) >> 4) & 255;
int LO = ((int)(uintptr_t)pargs) & 15;
DP a = pdisp[disp];
if (a->stop)
{
InterlockedDecrement(a->pnum_threads);
return 0;
}
if ((ss >= a->begin_ss) && (ss <= a->end_ss))
{
for (i = 0; i < a->size; i++)
{
(a->fft_in[ss][LO])[i] = a->window[i] * (double)((a->I_samples[ss][LO])[a->IQO_idx[ss][LO]]);
if(++a->IQO_idx[ss][LO] >= a->bsize)
a->IQO_idx[ss][LO] -= a->bsize;
}
if (a->stop)
{
InterlockedDecrement(a->pnum_threads);
return 0;
}
fftw_execute (a->plan[ss][LO]);
}
if (a->stop)
{
InterlockedDecrement(a->pnum_threads);
return 0;
}
EnterCriticalSection(&(a->EliminateSection[ss]));
if ((ss >= a->begin_ss) && (ss <= a->end_ss))
eliminate(disp, ss, LO);
a->spec_flag[ss] |= 1 << LO;
if (a->spec_flag[ss] == ((1 << a->num_fft) - 1))
{
a->spec_flag[ss] = 0;
LeaveCriticalSection (&(a->EliminateSection[ss]));
EnterCriticalSection (&a->StitchSection);
a->stitch_flag |= ((uint64_t)1) << ss;
if (a->stitch_flag == ((((uint64_t)1) << a->num_stitch) - 1))
{
a->stitch_flag = 0;
LeaveCriticalSection(&a->StitchSection);
for (j = 0; j < dMAX_STITCH; j++)
for (i = 0; i < dMAX_NUM_FFT; i++)
InterlockedBitTestAndReset(&(a->input_busy[j][i]), 0);
stitch(disp);
}
else
LeaveCriticalSection(&a->StitchSection);
}
else
LeaveCriticalSection (&(a->EliminateSection[ss]));
InterlockedDecrement(a->pnum_threads);
return 1;
}
/********************************************************************************************************
* *
* BEGIN CODE TO GET MAX FFT_BIN WITHIN A FREQ RANGE *
* *
********************************************************************************************************/
// Call in XCreateAnalyzer(...)
// This gets initialized for each 'disp' that is set up.
void Init_DetectMaxBin(int disp)
{
DP a = pdisp[disp];
InitializeCriticalSectionAndSpinCount(&a->cs_dmb, 2500);
a->dmb_run = 0;
a->dmb_disp = disp;
a->dmb_ss = 0;
a->dmb_LO = 0;
a->dmb_rate = 48000.0;
a->dmb_fLow = 0.0;
a->dmb_fHigh = 0.0;
a->dmb_begin0 = 0;
a->dmb_end0 = -1;
a->dmb_begin1 = 0;
a->dmb_end1 = -1;
a->dmb_max_dB = -400.0;
}
// Call in DestroyAnalyzer(...)
void Destroy_DetectMaxBin(int disp)
{
DP a = pdisp[disp];
DeleteCriticalSection(&a->cs_dmb);
}
// Called from SetupDetectMaxBin(...) AND anytime 'size' changes, e.g., in SetAnalyzer(...)
void calc_dmb(int disp, int size)
{
DP a = pdisp[disp];
double bin_spacing, min_freq, max_freq;
bin_spacing = a->dmb_rate / size;
min_freq = -a->dmb_rate * (double)((size / 2 - 1)) / (double)size;
max_freq = -min_freq;
if (a->dmb_fLow < min_freq) a->dmb_fLow = min_freq;
if (a->dmb_fLow > max_freq) a->dmb_fLow = max_freq;
if (a->dmb_fHigh < min_freq) a->dmb_fHigh = min_freq;
if (a->dmb_fHigh > max_freq) a->dmb_fHigh = max_freq;
if (a->dmb_fLow < 0.0 && a->dmb_fHigh < 0.0)
{
a->dmb_begin0 = size - (int)ceil(-a->dmb_fLow / bin_spacing);
// Example: 1024, 48K, -2700 => 966
// Example: 1024, 48K, -23953.125==min_freq => 513
// NOTE: since dmb_fLow < 0.0, begin0 will never be greater than a->size - 1
a->dmb_end0 = size - (int)ceil(-a->dmb_fHigh / bin_spacing);
// Example: 1024, 48K, -300 => 1017
a->dmb_begin1 = 0;
a->dmb_end1 = -1;
}
else if (a->dmb_fLow >= 0.0 && a->dmb_fHigh >= 0.0)
{
a->dmb_begin0 = (int)round(a->dmb_fLow / bin_spacing);
// Example: 1024, 48K, +300 => 6
a->dmb_end0 = (int)round(a->dmb_fHigh / bin_spacing);
// Example: 1024, 48K, +2700 => 58
// Example: 1024, 48K, 23953.125==max_freq => 511
a->dmb_begin1 = 0;
a->dmb_end1 = -1;
}
else // fLow < 0.0 && fHigh >= 0.0
{
a->dmb_begin0 = size - (int)ceil(-a->dmb_fLow / bin_spacing);
a->dmb_end0 = size - 1;
a->dmb_begin1 = 0;
a->dmb_end1 = (int)round(a->dmb_fHigh / bin_spacing);
}
if (a->dmb_tau > 0.0 && a->dmb_frame_rate > 0)
a->dmb_decay = exp(-1.0 / (a->dmb_tau * (double)a->dmb_frame_rate));
else
a->dmb_decay = 0.0;
a->dmb_max_dB = -400.0;
}
// Call from console, for each 'disp' for which this function is desired.
// Call for initial setup and anytime one of the parameters changes.
// run: Set to '1' if this is being used; cycles can be saved by setting to
// '0' if this is not currently in use.
// disp: Display identifier number.
// ss: Set to '0' for Thetis use.
// LO: Set to '0' for Thetis use.
// rate: Sample_rate of display data, e.g., '192000.0'.
// fLow: Lowest frequency of frequency range to evaluate, referenced to center_frequency
// to which DDC is tuned. For example, for LSB, not using CTUN, this might be
// '-3000.0'.
// fHigh: Highest frequency of frequency range to evaluate, referenced to center_frequency
// to which DDC is tuned. Example, LSB: '-300.0'.
// tau: Metering time constant in seconds. E.g., '0.5'.
// frame_rate: Display frame_rate currently in use. E.g., '60'.
PORT
void SetupDetectMaxBin(int run, int disp, int ss, int LO, double rate,
double fLow, double fHigh, double tau, int frame_rate)
{
// We only allow setup for one (ss,LO) pair at a time for each 'disp'.
DP a = pdisp[disp];
a->dmb_run = run;
a->dmb_disp = disp;
a->dmb_ss = ss;
a->dmb_LO = LO;
a->dmb_rate = rate;
a->dmb_fLow = fLow;
a->dmb_fHigh = fHigh;
a->dmb_tau = tau;
a->dmb_frame_rate = frame_rate;
calc_dmb(a->dmb_disp, a->size);
}
// Call this function in 'Cspectra(...)', after the FFT.
void DetectMaxBin(int disp, int ss, int LO)
{
DP a = pdisp[disp];
int i;
double mag, dmb_max;
double dmb_max_dB;
// If 'run' is set and the FFT Output is from the correct disp, ss, LO ...
if (a->dmb_run && disp == a->dmb_disp && ss == a->dmb_ss && LO == a->dmb_LO)
{
fftw_complex* fft_out = a->fft_out[ss][LO];
dmb_max = 1.0e-60;
EnterCriticalSection(&a->cs_dmb);
for (i = a->dmb_begin0; i <= a->dmb_end0; i++)
{
mag = fft_out[i][0] * fft_out[i][0] + fft_out[i][1] * fft_out[i][1];
if (mag > dmb_max) dmb_max = mag;
}
for (i = a->dmb_begin1; i <= a->dmb_end1; i++)
{
mag = fft_out[i][0] * fft_out[i][0] + fft_out[i][1] * fft_out[i][1];
if (mag > dmb_max) dmb_max = mag;
}
a->dmb_max_dB -= fabs((1.0 - a->dmb_decay) * a->dmb_max_dB);
dmb_max_dB = 10.0 * mlog10(a->scale * dmb_max);
if (dmb_max_dB > a->dmb_max_dB) a->dmb_max_dB = dmb_max_dB;
LeaveCriticalSection(&a->cs_dmb);
// for test only.
// printf("Max Bin = %.5e\n", a->dmb_max_dB);
}
}
// Call from console, for each 'disp' for which this function is desired.
// Always returns the value from the most recent display frame.
PORT
double GetDetectMaxBin(int disp)
{
DP a = pdisp[disp];
double dmb_max_dB;
EnterCriticalSection(&a->cs_dmb);
dmb_max_dB = a->dmb_max_dB;
LeaveCriticalSection(&a->cs_dmb);
return dmb_max_dB;
}
/********************************************************************************************************
* *
* END CODE TO GET MAXIMUM FFT_BIN WITHIN A FREQ RANGE *
* *
********************************************************************************************************/
DWORD WINAPI Cspectra (void *pargs)
{
int i, j;
int disp = ((int)(uintptr_t)pargs) >> 12;
int ss = (((int)(uintptr_t)pargs) >> 4) & 255;
int LO = ((int)(uintptr_t)pargs) & 15;
DP a = pdisp[disp];
int trans_size = a->size * sizeof(double);
if (a->stop)
{
InterlockedDecrement(a->pnum_threads);
return 0;
}
if ((ss >= a->begin_ss) && (ss <= a->end_ss))
{
for (i = 0; i < a->size; i++)
{
(a->Cfft_in[ss][LO])[i][0] = a->window[i] * (double)((a->I_samples[ss][LO])[a->IQO_idx[ss][LO]]);
(a->Cfft_in[ss][LO])[i][1] = a->window[i] * (double)((a->Q_samples[ss][LO])[a->IQO_idx[ss][LO]]);
if(++a->IQO_idx[ss][LO] >= a->bsize)
a->IQO_idx[ss][LO] -= a->bsize;
}
if (a->stop)
{
InterlockedDecrement(a->pnum_threads);
return 0;
}
fftw_execute (a->Cplan[ss][LO]);
// Detect value of Max FFT Bin in a freq range
DetectMaxBin(disp, ss, LO);
//
}
if (a->stop)
{
InterlockedDecrement(a->pnum_threads);
return 0;
}
if (InterlockedBitTestAndReset(&(a->snap[ss][LO]), 0))
{
memcpy((char *)(a->snap_buff[ss][LO]), (char *)(a->fft_out[ss][LO]) + trans_size, trans_size);
memcpy((char *)(a->snap_buff[ss][LO]) + trans_size, (char *)(a->fft_out[ss][LO]), trans_size);
SetEvent(a->hSnapEvent[ss][LO]);
}
EnterCriticalSection(&(a->EliminateSection[ss]));
if ((ss >= a->begin_ss) && (ss <= a->end_ss))
Celiminate(disp, ss, LO);
a->spec_flag[ss] |= 1 << LO;
if (a->spec_flag[ss] == ((1 << a->num_fft) - 1))
{
a->spec_flag[ss] = 0;
LeaveCriticalSection (&(a->EliminateSection[ss]));
EnterCriticalSection (&a->StitchSection);
a->stitch_flag |= ((uint64_t)1) << ss;
if (a->stitch_flag == ((((uint64_t)1) << a->num_stitch) - 1))
{
a->stitch_flag = 0;
LeaveCriticalSection(&a->StitchSection);
for (j = 0; j < dMAX_STITCH; j++)
for (i = 0; i < dMAX_NUM_FFT; i++)
InterlockedBitTestAndReset(&(a->input_busy[j][i]), 0);
stitch(disp);
}
else
LeaveCriticalSection(&a->StitchSection);
}
else
LeaveCriticalSection (&(a->EliminateSection[ss]));
InterlockedDecrement(a->pnum_threads);
return 1;
}
void interpolate(int disp, int set, double fmin, double fmax, int num_pixels)
{
DP a = pdisp[disp];
int i;
double f;
int n = a->n_freqs[set];
int k;
int kmin = 0;
int kmax = n - 1;
int kdelta;
double dx;
double mag;
for (i = 0; i < num_pixels; i++)
{
f = fmin + (double)i * (fmax - fmin) / (double)(num_pixels - 1);
if (f < (a->freqs[set])[0])
k = 0;
else if (f > (a->freqs[set])[n - 1])
k = n - 2;
else
{
kdelta = 1;
while (f < (a->freqs[set])[kmin])
{
kmin = max(0, kmin - kdelta);
kdelta += kdelta;
}
while (f > (a->freqs[set])[kmax])
{
kmax = min(n - 1, kmax + kdelta);
kdelta += kdelta;
}
k = (kmin + kmax) / 2;
while ((kmax - kmin) > 1)
{
k = (kmin + kmax) / 2;
if (f > (a->freqs[set])[k])
kmin = k;
else
kmax = k--;
}
}
dx = f - (a->freqs[set])[k];
mag = (((a->ac3[set][0])[k] * dx + (a->ac2[set][0])[k]) * dx + (a->ac1[set][0])[k]) * dx + (a->ac0[set][0])[k];
a->cd[i] = mag * mag;
}
return;
}
int build_interpolants(int disp, int set, int n, int m, double *x, double (*y)[dMAX_M])
{
DP a = pdisp[disp];
double dx[dMAX_N];
double idx[dMAX_N];
double dmain[dMAX_N];
double dsub[dMAX_N];
double dsup[dMAX_N];
double d[dMAX_N][dMAX_M];
double S[dMAX_N][dMAX_M];
double b[dMAX_N];
double v[dMAX_N][dMAX_M];
double tmp;
int i, j;
if (n < 3) return -1;
for (i = 0; i < n - 1; i++)
{
dx[i] = x[i + 1] - x[i];
if (dx[i] < 1e-30)
return -1;
idx[i] = 1.0 / dx[i];
}
for (i = 1; i <= n - 2; i++)
{
if (i == 1)
{
dsub[i] = 0.0;
dmain[i] = 3.0 * dx[i - 1] + 2.0 * dx[i];
dsup[i] = dx[i];
}
else if (i == (n - 2))
{
dsub[i] = dx[i - 1];
dmain[i] = 2.0 * dx[i - 1] + 3.0 * dx[i];
dsup[i] = 0.0;
}
else
{
dsub[i] = dx[i - 1];
dmain[i] = 2.0 * (dx[i - 1] + dx[i]);
dsup[i] = dx[i];
}
for (j = 0; j < m; j++)
d[i][j] = 6.0 * ((y[i + 1][j] - y[i][j]) * idx[i] - (y[i][j] - y[i - 1][j]) * idx[i - 1]);
}
b[1] = dmain[1];
for (j = 0; j < m; j++)
v[1][j] = d[1][j];
for (i = 2; i <= n - 2; i++)
{
tmp = dsub[i] / b[i - 1];
b[i] = dmain[i] - tmp * dsup[i - 1];
for (j = 0; j < m; j++)
v[i][j] = d[i][j] - tmp * v[i - 1][j];
}
for (j = 0; j < m; j++)
S[n - 2][j] = v[n - 2][j] / b[n - 2];
for (i = n - 3; i >= 1; i--)
for (j = 0; j < m; j++)
S[i][j] = (v[i][j] - dsup[i] * S[i + 1][j]) / b[i];
for (j = 0; j < m; j++)
{
S[0][j] = S[1][j];
S[n - 1][j] = S[n - 2][j];
}
for (i = 0; i < n - 1; i++)
for (j = 0; j < m; j++)
{
(a->ac3[set][j])[i] = (S[i + 1][j] - S[i][j]) / (6.0 * dx[i]);
(a->ac2[set][j])[i] = 0.5 * S[i][j];
(a->ac1[set][j])[i] = (y[i + 1][j] - y[i][j]) * idx[i] - (2.0 * dx[i] * S[i][j] + dx[i] * S[i + 1][j]) / 6.0;
(a->ac0[set][j])[i] = y[i][j];
}
return 0;
}
void __cdecl sendbuf(void *arg)
{
DP a = pdisp[(int)(uintptr_t)arg];
while(!a->end_dispatcher)
{
for (a->ss = 0; a->ss < a->num_stitch; a->ss++)
for (a->LO = 0; a->LO < a->num_fft; a->LO++)
{
if (!_InterlockedAnd(&(a->input_busy[a->ss][a->LO]), 1) && _InterlockedAnd(&(a->buff_ready[a->ss][a->LO]), 1))
{
InterlockedBitTestAndSet(&(a->input_busy[a->ss][a->LO]), 0);
a->IQO_idx[a->ss][a->LO] = a->IQout_index[a->ss][a->LO];
InterlockedIncrement(a->pnum_threads);
if (a->type == 0)
QueueUserWorkItem(spectra, (void *)(((uintptr_t)arg << 12) + (a->ss << 4) + a->LO), 0);
else
QueueUserWorkItem(Cspectra, (void *)(((uintptr_t)arg << 12) + (a->ss << 4) + a->LO), 0);
if((a->IQout_index[a->ss][a->LO] += a->incr) >= a->bsize)
a->IQout_index[a->ss][a->LO] -= a->bsize;
EnterCriticalSection(&(a->BufferControlSection[a->ss][a->LO]));
if ((a->have_samples[a->ss][a->LO] -= a->incr) < a->size)
InterlockedBitTestAndReset(&(a->buff_ready[a->ss][a->LO]), 0);
LeaveCriticalSection(&(a->BufferControlSection[a->ss][a->LO]));
}
}
//
// Block until a Spectrum*() call announces new samples. This used to be
// Sleep(1), i.e. 1000 wakeups per second spent re-reading the same flags
// -- 0.67% of a core even with no data arriving at all.
//
// Whoever sets end_dispatcher also signals the semaphore, so the wait
// below always has a way out.
//
WaitForSingleObject(a->Sem_BuffReady, INFINITE);
}
InterlockedBitTestAndReset(&a->dispatcher, 0);
_endthread();
}
void CalcBandwidthNormalization (DP a)
{
double bin_width;
if (a->sample_rate <= 0 || a->size <= 0) { a->norm_oneHz = 0.0; return; }
bin_width = (double)a->sample_rate / (double)a->size;
if (bin_width <= 0.0) { a->norm_oneHz = 0.0; return; }
a->norm_oneHz = 10.0 * mlog10 (1.0 / bin_width);
}
PORT
void ResetPixelBuffers(int disp)
{
//[2.10.2]MW0LGE reset all the pixel and average buffers
DP a = pdisp[disp];
int i, j, k;
EnterCriticalSection(&a->SetAnalyzerSection);
EnterCriticalSection(&a->ResampleSection);
for (i = 0; i < dMAX_PIXOUTS; i++)
{
for (j = 0; j < dMAX_PIXELS; j++)
a->t_pixels[i][j] = 0.0;
for (j = 0; j < dMAX_AVERAGE; j++)
for (k = 0; k < dMAX_PIXELS; k++)
a->av_buff[i][j][k] = 0.0f;
switch (a->av_mode[i])
{
case 1:
for (j = 0; j < dMAX_PIXELS; j++)
a->av_sum[i][j] = 1.0e-12;
break;
case 2:
//done below
//a->avail_frames[i] = 0;
//a->av_in_idx[i] = 0;
//a->av_out_idx[i] = 0;
break;
case 3:
for (j = 0; j < dMAX_PIXELS; j++)
a->av_sum[i][j] = -160.0;
break;
default:
memset((void*)a->av_sum[i], 0, sizeof(double) * dMAX_PIXELS);
break;
}
a->avail_frames[i] = 0;
a->av_in_idx[i] = 0;
a->av_out_idx[i] = 0;
EnterCriticalSection(&a->PB_ControlsSection[i]);
a->w_pix_buff[i] = 0;
a->r_pix_buff[i] = 0;
a->last_pix_buff[i] = 0;
for (j = 0; j < dNUM_PIXEL_BUFFS; j++)
a->pb_ready[i][j] = 0;
LeaveCriticalSection(&a->PB_ControlsSection[i]);
}
memset((void*)a->pre_av_out, 0, sizeof(double) * a->max_size * a->max_stitch);
LeaveCriticalSection(&a->ResampleSection);
EnterCriticalSection(&a->StitchSection);
for (i = 0; i < dMAX_STITCH; i++)
for (j = 0; j < dMAX_NUM_FFT; j++)
a->input_busy[i][j] = 0;
for (i = 0; i < dMAX_STITCH; i++)
a->spec_flag[i] = 0;
a->stitch_flag = 0;
a->ss = 0;
a->LO = 0;
for (i = 0; i < dMAX_STITCH; i++)
for (j = 0; j < dMAX_NUM_FFT; j++)
{
EnterCriticalSection(&(a->BufferControlSection[i][j]));
a->buff_ready[i][j] = 0;
a->have_samples[i][j] = 0;
a->IQin_index[i][j] = 0;
a->IQout_index[i][j] = 0;
LeaveCriticalSection(&(a->BufferControlSection[i][j]));
}
LeaveCriticalSection(&a->StitchSection);
LeaveCriticalSection(&a->SetAnalyzerSection);
}
PORT
void SetAnalyzer ( int disp, // display identifier
int n_pixout, // pixel output identifier
int n_fft, // number of LO frequencies = number of ffts used in elimination
int typ, // 0 for real input data (I only); 1 for complex input data (I & Q)
int *flp, // vector with one elt for each LO frequency, 1 if high-side LO, 0 otherwise
int sz, // size of the fft, i.e., number of input samples
int bf_sz, // number of samples transferred for each OpenBuffer()/CloseBuffer()
int win_type, // integer specifying which window function to use
double pi, // PiAlpha parameter for Kaiser window
int ovrlp, // number of samples each fft (other than the first) is to re-use from the previous
int clp, // number of fft output bins to be clipped from EACH side of each sub-span
double fscLin, // number of bins to clip from low end of entire span
double fscHin, // number of bins to clip from high end of entire span
int n_pix, // number of pixel values to return. may be either <= or > number of bins
int n_stch, // number of sub-spans to concatenate to form a complete span
int calset, // identifier of which set of calibration data to use
double fmin, // frequency at first pixel value
double fmax, // frequency at last pixel value
int max_w
)
{
WDSP_FPE_GUARD;
DP a = pdisp[disp];
int i, j;
EnterCriticalSection(&a->SetAnalyzerSection);
a->end_dispatcher = 1;
// wake the dispatcher out of its blocking wait so it can observe the flag;
// it does not take SetAnalyzerSection, so holding it here is safe
ReleaseSemaphore(a->Sem_BuffReady, 1, 0);
while (InterlockedAnd(&a->dispatcher, 1))
Sleep(1);
a->stop = 1;
while (_InterlockedAnd(a->pnum_threads, 1023))
Sleep(1);
a->num_pixout = n_pixout;
a->num_fft = n_fft;
a->type = typ;
a->buff_size = bf_sz;
for (i = 0; i < a->num_fft; i++)
a->flip[i] = *(flp + i);
a->overlap = ovrlp;
a->clip = clp;
a->fsclipL = fscLin;
a->fsclipH = fscHin;
a->num_stitch = n_stch;
if (sz != a->size)
{
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
{
if (a->plan[i][j]) fftw_destroy_plan (a->plan[i][j]);
if (a->Cplan[i][j]) fftw_destroy_plan (a->Cplan[i][j]);
a->plan[i][j] = fftw_plan_dft_r2c_1d(sz, a->fft_in[i][j], a->fft_out[i][j], FFTW_PATIENT);
a->Cplan[i][j] = fftw_plan_dft_1d(sz, a->Cfft_in[i][j], a->fft_out[i][j], FFTW_FORWARD, FFTW_PATIENT);
}
// Setup DetectMaxBin for a 'size' change.
calc_dmb(disp, sz);
//
}
if ((sz != a->size) || (win_type != a->window_type) || (pi != a->PiAlpha))
new_window(disp, win_type, sz, pi);
a->size = sz;
a->window_type = win_type;
a->PiAlpha = pi;
a->max_writeahead = max_w;
CalcBandwidthNormalization (a);
if (((fmin != a->f_min) || (fmax != a->f_max)) && ((fmin == 0.0) && (fmax == 0.0)))
for (i = 0; i < dMAX_PIXELS; i++)
a->cd[i] = 1.0;
if (((fmax != 0.0) || (fmin != 0.0)) && ((n_pix != a->num_pixels) || (fmin != a->f_min) || (fmax != a->f_max) || (calset != a->cal_set) || a->cal_changed))
interpolate(disp, calset, fmin, fmax, n_pix);
a->incr = a->size - a->overlap;
a->num_pixels = n_pix;
a->f_min = fmin;
a->f_max = fmax;
a->cal_set = calset;
a->cal_changed = 0;
if (a->type == 0)
{
a->out_size = a->size / 2 + 1;
a->scale = 4.0 / ((double)a->size * (double)a->size);
}
else
{
a->out_size = a->size;
a->scale = 1.0 / ((double)a->size * (double)a->size);
}
a->begin_ss = 0;
a->end_ss = a->num_stitch - 1;
a->fscL = (int)a->fsclipL;
a->fscH = (int)a->fsclipH;
while (a->fscL >= (a->out_size - 1 - 2 * a->clip))
{
a->fscL -= a->out_size - 1 - 2 * a->clip;
a->ss_bins[a->begin_ss] = 0;
a->begin_ss++;
}
while (a->fscH >= (a->out_size - 1 - 2 * a->clip))
{
a->fscH -= a->out_size - 1 - 2 * a->clip;
a->ss_bins[a->end_ss] = 0;
a->end_ss--;
}
a->pix_per_bin = (double)a->num_pixels / ((double)(a->num_stitch * (a->out_size - 1 - 2 * a->clip)) - a->fsclipL - a->fsclipH - 1.0);
a->det_offset = -a->pix_per_bin * (a->fsclipL - floor(a->fsclipL));
a->bin_per_pix = ((double)(a->num_stitch * (a->out_size - 1 - 2 * a->clip)) - 1.0 - a->fsclipL - a->fsclipH) / ((double)a->num_pixels - 1.0);
for (i = 0; i < dMAX_STITCH; i++)
for (j = 0; j < dMAX_NUM_FFT; j++)
a->input_busy[i][j] = 0;
for (i = 0; i < dMAX_STITCH; i++)
a->spec_flag[i] = 0;
a->stitch_flag = 0;
for (i = 0; i < dMAX_PIXOUTS; i++)
{
a->w_pix_buff[i] = 0;
a->r_pix_buff[i] = 0;
a->last_pix_buff[i] = 0;
for (j = 0; j < dNUM_PIXEL_BUFFS; j++)
a->pb_ready[i][j] = 0;
}
a->ss = 0;
a->LO = 0;
for (i = 0; i < dMAX_STITCH; i++)
for (j = 0; j < dMAX_NUM_FFT; j++)
{
a->buff_ready[i][j] = 0;
a->have_samples[i][j] = 0;
a->IQin_index[i][j] = 0;
a->IQout_index[i][j] = 0;
}
a->stop = 0;
a->end_dispatcher = 0;
LeaveCriticalSection(&a->SetAnalyzerSection);
}
PORT
void XCreateAnalyzer( int disp,
int *success,
int m_size,
int m_num_fft,
int m_stitch,
char *app_data_path
)
{
WDSP_FPE_GUARD;
int i, j;
DP a = (DP) malloc0 (sizeof(dp));
pdisp[disp] = a;
a->max_size = m_size;
a->max_num_fft = m_num_fft;
a->max_stitch = m_stitch;
a->pnum_threads = (LONG*) malloc0 (sizeof (LONG));
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
{
a->hSnapEvent[i][j] = CreateEvent(NULL, FALSE, FALSE, TEXT("snap"));
a->snap[i][j] = 0;
}
a->Sem_BuffReady = CreateSemaphore(0, 0, 1000, 0);
InitializeCriticalSectionAndSpinCount(&a->ResampleSection, 0);
InitializeCriticalSectionAndSpinCount(&a->SetAnalyzerSection, 0);
InitializeCriticalSectionAndSpinCount(&a->StitchSection, 0);
for (i = 0; i < dMAX_PIXOUTS; i++)
InitializeCriticalSectionAndSpinCount(&a->PB_ControlsSection[i], 0);
for (i = 0; i < dMAX_STITCH; i++)
{
InitializeCriticalSectionAndSpinCount(&(a->EliminateSection[i]), 0);
for (j = 0; j < dMAX_NUM_FFT; j++)
InitializeCriticalSectionAndSpinCount(&(a->BufferControlSection[i][j]), 0);
}
a->window = (double*) malloc0 (sizeof(double) * a->max_size);
for (i = 0; i < a->max_stitch; i++)
{
a->result[i] = (double*) malloc0 (sizeof(double) * a->max_size);
}
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
{
a->plan[i][j] = 0;
a->Cplan[i][j] = 0;
a->fft_in[i][j] = (double*) malloc0 (sizeof(double) * a->max_size);
a->Cfft_in[i][j] = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * a->max_size);
a->fft_out[i][j] = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * a->max_size);
}
a->pre_av_out = (double*) malloc0 (sizeof(double) * a->max_size * a->max_stitch);
for (i = 0; i < dMAX_PIXOUTS; i++)
{
a->det_type[i] = 0;
a->av_mode[i] = 0;
a->av_sum[i] = (double*) malloc0 (sizeof(double) * dMAX_PIXELS);
for (j = 0; j < dMAX_AVERAGE; j++)
a->av_buff[i][j] = (double*) malloc0 (sizeof(double) * dMAX_PIXELS);
a->t_pixels[i] = (double*) malloc0 (sizeof(double) * dMAX_PIXELS);
for (j = 0; j < dNUM_PIXEL_BUFFS; j++)
a->pixels[i][j] = (dOUTREAL*) malloc0 (sizeof(dOUTREAL) * dMAX_PIXELS);
}
a->cd = (double*) malloc0 (sizeof(double) * dMAX_PIXELS);
for (j = 0; j < dMAX_PIXELS; j++)
a->cd[j] = 1.0;
for (i = 0; i < dMAX_CAL_SETS; i++)
{
a->freqs[i] = (double*) malloc0 (sizeof(double) * dMAX_N);
for (j = 0; j < dMAX_M; j++)
{
a->ac3[i][j] = (double*) malloc0 (sizeof(double) * dMAX_N);
a->ac2[i][j] = (double*) malloc0 (sizeof(double) * dMAX_N);
a->ac1[i][j] = (double*) malloc0 (sizeof(double) * dMAX_N);
a->ac0[i][j] = (double*) malloc0 (sizeof(double) * dMAX_N);
}
}
a->size = -1;
a->window_type = -1;
a->num_pixels = -1;
a->cal_set = -1;
a->f_min = -1.0;
a->f_max = -1.0;
a->bsize = a->max_size * dSAMP_BUFF_MULT;
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
{
a->I_samples[i][j] = (dINREAL*) malloc0 (sizeof(dINREAL) * a->bsize);
a->Q_samples[i][j] = (dINREAL*) malloc0 (sizeof(dINREAL) * a->bsize);
}
// Initialize DetectMaxBin functionality
Init_DetectMaxBin(disp);
// for test only.
// SetupDetectMaxBin(1, 0, 0, 0, 192000.0, -3000.0, -300.0, 0.5, 60);
//
*success = 0;
}
PORT
void DestroyAnalyzer(int disp)
{
DP a = pdisp[disp];
int i, j;
a->end_dispatcher = 1;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0);
while (InterlockedAnd(&a->dispatcher, 1))
Sleep(1);
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
{
_aligned_free (a->I_samples[i][j]);
_aligned_free (a->Q_samples[i][j]);
}
for (i = 0; i < dMAX_CAL_SETS; i++)
{
_aligned_free (a->freqs[i]);
for (j = 0; j < dMAX_M; j++)
{
_aligned_free (a->ac3[i][j]);
_aligned_free (a->ac2[i][j]);
_aligned_free (a->ac1[i][j]);
_aligned_free (a->ac0[i][j]);
}
}
_aligned_free (a->cd);
for (i = 0; i < dMAX_PIXOUTS; i++)
{
for (j = 0; j < dNUM_PIXEL_BUFFS; j++)
_aligned_free (a->pixels[i][j]);
_aligned_free (a->t_pixels[i]);
for (j = 0; j < dMAX_AVERAGE; j++)
_aligned_free (a->av_buff[i][j]);
_aligned_free (a->av_sum[i]);
}
_aligned_free (a->pre_av_out);
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
{
fftw_destroy_plan (a->plan[i][j]);
fftw_destroy_plan (a->Cplan[i][j]);
fftw_free (a->Cfft_in[i][j]);
_aligned_free (a->fft_in[i][j]);
fftw_free (a->fft_out[i][j]);
}
for (i = 0; i < a->max_stitch; i++)
_aligned_free (a->result[i]);
_aligned_free (a->window);
for (i = 0; i < dMAX_STITCH; i++)
{
DeleteCriticalSection(&(a->EliminateSection[i]));
for (j = 0; j < dMAX_NUM_FFT; j++)
DeleteCriticalSection(&(a->BufferControlSection[i][j]));
}
for (i = 0; i < dMAX_PIXOUTS; i++)
DeleteCriticalSection(&a->PB_ControlsSection[i]);
DeleteCriticalSection(&a->StitchSection);
DeleteCriticalSection(&a->SetAnalyzerSection);
DeleteCriticalSection(&a->ResampleSection);
for (i = 0; i < a->max_stitch; i++)
for (j = 0; j < a->max_num_fft; j++)
CloseHandle(a->hSnapEvent[i][j]);
CloseHandle(a->Sem_BuffReady);
_aligned_free ((void *) a->pnum_threads);
// Destroy DetectMaxBin functionality.
Destroy_DetectMaxBin(disp);
//
_aligned_free (a);
}
PORT
void GetPixels ( int disp,
int pixout,
dOUTREAL *pix, //if new pixel values avail, copies to pix and sets flag = 1
int *flag //else, returns 0 (try again later)
)
{
DP a = pdisp[disp];
EnterCriticalSection(&a->PB_ControlsSection[pixout]);
a->r_pix_buff[pixout] = a->last_pix_buff[pixout];
LeaveCriticalSection(&a->PB_ControlsSection[pixout]);
if (_InterlockedAnd(&(a->pb_ready[pixout][a->r_pix_buff[pixout]]), 1))
{
memcpy (pix, a->pixels[pixout][a->r_pix_buff[pixout]], a->num_pixels * sizeof(dOUTREAL));
*flag = 1;
InterlockedBitTestAndReset(&(a->pb_ready[pixout][a->r_pix_buff[pixout]]), 0);
}
else
*flag = 0;
}
PORT
void SnapSpectrum( int disp,
int ss,
int LO,
double *snap_buff)
{
DP a = pdisp[disp];
a->snap_buff[ss][LO] = snap_buff;
InterlockedBitTestAndSet(&(a->snap[ss][LO]), 0);
WaitForSingleObject(a->hSnapEvent[ss][LO], INFINITE);
}
PORT
void SnapSpectrumTimeout( int disp,
int ss,
int LO,
double* snap_buff,
DWORD timeout,
int* flag)
{
DP a = pdisp[disp];
a->snap_buff[ss][LO] = snap_buff;
ResetEvent(a->hSnapEvent[ss][LO]);
InterlockedBitTestAndSet(&(a->snap[ss][LO]), 0);
if (!WaitForSingleObject(a->hSnapEvent[ss][LO], timeout))
*flag = 1;
else
{
InterlockedBitTestAndReset(&(a->snap[ss][LO]), 0);
*flag = 0;
}
}
int calcompare (const void * a, const void * b)
{
if (*(double*)a < *(double*)b)
return -1;
else if (*(double*)a == *(double*)b)
return 0;
else
return 1;
}
PORT
void SetCalibration ( int disp,
int set_num, //identifier for this calibration data set
int n_points, //number of calibration points in the set
double (*cal)[dMAX_M+1] //pointer to the calibration table, first
) // column is frequency, add'l columns are
// data for variables being calibrated
{
DP a = pdisp[disp];
int i, j;
int k = 0;
double y [dMAX_N][dMAX_M];
qsort (cal, n_points, (dMAX_M+1) * sizeof(double), calcompare);
for (i = 0; i < n_points; i++)
{
if ((i == n_points - 1) || (cal[i][0] != cal[i + 1][0]))
{
(a->freqs[set_num])[k] = cal[i][0];
for (j = 0; j < dMAX_M; j++)
y[k][j] = cal[i][j + 1];
k++;
}
}
a->n_freqs[set_num] = k;
build_interpolants(disp, set_num, a->n_freqs[set_num], dMAX_M, a->freqs[set_num], y);
a->cal_changed = 1;
}
PORT
void OpenBuffer(int disp, int ss, int LO, void **Ipointer, void **Qpointer)
{
DP a = pdisp[disp];
EnterCriticalSection(&a->SetAnalyzerSection);
*Ipointer = &((a->I_samples[ss][LO])[a->IQin_index[ss][LO]]);
*Qpointer = &((a->Q_samples[ss][LO])[a->IQin_index[ss][LO]]);
LeaveCriticalSection(&a->SetAnalyzerSection);
}
PORT
void CloseBuffer(int disp, int ss, int LO)
{
DP a = pdisp[disp];
EnterCriticalSection(&a->SetAnalyzerSection);
EnterCriticalSection(&(a->BufferControlSection[ss][LO]));
if (a->have_samples[ss][LO] > a->max_writeahead)
{
//if we're receiving samples too much faster than we're consuming them, skip some
if ((a->IQout_index[ss][LO] += a->have_samples[ss][LO] - a->max_writeahead) >= a->bsize)
a->IQout_index[ss][LO] -= a->bsize;
a->have_samples[ss][LO] = a->max_writeahead;
}
if ((a->have_samples[ss][LO] += a->buff_size) >= a->size)
InterlockedBitTestAndSet(&(a->buff_ready[ss][LO]), 0);
LeaveCriticalSection(&(a->BufferControlSection[ss][LO]));
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet (&a->dispatcher, 0);
LeaveCriticalSection(&a->SetAnalyzerSection);
_beginthread(sendbuf, 0, (void *)(uintptr_t)disp);
}
else
LeaveCriticalSection(&a->SetAnalyzerSection);
}
PORT
void Spectrum(int disp, int ss, int LO, dINREAL* pI, dINREAL* pQ)
{
dINREAL *Ipointer;
dINREAL *Qpointer;
DP a = pdisp[disp];
EnterCriticalSection(&a->SetAnalyzerSection);
Ipointer = &((a->I_samples[ss][LO])[a->IQin_index[ss][LO]]);
Qpointer = &((a->Q_samples[ss][LO])[a->IQin_index[ss][LO]]);
LeaveCriticalSection(&a->SetAnalyzerSection);
memcpy(Ipointer, pI, a->buff_size * sizeof(dINREAL));
memcpy(Qpointer, pQ, a->buff_size * sizeof(dINREAL));
EnterCriticalSection(&a->SetAnalyzerSection);
EnterCriticalSection(&(a->BufferControlSection[ss][LO]));
if (a->have_samples[ss][LO] > a->max_writeahead)
{
//if we're receiving samples too much faster than we're consuming them, skip some
if ((a->IQout_index[ss][LO] += a->have_samples[ss][LO] - a->max_writeahead) >= a->bsize)
a->IQout_index[ss][LO] -= a->bsize;
a->have_samples[ss][LO] = a->max_writeahead;
}
if ((a->have_samples[ss][LO] += a->buff_size) >= a->size)
InterlockedBitTestAndSet(&(a->buff_ready[ss][LO]), 0);
LeaveCriticalSection(&(a->BufferControlSection[ss][LO]));
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet(&a->dispatcher, 0);
LeaveCriticalSection(&a->SetAnalyzerSection);
_beginthread(sendbuf, 0, (void *)(uintptr_t)disp);
}
else
LeaveCriticalSection(&a->SetAnalyzerSection);
}
PORT
void Spectrum2(int run, int disp, int ss, int LO, dINREAL* pbuff)
{
if (run)
{
int i;
dINREAL *Ipointer;
dINREAL *Qpointer;
DP a = pdisp[disp];
EnterCriticalSection(&a->SetAnalyzerSection);
Ipointer = &((a->I_samples[ss][LO])[a->IQin_index[ss][LO]]);
Qpointer = &((a->Q_samples[ss][LO])[a->IQin_index[ss][LO]]);
LeaveCriticalSection(&a->SetAnalyzerSection);
for (i = 0; i < a->buff_size; i++)
{
Ipointer[i] = pbuff[2 * i + 1];
Qpointer[i] = pbuff[2 * i + 0];
}
EnterCriticalSection(&a->SetAnalyzerSection);
EnterCriticalSection(&(a->BufferControlSection[ss][LO]));
if (a->have_samples[ss][LO] > a->max_writeahead)
{
//if we're receiving samples too much faster than we're consuming them, skip some
if ((a->IQout_index[ss][LO] += a->have_samples[ss][LO] - a->max_writeahead) >= a->bsize)
a->IQout_index[ss][LO] -= a->bsize;
a->have_samples[ss][LO] = a->max_writeahead;
}
if ((a->have_samples[ss][LO] += a->buff_size) >= a->size)
InterlockedBitTestAndSet(&(a->buff_ready[ss][LO]), 0);
LeaveCriticalSection(&(a->BufferControlSection[ss][LO]));
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet(&a->dispatcher, 0);
LeaveCriticalSection(&a->SetAnalyzerSection);
_beginthread(sendbuf, 0, (void *)(uintptr_t)disp);
}
else
LeaveCriticalSection(&a->SetAnalyzerSection);
}
}
PORT
void Spectrum0(int run, int disp, int ss, int LO, double* pbuff)
{
if (run)
{
int i;
dINREAL *Ipointer;
dINREAL *Qpointer;
DP a = pdisp[disp];
EnterCriticalSection(&a->SetAnalyzerSection);
Ipointer = &((a->I_samples[ss][LO])[a->IQin_index[ss][LO]]);
Qpointer = &((a->Q_samples[ss][LO])[a->IQin_index[ss][LO]]);
LeaveCriticalSection(&a->SetAnalyzerSection);
for (i = 0; i < a->buff_size; i++)
{
Ipointer[i] = (dINREAL)pbuff[2 * i + 1];
Qpointer[i] = (dINREAL)pbuff[2 * i + 0];
}
EnterCriticalSection(&a->SetAnalyzerSection);
EnterCriticalSection(&(a->BufferControlSection[ss][LO]));
if (a->have_samples[ss][LO] > a->max_writeahead)
{
//if we're receiving samples too much faster than we're consuming them, skip some
if ((a->IQout_index[ss][LO] += a->have_samples[ss][LO] - a->max_writeahead) >= a->bsize)
a->IQout_index[ss][LO] -= a->bsize;
a->have_samples[ss][LO] = a->max_writeahead;
}
if ((a->have_samples[ss][LO] += a->buff_size) >= a->size)
InterlockedBitTestAndSet(&(a->buff_ready[ss][LO]), 0);
LeaveCriticalSection(&(a->BufferControlSection[ss][LO]));
if((a->IQin_index[ss][LO] += a->buff_size) >= a->bsize) //REQUIRES buff_size IS A SUB-MULTIPLE OF SIZE OF INPUT SAMPLE BUFFS!
a->IQin_index[ss][LO] = 0;
ReleaseSemaphore(a->Sem_BuffReady, 1, 0); // new samples: let the dispatcher run
if (!InterlockedAnd(&a->dispatcher, 1))
{
InterlockedBitTestAndSet(&a->dispatcher, 0);
LeaveCriticalSection(&a->SetAnalyzerSection);
_beginthread(sendbuf, 0, (void *)(uintptr_t)disp);
}
else
LeaveCriticalSection(&a->SetAnalyzerSection);
}
}
PORT
void SetDisplayDetectorMode (int disp, int pixout, int mode)
{
DP a = pdisp[disp];
if (a->det_type[pixout] != mode)
{
EnterCriticalSection (&a->ResampleSection);
a->det_type[pixout] = mode;
LeaveCriticalSection (&a->ResampleSection);
}
}
PORT
void SetDisplayAverageMode (int disp, int pixout, int mode)
{
int i;
DP a = pdisp[disp];
if (a->av_mode[pixout] != mode)
{
EnterCriticalSection (&a->ResampleSection);
a->av_mode[pixout] = mode;
switch (mode)
{
case 1:
for (i = 0; i < dMAX_PIXELS; i++)
a->av_sum[pixout][i] = 1.0e-12;
break;
case 2:
a->avail_frames[pixout] = 0;
a->av_in_idx[pixout] = 0;
a->av_out_idx[pixout] = 0;
break;
case 3:
for (i = 0; i < dMAX_PIXELS; i++)
a->av_sum[pixout][i] = -160.0;
break;
default:
memset ((void *)a->av_sum[pixout], 0, sizeof(double) * dMAX_PIXELS);
break;
}
LeaveCriticalSection (&a->ResampleSection);
}
}
PORT
void SetDisplayNumAverage (int disp, int pixout, int num)
{
DP a = pdisp[disp];
if (a->num_average[pixout] != num)
{
EnterCriticalSection (&a->ResampleSection);
a->num_average[pixout] = num;
a->avail_frames[pixout] = 0;
a->av_in_idx[pixout] = 0;
a->av_out_idx[pixout] = 0;
LeaveCriticalSection (&a->ResampleSection);
}
}
PORT
void SetDisplayAvBackmult (int disp, int pixout, double mult)
{
DP a = pdisp[disp];
if (a->av_backmult[pixout] != mult)
{
EnterCriticalSection (&a->ResampleSection);
a->av_backmult[pixout] = mult;
LeaveCriticalSection (&a->ResampleSection);
}
}
PORT
void SetDisplaySampleRate (int disp, int rate)
{
DP a = pdisp[disp];
if (a->sample_rate != rate)
{
EnterCriticalSection (&a->ResampleSection);
a->sample_rate = rate;
CalcBandwidthNormalization (a);
LeaveCriticalSection (&a->ResampleSection);
}
}
PORT
void SetDisplayNormOneHz (int disp, int pixout, int norm)
{
DP a = pdisp[disp];
if (a->normalize[pixout] != norm)
{
EnterCriticalSection (&a->ResampleSection);
a->normalize[pixout] = norm;
LeaveCriticalSection (&a->ResampleSection);
}
}
PORT
double GetDisplayENB (int disp)
{
DP a = pdisp[disp];
double enb;
EnterCriticalSection(&a->SetAnalyzerSection);
enb = 1.0 / a->inv_enb;
LeaveCriticalSection(&a->SetAnalyzerSection);
return enb;
}
Reference in New Issue View Git Blame Copy Permalink