cml/libcw
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

cwDsp.h/cpp: Many additions and changes.

Browse Source
This commit is contained in:
kevin 2021-08-15 16:01:27 -04:00
parent a5c7f41ecf
commit fff1af607e
2 changed files with 87 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
cwDsp.cpp
View File

@ -3,6 +3,7 @@
#include "cwCommonImpl.h" #include "cwCommonImpl.h"
#include "cwMem.h" #include "cwMem.h"
#include "cwUtility.h" #include "cwUtility.h"
#include "cwMath.h"
#include "cwVectOps.h" #include "cwVectOps.h"
#include "cwDsp.h" #include "cwDsp.h"
@ -21,7 +22,9 @@ cw::rc_t cw::dsp::fft::test()
real_t srate = xN; real_t srate = xN;
real_t hz = 1; real_t hz = 1;
real_t xV[xN]; real_t xV[xN];
struct ptr_str<real_t>* p = create<real_t>(xN,flags); struct obj_str<real_t>* p = nullptr;
create<real_t>(p,xN,flags);
if(p != nullptr ) if(p != nullptr )
{ {
@ -54,21 +57,21 @@ cw::rc_t cw::dsp::ifft::test()
{ {
typedef double real_t; typedef double real_t;
struct fft::ptr_str<real_t>* ft = nullptr; struct fft::obj_str<real_t>* ft = nullptr;
struct ptr_str<real_t>* ift = nullptr; struct obj_str<real_t>* ift = nullptr;
rc_t rc = kOkRC; rc_t rc = kOkRC;
unsigned xN = 16; unsigned xN = 16;
real_t xV[xN]; real_t xV[xN];
if( (ft = fft::create<real_t>(xN,fft::kToPolarFl)) == nullptr ) if( (rc = fft::create<real_t>(ft,xN,fft::kToPolarFl)) != kOkRC )
{ {
rc = cwLogError(kOpFailRC,"FFT procedure allocation failed."); rc = cwLogError(rc,"FFT procedure allocation failed.");
goto errLabel; goto errLabel;
} }
if((ift = create<real_t>(fft::bin_count(ft))) == nullptr ) if((rc = create<real_t>(ift, fft::bin_count(ft))) != kOkRC )
{ {
rc = cwLogError(kOpFailRC,"IFFT procedure allocation failed."); rc = cwLogError(rc,"IFFT procedure allocation failed.");
goto errLabel; goto errLabel;
} }
@ -78,7 +81,7 @@ cw::rc_t cw::dsp::ifft::test()
fft::exec(ft,xV,xN); fft::exec(ft,xV,xN);
exec(ift, fft::magn(ft), fft::phase(ft) ); exec_polar(ift, fft::magn(ft), fft::phase(ft) );
vop::print( xV, xN, "%f ", "sin " ); vop::print( xV, xN, "%f ", "sin " );
vop::print( magn(ft), fft::bin_count(ft), "%f ", "mag " ); vop::print( magn(ft), fft::bin_count(ft), "%f ", "mag " );
@ -111,7 +114,9 @@ cw::rc_t cw::dsp::convolve::test()
// correct output from apply: real_t zV[] = { 1., 0.5, 0.25,0.1, 1.05,0.5, 0.25,0.1, 1.05,0.5, 0.25,0.1, 0.05,0.0, 0.0 }; // correct output from apply: real_t zV[] = { 1., 0.5, 0.25,0.1, 1.05,0.5, 0.25,0.1, 1.05,0.5, 0.25,0.1, 0.05,0.0, 0.0 };
ptr_str<real_t>* p = create(hV,hN,xN); obj_str<real_t>* p = nullptr;
create(p,hV,hN,xN);
exec(p,xV,xN); exec(p,xV,xN);

114
cwDsp.h
View File

@ -24,17 +24,18 @@ namespace cw
T kaiser_beta_from_sidelobe_reject( const T& sidelobeRejectDb ) T kaiser_beta_from_sidelobe_reject( const T& sidelobeRejectDb )
{ {
double beta; double beta;
T slrDb = sidelobeRejectDb;
if( sidelobeRejectDb < 13.26 ) if( slrDb < 13.26 )
sidelobeRejectDb = 13.26; slrDb = 13.26;
else else
if( sidelobeRejectDb > 120.0) if( slrDb > 120.0)
sidelobeRejectDb = 120.0; slrDb = 120.0;
if( sidelobeRejectDb < 60.0 ) if( slrDb < 60.0 )
beta = (0.76609 * pow(sidelobeRejectDb - 13.26,0.4)) + (0.09834*(sidelobeRejectDb-13.26)); beta = (0.76609 * pow(slrDb - 13.26,0.4)) + (0.09834*(slrDb-13.26));
else else
beta = 0.12438 * (sidelobeRejectDb + 6.3); beta = 0.12438 * (slrDb + 6.3);
return beta; return beta;
} }
@ -49,7 +50,7 @@ namespace cw
{ {
bool zeroFl = false; bool zeroFl = false;
int M = 0; int M = 0;
double den = cmBessel0(beta); // wnd func denominator double den = math::bessel0<T>(beta); // wnd func denominator
int cnt = n; int cnt = n;
int i; int i;
@ -73,7 +74,7 @@ namespace cw
{ {
double v0 = (double)(i - M); double v0 = (double)(i - M);
double num = cmBessel0(beta * sqrt(1.0 - ((v0*v0)/Msqrd))); double num = math::bessel0(beta * sqrt(1.0 - ((v0*v0)/Msqrd)));
dbp[i] = (T)(num/den); dbp[i] = (T)(num/den);
} }
@ -158,9 +159,11 @@ namespace cw
unsigned n = dn/2; unsigned n = dn/2;
T incr = 1.0/n; T incr = 1.0/n;
seq(dbp,n,0,incr); T v0 = 0;
T v1 = 1;
vop::seq(dbp,n,v0,incr);
seq(dbp+n,dn-n,1,-incr); vop::seq(dbp+n,dn-n,v1,-incr);
return dbp; return dbp;
} }
@ -200,7 +203,7 @@ namespace cw
}; };
template< typename T > template< typename T >
struct ptr_str struct obj_str
{ {
unsigned flags; unsigned flags;
@ -222,9 +225,9 @@ namespace cw
}; };
template< typename T > template< typename T >
struct ptr_str<T>* create( unsigned xN, unsigned flags=kToPolarFl ) rc_t create( struct obj_str<T>*& p, unsigned xN, unsigned flags=kToPolarFl )
{ {
struct ptr_str<T>* p = mem::allocZ< ptr_str<T> >(1); p = mem::allocZ< obj_str<T> >(1);
p->flags = flags; p->flags = flags;
p->inN = xN; p->inN = xN;
@ -246,11 +249,11 @@ namespace cw
p->u.dplan = fftw_plan_dft_r2c_1d((int)xN, (double*)p->inV, reinterpret_cast<fftw_complex*>(p->cplxV), FFTW_MEASURE ); p->u.dplan = fftw_plan_dft_r2c_1d((int)xN, (double*)p->inV, reinterpret_cast<fftw_complex*>(p->cplxV), FFTW_MEASURE );
} }
return p; return kOkRC;
} }
template< typename T > template< typename T >
rc_t destroy( struct ptr_str<T>*& p ) rc_t destroy( struct obj_str<T>*& p )
{ {
if( p == nullptr ) if( p == nullptr )
return kOkRC; return kOkRC;
@ -278,7 +281,7 @@ namespace cw
template< typename T > template< typename T >
rc_t exec( struct ptr_str<T>* p, const T* xV, unsigned xN ) rc_t exec( struct obj_str<T>* p, const T* xV, unsigned xN )
{ {
rc_t rc = kOkRC; rc_t rc = kOkRC;
@ -326,13 +329,13 @@ namespace cw
} }
template< typename T > template< typename T >
unsigned bin_count( ptr_str<T>* p ) { return p->binN; } unsigned bin_count( obj_str<T>* p ) { return p->binN; }
template< typename T > template< typename T >
const T* magn( ptr_str<T>* p ) { return p->magV; } const T* magn( obj_str<T>* p ) { return p->magV; }
template< typename T > template< typename T >
const T* phase( ptr_str<T>* p ) { return p->phsV; } const T* phase( obj_str<T>* p ) { return p->phsV; }
rc_t test(); rc_t test();
} }
@ -344,7 +347,7 @@ namespace cw
{ {
template< typename T > template< typename T >
struct ptr_str struct obj_str
{ {
T *outV; T *outV;
std::complex<T> *cplxV; std::complex<T> *cplxV;
@ -360,9 +363,9 @@ namespace cw
}; };
template< typename T > template< typename T >
struct ptr_str<T>* create( unsigned binN ) rc_t create( struct obj_str<T>*& p, unsigned binN )
{ {
struct ptr_str<T>* p = mem::allocZ< ptr_str<T> >(1); p = mem::allocZ< obj_str<T> >(1);
p->binN = binN; p->binN = binN;
p->outN = (binN-1)*2; p->outN = (binN-1)*2;
@ -380,11 +383,11 @@ namespace cw
p->u.dplan = fftw_plan_dft_c2r_1d((int)p->outN, reinterpret_cast<fftw_complex*>(p->cplxV), (double*)p->outV, FFTW_BACKWARD | FFTW_MEASURE ); p->u.dplan = fftw_plan_dft_c2r_1d((int)p->outN, reinterpret_cast<fftw_complex*>(p->cplxV), (double*)p->outV, FFTW_BACKWARD | FFTW_MEASURE );
} }
return p; return kOkRC;;
} }
template< typename T > template< typename T >
rc_t destroy( struct ptr_str<T>*& p ) rc_t destroy( struct obj_str<T>*& p )
{ {
if( p == nullptr ) if( p == nullptr )
return kOkRC; return kOkRC;
@ -410,9 +413,8 @@ namespace cw
template< typename T > template< typename T >
rc_t exec( struct ptr_str<T>* p, const T* magV, const T* phsV ) rc_t exec_polar( struct obj_str<T>* p, const T* magV, const T* phsV )
{ {
rc_t rc = kOkRC; rc_t rc = kOkRC;
if( magV != nullptr && phsV != nullptr ) if( magV != nullptr && phsV != nullptr )
@ -433,10 +435,35 @@ namespace cw
} }
template< typename T > template< typename T >
unsigned out_count( struct ptr_str<T>* p ) { return p->outN; } rc_t exec_rect( struct obj_str<T>* p, const T* iV, const T* rV )
{
rc_t rc = kOkRC;
if( iV != nullptr && rV != nullptr )
{
unsigned i,j;
for(i=0; i<p->binN; ++i)
p->cplxV[i] = std::complex(rV[i], iV[i]);
for(i=p->outN-1,j=1; j<p->binN-1; --i,++j)
p->cplxV[i] = std::complex(rV[j], iV[j]);
if( std::is_same<T,float>::value )
fftwf_execute(p->u.fplan);
else
fftw_execute(p->u.dplan);
}
return rc;
}
template< typename T > template< typename T >
const T* out( struct ptr_str<T>* p ) { return p->outV; } unsigned out_count( struct obj_str<T>* p ) { return p->outN; }
template< typename T >
const T* out( struct obj_str<T>* p ) { return p->outV; }
rc_t test(); rc_t test();
@ -449,10 +476,10 @@ namespace cw
{ {
template< typename T > template< typename T >
struct ptr_str struct obj_str
{ {
struct fft::ptr_str<T>* ft; struct fft::obj_str<T>* ft;
struct ifft::ptr_str<T>* ift; struct ifft::obj_str<T>* ift;
std::complex<T>* hV; std::complex<T>* hV;
unsigned hN; unsigned hN;
@ -464,17 +491,18 @@ namespace cw
}; };
template< typename T > template< typename T >
struct ptr_str<T>* create(const T* hV, unsigned hN, unsigned procSmpN, T hScale=1 ) rc_t create(struct obj_str<T>*& p, const T* hV, unsigned hN, unsigned procSmpN, T hScale=1 )
{ {
struct ptr_str<T>* p = mem::allocZ<struct ptr_str<T>>(1); p = mem::allocZ<struct obj_str<T>>(1);
unsigned cN = nextPowerOfTwo( hN + procSmpN - 1 ); unsigned cN = math::nextPowerOfTwo( hN + procSmpN - 1 );
p->ft = fft::create<T>(cN,0); fft::create<T>(p->ft,cN,0);
unsigned binN = fft::bin_count( p->ft ); unsigned binN = fft::bin_count( p->ft );
p->ift = ifft::create<T>(binN); ifft::create<T>(p->ift, binN);
p->hN = hN; p->hN = hN;
p->hV = mem::allocZ< std::complex<T> >(binN); p->hV = mem::allocZ< std::complex<T> >(binN);
p->outV = mem::allocZ<T>( cN ); p->outV = mem::allocZ<T>( cN );
@ -489,11 +517,11 @@ namespace cw
printf("procN:%i cN:%i hN:%i binN:%i outN:%i\n", procSmpN, cN, hN, binN, p->outN ); printf("procN:%i cN:%i hN:%i binN:%i outN:%i\n", procSmpN, cN, hN, binN, p->outN );
return p; return kOkRC;
} }
template< typename T > template< typename T >
rc_t destroy( struct ptr_str<T>*& pRef ) rc_t destroy( struct obj_str<T>*& pRef )
{ {
if( pRef == nullptr ) if( pRef == nullptr )
return kOkRC; return kOkRC;
@ -507,7 +535,7 @@ namespace cw
} }
template< typename T > template< typename T >
rc_t exec( struct ptr_str<T>* p, const T* xV, unsigned xN ) rc_t exec( struct obj_str<T>* p, const T* xV, unsigned xN )
{ {
// take FT of input signal // take FT of input signal
fft::exec( p->ft, xV, xN ); fft::exec( p->ft, xV, xN );
@ -517,7 +545,7 @@ namespace cw
p->ift->cplxV[i] = p->hV[i] * p->ft->cplxV[i]; p->ift->cplxV[i] = p->hV[i] * p->ft->cplxV[i];
// take the IFFT of the convolved spectrum // take the IFFT of the convolved spectrum
ifft::exec<T>(p->ift,nullptr,nullptr); ifft::exec_polar<T>(p->ift,nullptr,nullptr);
// sum with previous impulse response tail // sum with previous impulse response tail
vop::add( p->outV, (const T*)p->olaV, (const T*)p->ift->outV, p->outN-1 ); vop::add( p->outV, (const T*)p->olaV, (const T*)p->ift->outV, p->outN-1 );
@ -535,9 +563,11 @@ namespace cw
rc_t apply( const T* xV, unsigned xN, const T* hV, unsigned hN, T* yV, unsigned yN, T hScale=1 ) rc_t apply( const T* xV, unsigned xN, const T* hV, unsigned hN, T* yV, unsigned yN, T hScale=1 )
{ {
unsigned procSmpN = std::min(xN,hN); unsigned procSmpN = std::min(xN,hN);
ptr_str<T> *p = create(hV,hN,procSmpN,hScale); obj_str<T> *p = nullptr;
unsigned yi = 0; unsigned yi = 0;
create(p,hV,hN,procSmpN,hScale);
//printf("procSmpN:%i\n",procSmpN); //printf("procSmpN:%i\n",procSmpN);
for(unsigned xi=0; xi<xN && yi<yN; xi+=procSmpN ) for(unsigned xi=0; xi<xN && yi<yN; xi+=procSmpN )