All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines
Brief description of class still missing. More...
#include <HRFK.h>
Public Member Functions | |
| Complex * | crossCorrelationMatrix (int component) |
| Complex * | crossCorrelationMatrix (double angle) |
| HRFK (QList< FKStationSignals * > array) | |
| void | initOperator (Complex *covmat, double dampingFactor, QString *log=0) |
| virtual double | value (double x, double y) const |
| virtual double | value (double x, double y, int index) const |
| ~HRFK () | |
Brief description of class still missing.
Full description of class still missing
| ArrayCore::HRFK::HRFK | ( | QList< FKStationSignals * > | array | ) |
Description of constructor still missing
References ArrayCore::FK::_array, ArrayCore::FKStationSignals::isSelected(), and TRACE.
: FK(array) { TRACE; _selectedStationCount=0; _stationIndexes=new int[_array.count()]; int n=_array.count(); for(int i=0; i<n; i++ ) { FKStationSignals * s=static_cast<FKStationSignals *>(_array.at(i)); if(s->isSelected()) { _stationIndexes[_selectedStationCount]=i; _selectedStationCount++; } } // Allocation of cross corelation matrix called R _Rmatrix=new Complex[ _selectedStationCount * _selectedStationCount ]; // Proper identification of horizontal components if any switch (array.first()->nComponents()) { case 2: _northIndex=0; _eastIndex=1; break; case 3: _northIndex=1; _eastIndex=2; break; default: _northIndex=0; _eastIndex=0; break; } }
| Complex * ArrayCore::HRFK::crossCorrelationMatrix | ( | int | component | ) |
References ArrayCore::FK::_array, ArrayCore::FK::_gaussianPtr, ArrayCore::FK::_iFreqMax, ArrayCore::FK::_iFreqMin, QGpCoreTools::conjugate(), and ArrayCore::FKStationSignals::getSignalSpectrum().
Referenced by HRFKLoopTask::initGridValues(), and HRFKTimeWindows::setComponent().
{
Complex * covmat=new Complex[ _selectedStationCount * _selectedStationCount ];
// Cache for signals
Complex * sig=new Complex[ _selectedStationCount ];
/*
Filling in the upper part of the matrix and the diagonal elements with the cross products
Summation over all frequency samples
*/
Complex tmp, specRow, specCol;
QList<FKStationSignals *>::iterator itrow, itcol;
for(register int iFreq=_iFreqMin;iFreq <= _iFreqMax;iFreq++ ) {
// Fill in the cache with signals
for(register int s=0;s < _selectedStationCount;s++ ) {
FKStationSignals * stat=_array.at(_stationIndexes[s] );
sig[s]=stat->getSignalSpectrum(component, iFreq);
}
for(register int row=0;row < _selectedStationCount;row++ ) {
for(register int col=row;col < _selectedStationCount;col++ ) {
tmp=sig[col];
tmp=conjugate(tmp);
tmp *= sig[row];
tmp *= _gaussianPtr[iFreq];
covmat[ col * _selectedStationCount + row ] += tmp;
}
}
}
delete [] sig;
return covmat;
}
| Complex * ArrayCore::HRFK::crossCorrelationMatrix | ( | double | angle | ) |
angle is counted from North clockwize
References ArrayCore::FK::_array, ArrayCore::FK::_gaussianPtr, ArrayCore::FK::_iFreqMax, ArrayCore::FK::_iFreqMin, QGpCoreTools::conjugate(), QGpCoreTools::Angle::cos(), ArrayCore::FKStationSignals::getSignalSpectrum(), QGpCoreTools::Angle::setDegreeAzimuth(), and QGpCoreTools::Angle::sin().
{
Angle a;
a.setDegreeAzimuth(angle);
Complex * covmat=new Complex[ _selectedStationCount * _selectedStationCount ];
// Cache for rotated signals
Complex * sig=new Complex[ _selectedStationCount ];
/*
Filling in the upper part of the matrix and the diagonal elements with the cross products
Summation over all frequency samples
*/
Complex tmp, specRow, specCol;
QList<FKStationSignals *>::iterator itrow, itcol;
for(register int iFreq=_iFreqMin;iFreq <= _iFreqMax;iFreq++ ) {
// Fill in the cache with rotated signals
for(register int s=0;s < _selectedStationCount;s++ ) {
FKStationSignals * stat=_array.at(_stationIndexes[s] );
Complex& c=sig[s];
c=stat->getSignalSpectrum(_eastIndex, iFreq);
c *= a.cos();
tmp=stat->getSignalSpectrum(_northIndex, iFreq);
tmp *= a.sin();
c += tmp;
}
for(register int row=0;row < _selectedStationCount;row++ ) {
for(register int col=row;col < _selectedStationCount;col++ ) {
tmp=sig[col];
tmp=conjugate(tmp);
tmp *= sig[row];
tmp *= _gaussianPtr[iFreq];
covmat[ col * _selectedStationCount + row ] += tmp;
}
}
}
delete [] sig;
return covmat;
}
| void ArrayCore::HRFK::initOperator | ( | Complex * | covmat, |
| double | dampingFactor, | ||
| QString * | log = 0 |
||
| ) |
References ArrayCore::FK::_frequency, ArrayCore::FK::_iFreqMax, ArrayCore::FK::_iFreqMin, QGpCoreTools::conjugate(), QGpCoreTools::Complex::re, QGpCoreTools::sqrt(), QGpCoreTools::tr(), and ArrayCore::zgesvd_().
Referenced by HRFKLoopTask::initGridValues(), and HRFKTimeWindows::setComponent().
{
/*
LAPACK is assuming matrix stored in vectors column by column
For a MxN matrix:
a(1,1) ... a(m,1) a(1,2) ... a(m,2) ... a(1,n) ... a(m,n)
On output for _Rmatrix it is the contrary:
For a MxN matrix:
a(1,1) ... a(1,n) a(2,1) ... a(2,n) ... a(m,1) ... a(m,n)
*/
int nostat2=_selectedStationCount * _selectedStationCount;
//printf("Upper\n");
//printMatrix(covmat, _selectedStationCount, _selectedStationCount);
/*
Computes the auto-power for normalizing, normalizing factor also include the division by the number
of frequency samples
*/
double * scale=new double[ _selectedStationCount ];
for(register int row=0;row < _selectedStationCount;row++ ) {
double s=covmat[ row * _selectedStationCount + row ].re();
if(s==0) { // Flat spectrum for one station
return initOperatorError( "Flat spectrum, check signal", log);
}
scale[ row ]=1.0/sqrt(s);
}
/*
Filling in the lower part of the matrix with conjugates of the upper part and normalizing
*/
for(register int row=0;row < _selectedStationCount;row++ ) {
for(register int col=row;col < _selectedStationCount;col++ ) {
Complex& upper=covmat[ col * _selectedStationCount + row ];
Complex& lower=covmat[ row * _selectedStationCount + col ];
upper *= scale[ row ] * scale[ col ];
lower=conjugate(upper);
}
}
delete [] scale;
//printf("Full:\n");
//printMatrix(covmat, _selectedStationCount, _selectedStationCount);
// Damping
if(dampingFactor>0) {
dampingFactor=1.0-dampingFactor;
for(register int row=0;row < _selectedStationCount;row++ ) {
for(register int col=row;col < _selectedStationCount;col++ ) {
covmat[ row * _selectedStationCount + col ] *= dampingFactor;
}
}
/*
Add white noise to avoid singular matrix
*/
dampingFactor=1.0-dampingFactor;
for(register int el=0; el < _selectedStationCount;el++ ) {
covmat[ el + _selectedStationCount * el ] += dampingFactor;
}
}
//printf("With noise:\n");
//printMatrix(covmat, _selectedStationCount, _selectedStationCount);
/*
need to allocate the matrices, etc.
I'll use now code from CLAPACK!
matrices are allocated in 1 dim.-array!
*/
Complex *a; /* copy of covmat */
Complex *vt; /* nullspace vectors */
Complex *u; /* eigenvectors */
double *s; /* eigenvalues */
Complex *work; /* work space */
double *rwork; /* another workspace */
char jobu, jobvt; /* flags what to do in zgesvd_ */
int lda, ldu, ldvt; /* dimensions of SVD problem */
int lwork, info; /* dimension of workspaces */
int mm, nn;
a=new Complex[ nostat2 ];
u=new Complex[ nostat2 ];
vt=new Complex[ nostat2 ];
work=new Complex[ 50 * _selectedStationCount + 1 ];
rwork=new double[ 5 * _selectedStationCount ];
s=new double [ _selectedStationCount ];
/*
Copy covmat to a
*/
for(register int i=nostat2 - 1;i >= 0;i-- ) {
a[ i ]=covmat[ i ];
}
//printf("a\n");
//printMatrix(a, _selectedStationCount, _selectedStationCount);
/*********************************************/
/* now find the eigenvalues and eigenvectors */
/* driver ZGESVD already sorts the basis! */
/*********************************************/
// CLAPACK is probalby not thread safe... lock it
jobu='A';
jobvt='A';
lda=ldu=ldvt=mm=nn=_selectedStationCount;
lwork=50 * _selectedStationCount + 1;
static QMutex mutex;
mutex.lock();
zgesvd_(&jobu, &jobvt, &mm, &nn, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, rwork, &info);
mutex.unlock();
/* Output is U, SIGMA and VT
A=U * SIGMA * VT
VT being the conjugate-transpose matrix of V
A^-1=V * SIGMA^-1 * UT
UT being the conjugate-transpose matrix of U
*/
//printf("u\n");
//printMatrix(u, _selectedStationCount, _selectedStationCount);
//printf("s\n");
//printMatrix(s, 1, _selectedStationCount
//printf("vt\n");
//printMatrix(vt, _selectedStationCount, _selectedStationCount);*/
// Check if matrix has a true inverse (good rank)
delete [] work;
delete [] rwork;
delete [] covmat;
delete [] a;
double sRef=1e-12 * s[ 0 ];
for(register int row=0;row < _selectedStationCount;row++ ) {
if(s[ row ] < sRef) {
static const char * msg=QT_TR_NOOP( "At %1 Hz, effective rank %2 (%3 frequency samples), "
"singular cross spectral matrix. The effective "
"rank should be at least equal to the number of stations (%4). Either "
"increase the time window lenght or the frequency band width. "
"An empy fk map will be produced." );
delete [] vt;
delete [] u;
delete [] s;
return initOperatorError(tr( msg).arg(_frequency)
.arg(row)
.arg(_iFreqMax-_iFreqMin+1)
.arg(_selectedStationCount), log);
}
}
// allocate a new temporatry matrix for multiplication of S and U
Complex * invSUT=new Complex[ nostat2 ];
for(register int row=0;row < _selectedStationCount;row++ ) {
double si=1./s[ row ];
for(register int col=0;col < _selectedStationCount;col++ ) {
Complex& invSUTElement=invSUT[ col * _selectedStationCount + row ];
Complex& uElement=u[ row * _selectedStationCount + col ];
invSUTElement=uElement;
invSUTElement=conjugate(invSUTElement);
invSUTElement *= si;
}
}
delete [] u;
delete [] s;
//printf("invSUT\n");
//printMatrix(invSUT, _selectedStationCount, _selectedStationCount);
// proceed with V * invSUT
for(register int row=0;row < _selectedStationCount;row++ ) {
for(register int col=0;col < _selectedStationCount;col++ ) {
Complex& RElement=_Rmatrix[ row * _selectedStationCount + col ];
RElement=0;
for(register int k=0;k < _selectedStationCount;k++ ) {
Complex& vtElement=vt[ row * _selectedStationCount + k ];
Complex tmp=vtElement;
tmp=conjugate(tmp);
tmp *= invSUT[ col * _selectedStationCount + k ];
RElement += tmp;
}
}
}
delete [] invSUT;
delete [] vt;
//printf("R\n");
//printMatrix(_Rmatrix, _selectedStationCount, _selectedStationCount);
// Check matrix inversion
/*for(register int row=0;row < _selectedStationCount;row++ ) {
for(register int col=0;col < _selectedStationCount;col++ ) {
Complex sum;
for(int k=0;k<_selectedStationCount;k++) {
Complex tmp=covmat[row * _selectedStationCount + k];
tmp*=_Rmatrix[col * _selectedStationCount + k];
sum += tmp;
}
printf("%12lg+i*%12lg ",sum.re(), sum.im());
}
printf("\n");
}*/
}
| double ArrayCore::HRFK::value | ( | double | x, |
| double | y | ||
| ) | const [virtual] |
Implemement this function to calculate the 2D function at x and y.
Reimplemented from ArrayCore::FK.
References ArrayCore::FK::_array, QGpCoreTools::Complex::abs(), QGpCoreTools::conjugate(), and ArrayCore::FK::maximumK2().
{
double k2=kx * kx + ky * ky;
if(k2 > maximumK2()) return -1;
for(register int i=0;i < _selectedStationCount;i++ ) {
static_cast<FKStationSignals *>(_array.at(_stationIndexes[i])) ->setCurrentShift(kx, ky);
}
Complex sum;
for(register int i=0;i < _selectedStationCount;i++ ) {
Complex shiftStat1=_array.at(_stationIndexes[i] ) ->getShift();
shiftStat1=conjugate(shiftStat1);
for(register int j=0;j < _selectedStationCount;j++ ) {
Complex tmp(shiftStat1);
tmp *= _array.at(_stationIndexes[j] ) ->getShift();
tmp *= _Rmatrix[ j * _selectedStationCount + i ];
sum += tmp;
}
}
return 1.0/sum.abs();
}
| double ArrayCore::HRFK::value | ( | double | x, |
| double | y, | ||
| int | index | ||
| ) | const [virtual] |
Re-mplemement this function to calculate the 2D function for x and y aligned to a grid. x and y correspond to index in the grid: index=iy * nx + ix
If nothing is cached re-implementation of this function is useless.
Reimplemented from ArrayCore::FK.
References ArrayCore::FK::_array, QGpCoreTools::Complex::abs(), QGpCoreTools::conjugate(), and ArrayCore::FK::maximumK2().
{
double k2=kx * kx + ky * ky;
if(k2 > maximumK2()) return -1;
Complex sum;
for(register int i=0;i < _selectedStationCount;i++ ) {
Complex shiftStat1=_array.at(_stationIndexes[i] ) ->getShift(index);
shiftStat1=conjugate(shiftStat1);
for(register int j=0;j < _selectedStationCount;j++ ) {
Complex tmp(shiftStat1);
tmp *= _array.at(_stationIndexes[j] ) ->getShift(index);
tmp *= _Rmatrix[ j * _selectedStationCount + i ];
sum += tmp;
}
}
return 1.0/sum.abs();
}
