📄 analysis.cpp
字号:
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeBackward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
Matrix<ComplexDouble> IFFT2(Matrix<ComplexDouble>& x)
{
int rows = x.Rows();
int cols = x.Columns();
Matrix<ComplexDouble> y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return y;
}
Matrix<ComplexDouble>& IFFT2I(Matrix<ComplexDouble>& x)
{
int rows = x.Rows();
int cols = x.Columns();
//Matrix<ComplexDouble> y(rows,cols);
int dims[2] = {cols,rows};
int strides[3] = {0,rows,1};
DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Strides
status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Not Inplace
//status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
//if(status != DFTI_NO_ERROR)
//{
// StatusDisplay(status);
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Problem with FFT!");
//}
// Backward scale
status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Commit Dfti descriptor
status = DftiCommitDescriptor(DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
// Compute Forward transform
status = DftiComputeBackward( DescHandle, x.Data());
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Problem with FFT!");
}
status = DftiFreeDescriptor(&DescHandle);
if(status != DFTI_NO_ERROR)
{
StatusDisplay(status);
Utility::Warning("Problem when trying to free the FFT descriptor handle!");
}
return x;
}
//=========================================================================
// Convolution
//=========================================================================
Vector<float> Conv(Vector<float>& input, Vector<float>& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
}
Vector<double> Conv(Vector<double>& input, Vector<double>& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
}
Vector<ComplexFloat> Conv(Vector<ComplexFloat>& input, Vector<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter size is smaller than the input size
if(filter.Length() > input.Length())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter length should not be larger than the input signal!");
}
int fftlength = input.Length()+filter.Length()-1;
if(PowerOfTwo)
{
fftlength = NextPow2(fftlength);
}
Vector<ComplexFloat> temp1;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp1 = Vector<ComplexFloat>(fftlength);
for(int i=0; i<input.Length(); i++)
{
temp1[i] = input[i];
}
}
else if(b == Circular)
{
temp1 = input;
}
else
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Invalid Border type!");
}
if(b == Symmetric)
{
int wid1 = filter.Length()/2;
int wid2 = filter.Length() - wid1 - 1;
if(wid1 != 0)
{
Vector<ComplexFloat> c1 = Reverse( input.Slice(0, wid1-1) );
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector<ComplexFloat> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
temp1.ReadFromVector(c2, input.Length());
}
}
else if(b == Replicate)
{
int wid1 = filter.Length()/2;
int wid2 = filter.Length() - wid1 - 1;
if(wid1 != 0)
{
Vector<ComplexFloat> c1(wid1, input.ElemNC(0));
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector<ComplexFloat> c2(wid2, input.ElemNC(input.Length()-1));
temp1.ReadFromVector(c2, input.Length());
}
}
Vector<ComplexFloat> temp2;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp2 = Vector<ComplexFloat>(fftlength);
for(int i=0; i<filter.Length(); i++)
{
temp2[i] = filter[i];
}
}
else if(b == Circular)
{
temp2 = Vector<ComplexFloat>(input.Length());
for(int i=0; i<filter.Length(); i++)
{
temp2[i] = filter[i];
}
}
else
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Invalid Border type!");
}
Vector<ComplexFloat> o = IFFT( FFT(temp1) * FFT(temp2) );
if(fa == Same && b != Circular)
{
o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
}
else if(fa == Valid)
{
o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
}
else if(fa == Full && PowerOfTwo)
{
o = o.Slice(0, input.Length()+filter.Length()-2);
}
//else if(fa != Full)
//{
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Invalid FilterArea type!");
//}
return o;
}
Vector<ComplexDouble> Conv(Vector<ComplexDouble>& input, Vector<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter size is smaller than the input size
if(filter.Length() > input.Length())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter length should not be larger than the input signal!");
}
int fftlength = input.Length()+filter.Length()-1;
if(PowerOfTwo)
{
fftlength = NextPow2(fftlength);
}
Vector<ComplexDouble> temp1;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp1 = Vector<ComplexDouble>(fftlength);
for(int i=0; i<input.Length(); i++)
{
temp1[i] = input[i];
}
}
else if(b == Circular)
{
temp1 = input;
}
else
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Invalid Border type!");
}
if(b == Symmetric)
{
int wid1 = filter.Length()/2;
int wid2 = filter.Length() - wid1 - 1;
if(wid1 != 0)
{
Vector<ComplexDouble> c1 = Reverse( input.Slice(0, wid1-1) );
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector<ComplexDouble> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
temp1.ReadFromVector(c2, input.Length());
}
}
else if(b == Replicate)
{
int wid1 = filter.Length()/2;
int wid2 = filter.Length() - wid1 - 1;
if(wid1 != 0)
{
Vector<ComplexDouble> c1(wid1, input.ElemNC(0));
temp1.ReadFromVector(c1, temp1.Length()-wid1);
}
if(wid2 != 0)
{
Vector<ComplexDouble> c2(wid2, input.ElemNC(input.Length()-1));
temp1.ReadFromVector(c2, input.Length());
}
}
Vector<ComplexDouble> temp2;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp2 = Vector<ComplexDouble>(fftlength);
for(int i=0; i<filter.Length(); i++)
{
temp2[i] = filter[i];
}
}
else if(b == Circular)
{
temp2 = Vector<ComplexDouble>(input.Length());
for(int i=0; i<filter.Length(); i++)
{
temp2[i] = filter[i];
}
}
else
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Invalid Border type!");
}
Vector<ComplexDouble> o = IFFT( FFT(temp1) * FFT(temp2) );
if(fa == Same && b != Circular)
{
o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
}
else if(fa == Valid)
{
o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
}
else if(fa == Full && PowerOfTwo)
{
o = o.Slice(0, input.Length()+filter.Length()-2);
}
//else if(fa != Full)
//{
// cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
// Utility::RunTimeError("Invalid FilterArea type!");
//}
return o;
}
// =============
// Conv2
// =============
Matrix<float> Conv2(Matrix<float>& input, Matrix<float>& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
}
Matrix<float> Conv2(Vector<float>& filter1, Vector<float>& filter2, Matrix<float>& input, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexFloat(filter1), ToComplexFloat(filter2), ToComplexFloat(input), fa, b, PowerOfTwo) );
}
Matrix<double> Conv2(Matrix<double>& input, Matrix<double>& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
}
Matrix<double> Conv2(Vector<double>& filter1, Vector<double>& filter2, Matrix<double>& input, FilterArea fa, Border b, bool PowerOfTwo)
{
return Real( Conv2(ToComplexDouble(filter1), ToComplexDouble(filter2), ToComplexDouble(input), fa, b, PowerOfTwo) );
}
Matrix<ComplexFloat> Conv2(Matrix<ComplexFloat>& input, Matrix<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo)
{
// Check if the filter size is smaller than the input size
if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
}
// Allocate area for input using the border
int fftrows = input.Rows()+filter.Rows()-1;
int fftcols = input.Columns()+filter.Columns()-1;
if(PowerOfTwo)
{
fftrows = NextPow2(fftrows);
fftcols = NextPow2(fftcols);
}
Matrix<ComplexFloat> temp1;
if(b == ZeroPad || b == Symmetric || b == Replicate)
{
temp1 = Matrix<ComplexFloat>(fftrows, fftcols,0);
for(int j=0; j<input.Columns(); j++)
{
for(int i=0; i<input.Rows(); i++)
{
temp1.ElemNC(i,j) = input.ElemNC(i,j);
}
}
}
else if(b == Circular)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -