Does anyone have an equivilent to the matlab/octave/sci.py 'fftshift' routine that is written in C? I am using FFTW directly to compute rank 1 and 2 FFTs but need to 'shift' the zero frequency component. -- FFTSHIFT Shift zero-frequency component to center of spectrum. For vectors, FFTSHIFT(X) swaps the left and right halves of X. For matrices, FFTSHIFT(X) swaps the first and third quadrants and the second and fourth quadrants. For N-D arrays, FFTSHIFT(X) swaps "half-spaces" of X along each dimension. -- Thanks in advance, Graham

# fftshift equivilent in C

Started by ●August 23, 2003

Reply by ●August 23, 20032003-08-23

> Does anyone have an equivilent to the matlab/octave/sci.py 'fftshift' > routine that is written in C? I am using FFTW directly to compute > rank 1 and 2 FFTs but need to 'shift' the zero frequency component. > > -- > FFTSHIFT > Shift zero-frequency component to center of spectrum. > For vectors, FFTSHIFT(X) swaps the left and right halves of > X. For matrices, FFTSHIFT(X) swaps the first and third > quadrants and the second and fourth quadrants. For N-D > arrays, FFTSHIFT(X) swaps "half-spaces" of X along each > dimension. > -- > > Thanks in advance, > GrahamGraham, to reorder in 1-D: int n; // FFT output vector length int n2; int i; complex x[n]; // you might typedef your own complex type complex tmp; n2 = n / 2; // half of vector length for (i = 0; i < n2; i++) { tmp = x[i]; x[i] = x[i+n2]; x[i+n2] = tmp; } to reorder in 2-D: int m, n; // FFT row and column dimensions might be different int m2, n2; int i, k; complex x[m][n]; complex tmp13, tmp24; m2 = m / 2; // half of row dimension n2 = n / 2; // half of column dimension // interchange entries in 4 quadrants, 1 <--> 3 and 2 <--> 4 for (i = 0; i < m2; i++) { for (k = 0; k < n2; k++) { tmp13 = x[i][k]; x[i][k] = x[i+m2][k+n2]; x[i+m2][k+n2] = tmp13; tmp24 = x[i+m2][k]; x[i+m2][k] = x[i][k+n2]; x[i][k+n2] = tmp24; } } Hope this helps, Andrew

Reply by ●August 23, 20032003-08-23

On 23 Aug 2003, Andrew Nesterov wrote:> int n; // FFT output vector length > int n2; > int i; > complex x[n]; // you might typedef your own complex type > complex tmp; > > n2 = n / 2; // half of vector length > > for (i = 0; i < n2; i++) > { > tmp = x[i]; > x[i] = x[i+n2]; > x[i+n2] = tmp; > }This is close to what I came up with but it only works for even n. Matlab gives:>> fftshift([1,2,3,4,5,6,7,8])ans = 5 6 7 8 1 2 3 4>> fftshift([1,2,3,4,5,6,7,8,9])ans = 6 7 8 9 1 2 3 4 5 Above gives: n = 8 n2 = 4 in = 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 out = 5.00 6.00 7.00 8.00 1.00 2.00 3.00 4.00 n = 9 n2 = 4 in = 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 out = 5.00 6.00 7.00 8.00 1.00 2.00 3.00 4.00 9.00 Any ideas on how to get the same answer as matlab for both even and odd n? Thanks, Graham

Reply by ●August 24, 20032003-08-24

Graham, It can be done out of place rather simply, but note a new array, z[n] was introduced. It also works for even n: for (i = 0; i < n; i++) { // find circulant index k = (i+n2) % n; // or, without modulo call k = (i+n2); if (k >= n) k = k - n; // move z[k] = x[i]; } For this particular case, an odd n and shift by [n/2], it might be done in place, but more complicated. As it is seen, the only problem is the overlap of i-th and (i+n2)-th entries: 0 1 2 3 4 5 6 7 8 5 6 7 8 0 1 2 3 4 ^ 0-th entry is moved into 4-th location, but 4-th location instead of been moved to 0-th position, as in an even length case, is moved into the last location, etc. This happens due to curcular addressing mode: e.g. n = 9, n2 = 4 and 0+4 mod 9 = 4, 4+4 mod 9 = 8, but 5+4 mod 9 = 0 and the 5th entry eventually goes into the first location. In an even length case, e.g. n = 8, the 4-th location is moved into 4+4 mod 8 = 0th location, thus the circulant shift might be accomplished by swap operations. The code works for even n as well (again, n2 = n/2): for (i = 0; i < n2; i++) // swap what can be swapped { tmp = x[i]; x[i] = x[i+n2]; x[i+n2] = tmp; } if (n & 1) // odd n, shift the rest { tmp = x[n-1]; x[n-1] = x[0]; for (i = 1; i < n2; i++) { x[i-1] = x[i]; } x[n2-1] = tmp; } It might be done in place more simply, with 4 load/store ops less, but works only for odd n. For even n, (ip1+n2) is n when i becomes n2-1, thus load/store is being done out of array bound on the last iteration. For odd n, (ip1+n2) is n-1 on the last iteration: for (i = 0; i < n2; i++) { ip1 = i + 1 tmp = x[i+n2]; x[i+n2] = x[i]; x[i] = x[ip1+n2]; x[ip1+n2] = tmp; } Regards, Andrew> > int n; // FFT output vector length > > int n2; > > int i; > > complex x[n]; // you might typedef your own complex type > > complex tmp; > > > > n2 = n / 2; // half of vector length > > > > for (i = 0; i < n2; i++) > > { > > tmp = x[i]; > > x[i] = x[i+n2]; > > x[i+n2] = tmp; > > } > > This is close to what I came up with but it only works for even n. > > Matlab gives: > > >> fftshift([1,2,3,4,5,6,7,8]) > ans = 5 6 7 8 1 2 3 4 > > >> fftshift([1,2,3,4,5,6,7,8,9]) > ans = 6 7 8 9 1 2 3 4 5 > > Above gives: > > n = 8 > n2 = 4 > in = 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 > out = 5.00 6.00 7.00 8.00 1.00 2.00 3.00 4.00 > > n = 9 > n2 = 4 > in = 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 > out = 5.00 6.00 7.00 8.00 1.00 2.00 3.00 4.00 9.00 > > Any ideas on how to get the same answer as matlab for both even and odd n? > > Thanks, > Graham