Кто-нибудь знает, есть ли бесплатная библиотека с открытым исходным кодом, которая реализовала эти две функции так, как они определены в Matlab?
Спасибо
Исходный код fftshift/ifftshift C/С++
Ответ 1
Или вы можете сделать это сами, набрав type fftshift и перекодировку, что в C++. Это не тот сложный код Matlab.
Изменение: я заметил, что этот ответ несколько раз проголосовал несколько раз и прокомментировал это отрицательно. Я помню время, когда type fftshift был более показательным, чем текущая реализация, но я мог ошибаться. Если бы я мог удалить ответ, я бы, как кажется, больше не релевант.
Вот версия (любезность Октава), которая реализует ее без circshift.
Ответ 2
FFTHIFT/IFFTSHIFT - это причудливый способ сделать CIRCSHIFT. Вы можете проверить, что FFTSHIFT можно переписать как CIRCSHIFT следующим образом. Вы можете определить макросы в C/C++, чтобы выровнять FFTSHIFT до CIRCSHIFT.
A = rand(m, n);
mm = floor(m / 2);
nn = floor(n / 2);
% All three of the following should provide zeros.
circshift(A,[mm, nn]) - fftshift(A)
circshift(A,[mm, 0]) - fftshift(A, 1)
circshift(A,[ 0, nn]) - fftshift(A, 2)
Аналогичные эквиваленты можно найти для IFFTSHIFT.
Циклический сдвиг может быть реализован очень просто со следующим кодом (может быть улучшен с параллельными версиями курса).
template<class ty>
void circshift(ty *out, const ty *in, int xdim, int ydim, int xshift, int yshift)
{
for (int i = 0; i < xdim; i++) {
int ii = (i + xshift) % xdim;
for (int j = 0; j < ydim; j++) {
int jj = (j + yshift) % ydim;
out[ii * ydim + jj] = in[i * ydim + j];
}
}
}
А потом
#define fftshift(out, in, x, y) circshift(out, in, x, y, (x/2), (y/2))
#define ifftshift(out, in, x, y) circshift(out, in, x, y, ((x+1)/2), ((y+1)/2))
Это было сделано немного экспромтом. Потерпите меня, если есть какие-либо проблемы форматирования/синтаксиса.
Ответ 3
Обычно центрирование БПФ выполняется с помощью v (k) = v (k) * (- 1) ** k в временной области. Смещение в частотной области является плохой заменой, поскольку математических причин и вычислительной эффективности. См. Стр. 27 из: http://show.docjava.com/pub/document/jot/v8n6.pdf
Я не уверен, почему документация Matlab делает это так, как они делают, они не дают технической справки.
Ответ 4
Возможно, этот код может помочь. Он выполняет fftshift/ifftshift только для 1D массива в одном буфере. Алгоритм прямого и обратного fftshift для четного числа элементов полностью идентичен.
void swap(complex *v1, complex *v2)
{
complex tmp = *v1;
*v1 = *v2;
*v2 = tmp;
}
void fftshift(complex *data, int count)
{
int k = 0;
int c = (int) floor((float)count/2);
// For odd and for even numbers of element use different algorithm
if (count % 2 == 0)
{
for (k = 0; k < c; k++)
swap(&data[k], &data[k+c]);
}
else
{
complex tmp = data[0];
for (k = 0; k < c; k++)
{
data[k] = data[c + k + 1];
data[c + k + 1] = data[k + 1];
}
data[c] = tmp;
}
}
void ifftshift(complex *data, int count)
{
int k = 0;
int c = (int) floor((float)count/2);
if (count % 2 == 0)
{
for (k = 0; k < c; k++)
swap(&data[k], &data[k+c]);
}
else
{
complex tmp = data[count - 1];
for (k = c-1; k >= 0; k--)
{
data[c + k + 1] = data[k];
data[k] = data[c + k];
}
data[c] = tmp;
}
}
ОБНОВЛЕНО: Также библиотека FFT (включая операции fftshift) для произвольного номера точек можно найти в Optolithium (под OptolithiumC/libs/fourier)
Ответ 5
Я проверил приведенный здесь код и сделал пример проекта для их проверки. Для 1D-кода можно просто использовать std::rotate
template <typename _Real>
static inline
void rotshift(complex<_Real> * complexVector, const size_t count)
{
int center = (int) floor((float)count/2);
if (count % 2 != 0) {
center++;
}
// odd: 012 34 changes to 34 012
std::rotate(complexVector,complexVector + center,complexVector + count);
}
template <typename _Real>
static inline
void irotshift(complex<_Real> * complexVector, const size_t count)
{
int center = (int) floor((float)count/2);
// odd: 01 234 changes to 234 01
std::rotate(complexVector,complexVector +center,complexVector + count);
}
Я предпочитаю использовать std::rotate по коду Алексея из-за его простоты.
Для 2D это усложняется. Для четных чисел это в основном перевернутый левый правый и переворачивается вверх ногами. Для нечетного это алгоритм circshift:
// A =
// 1 2 3
// 4 5 6
// 7 8 9
// fftshift2D(A)
// 9 | 7 8
// --------------
// 3 | 1 2
// 6 | 4 5
// ifftshift2D(A)
// 5 6 | 4
// 8 9 | 7
// --------------
// 2 3 | 1
Здесь я реализовал код circshift с интерфейсом, используя только один массив для ввода и вывода. Для четных чисел требуется только один массив, для нечетных чисел второй массив временно создается и копируется обратно во входной массив. Это приводит к снижению производительности из-за дополнительного времени для копирования массива.
template<class _Real>
static inline
void fftshift2D(complex<_Real> *data, size_t xdim, size_t ydim)
{
size_t xshift = xdim / 2;
size_t yshift = ydim / 2;
if ((xdim*ydim) % 2 != 0) {
// temp output array
std::vector<complex<_Real> > out;
out.resize(xdim * ydim);
for (size_t x = 0; x < xdim; x++) {
size_t outX = (x + xshift) % xdim;
for (size_t y = 0; y < ydim; y++) {
size_t outY = (y + yshift) % ydim;
// row-major order
out[outX + xdim * outY] = data[x + xdim * y];
}
}
// copy out back to data
copy(out.begin(), out.end(), &data[0]);
}
else {
// in and output array are the same,
// values are exchanged using swap
for (size_t x = 0; x < xdim; x++) {
size_t outX = (x + xshift) % xdim;
for (size_t y = 0; y < yshift; y++) {
size_t outY = (y + yshift) % ydim;
// row-major order
swap(data[outX + xdim * outY], data[x + xdim * y]);
}
}
}
}
template<class _Real>
static inline
void ifftshift2D(complex<_Real> *data, size_t xdim, size_t ydim)
{
size_t xshift = xdim / 2;
if (xdim % 2 != 0) {
xshift++;
}
size_t yshift = ydim / 2;
if (ydim % 2 != 0) {
yshift++;
}
if ((xdim*ydim) % 2 != 0) {
// temp output array
std::vector<complex<_Real> > out;
out.resize(xdim * ydim);
for (size_t x = 0; x < xdim; x++) {
size_t outX = (x + xshift) % xdim;
for (size_t y = 0; y < ydim; y++) {
size_t outY = (y + yshift) % ydim;
// row-major order
out[outX + xdim * outY] = data[x + xdim * y];
}
}
// copy out back to data
copy(out.begin(), out.end(), &data[0]);
}
else {
// in and output array are the same,
// values are exchanged using swap
for (size_t x = 0; x < xdim; x++) {
size_t outX = (x + xshift) % xdim;
for (size_t y = 0; y < yshift; y++) {
size_t outY = (y + yshift) % ydim;
// row-major order
swap(data[outX + xdim * outY], data[x + xdim * y]);
}
}
}
}
Ответ 6
Примечание: Есть лучшие ответы, я просто держу это здесь некоторое время... Я не знаю что.
Попробуйте следующее:
template<class T> void ifftShift(T *out, const T* in, size_t nx, size_t ny)
{
const size_t hlen1 = (ny+1)/2;
const size_t hlen2 = ny/2;
const size_t shft1 = ((nx+1)/2)*ny + hlen1;
const size_t shft2 = (nx/2)*ny + hlen2;
const T* src = in;
for(T* tgt = out; tgt < out + shft1 - hlen1; tgt += ny, src += ny) { // (nx+1)/2 times
copy(src, src+hlen1, tgt + shft2); //1->4
copy(src+hlen1, src+ny, tgt+shft2-hlen2); } //2->3
src = in;
for(T* tgt = out; tgt < out + shft2 - hlen2; tgt += ny, src += ny ){ // nx/2 times
copy(src+shft1, src+shft1+hlen2, tgt); //4->1
copy(src+shft1-hlen1, src+shft1, tgt+hlen2); } //3->2
};
Для матриц с четными размерами вы можете сделать это на месте, просто передавая один и тот же указатель в параметры ввода и вывода.
Также обратите внимание, что для 1D-массивов fftshift просто std:: rotate.
Ответ 7
Вы также можете использовать функцию arrayfire shift в качестве замены Matlab circshift и повторно реализовать остальную часть кода. Это может быть полезно, если вы заинтересованы в любых других функциях AF в любом случае (например, переносимость на GPU, просто изменив флаг компоновщика).
Однако, если весь ваш код предназначен для запуска на CPU и достаточно сложный, или вы не хотите использовать какой-либо другой формат данных (AF требует af:: массивы), придерживайтесь одной из других опций.
В итоге я перешел на AF, потому что мне пришлось бы повторно реализовать fftshift как ядро OpenCL, если не было обратно.
Ответ 8
Octave использует fftw для реализации (i) fftshift.
Ответ 9
Это даст эквивалентный результат ifftshift в matlab
ifftshift(vector< vector <double> > Hlow,int RowLineSpace, int ColumnLineSpace)
{
int pivotRow=floor(RowLineSpace/2);
int pivotCol=floor(ColumnLineSpace/2);
for(int i=pivotRow;i<RowLineSpace;i++){
for(int j=0;j<ColumnLineSpace;j++){
double temp=Hlow.at(i).at(j);
second.push_back(temp);
}
ifftShiftRow.push_back(second);
second.clear();
}
for(int i=0;i<pivotRow;i++){
for(int j=0;j<ColumnLineSpace;j++){
double temp=Hlow.at(i).at(j);
first.push_back(temp);
}
ifftShiftRow.push_back(first);
first.clear();
}
double** arr = new double*[RowLineSpace];
for(int i = 0; i < RowLineSpace; ++i)
arr[i] = new double[ColumnLineSpace];
int i1=0,j1=0;
for(int j=pivotCol;j<ColumnLineSpace;j++){
for(int i=0;i<RowLineSpace;i++){
double temp2=ifftShiftRow.at(i).at(j);
arr[i1][j1]=temp2;
i1++;
}
j1++;
i1=0;
}
for(int j=0;j<pivotCol;j++){
for(int i=0;i<RowLineSpace;i++){
double temp1=ifftShiftRow.at(i).at(j);
arr[i1][j1]=temp1;
i1++;
}
j1++;
i1=0;
}
for(int i=0;i<RowLineSpace;i++){
for(int j=0;j<ColumnLineSpace;j++){
double value=arr[i][j];
temp.push_back(value);
}
ifftShiftLow.push_back(temp);
temp.clear();
}
return ifftShiftLow;
}
Ответ 10
Вы можете использовать kissfft. Это разумно быстро, чрезвычайно просто в использовании и бесплатно. Для организации вывода, как вы хотите, требуется только:
a) сдвиг на (-dim_x/2, -dim_y/2,...) с периодическими граничными условиями
b) FFT или IFFT
c) сдвиг назад на (dim_x/2, dim_y/2,...) с периодическими граничными условиями
d) масштаб? (согласно вашим потребностям IFFT * FFT будет масштабировать функцию по dim_x * dim_y *... по умолчанию)