effects of ifft2(x, 'symmetric') while x is not conjugate symmetric

Hi all,

I understand that 'symmetric' in ifft2 assumes data is conjugate symmetric and will output real values, but I want to know how exactly is this done. Specifically:

1. with 'symmetric', will only half of the data be used or still the whole matrix is used, but just output real values?
2. if x is not conjugate symmetric but I still pass in 'symmetric', what will happen?

Thanks in advance.

Answer 1

[Matt J]

When you specify "symmetric", the upper half of the fft input array is ignored, e.g.

x=ones(1,7);

ifft(x)

ans = 1×7

1.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x(5:end)=rand(1,3); %write random junk into the upper half of x

ifft(x)

ans = 
   0.8258 + 0.0000i   0.0750 + 0.0982i  -0.0516 - 0.0461i   0.0637 + 0.1255i   0.0637 - 0.1255i  -0.0516 + 0.0461i   0.0750 - 0.0982i

ifft(x,'symmetric')

ans = 1×7

1.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Justin]

Thanks for your quick reply, that is the case!

Just wanted to add a little more details for 2D case, only the first (N/2+1) cols data, i.e., x(:, 1:N/2+1), will be considered with 'symmetric'

[Matt J]

Only the first M columns and also only the first M rows where,

M=floor((N+1)/2)

[Paul]

I don't think that's true for the rows that aren't in the first column.

X = ifft2(rand(10));
x1 = ifft2(X,'symmetric');
X(end,3) = 0;
x2 = ifft2(X,'symmetric');
isequal(x1,x2)

ans = logical
   0

[Matt J]

Yep, my bad. In higher dimensions, conjugate symmetry means the image is symmetric across the origin. So, if you view the array after fftshift(), the cells it would be using are the right half, minus half of one of the DC axes. In other words, it needs cells from two quadrants to cover the whole array by symmetry, and also only the non-negative DC axes are needed.

If you re-organize this in terms of the original ordering of X(i,j) it means the cells that can be discarded are the light-shaded ones below, which is consistent with what @Paul found in his tests:

Answer 2

[Paul]

It appears that only portions of the input matrix are used when 'symmetric' is specified on input to ifft2 and the input is 2D

rng(100);
X = rand(10,12) + 1j*rand(10,12);  % non-symmetric input
x1 = ifft2(X,'symmetric');
X(:,8:12) = 0;  % zero out the entire right side
X(7:10,1) = 0;  % zero out the top half of the first column
x2 = ifft2(X,'symmetric');
isequal(x1,x2)

ans = logical
   1

[Paul]

See this comment for visualization of why the zero'd out cells aren't needed.

[Paul]

ifft2 uses different operations depending on whether or not the input is a 2D matrix, i.e., ismatrix(X) is true. If true, then ifft2 calls ifftn, otherwise it uses two calls to ifft.

I was wondering if the conclusions in this thread apply for the N-D case as well, and in so checking came across an interesting behavior.

Case 1: the transform in the frequency domain is conjugate symmetric

rng(100)
X1 = fft2(rand(10,12));
x1 = ifft2(X1,'symmetric');  % garden variety 2D case, uses ifftn
X2 = cat(3,X1,X1);
x2 = ifft2(X2,'symmetric');  % ND case, uses ifft twice
isequal(x1,x2(:,:,1))

ans = logical
   0

In this case, where X1 is supposed to be conjugate symmetric by construction, the difference between x1 and x2 is small

norm(x1-x2(:,:,1),'fro')

ans = 1.2937e-15

But the results are not identical as might be expected.

Case 2: the transform in the frequency domain is not conjugate symmetric, but we treat is as such, as in the original example in this Answer.

X1 = rand(10,12) + 1j*rand(10,12);
x1 = ifft2(X1,'symmetric');
X2 = cat(3,X1,X1);
x2 = ifft2(X2,'symmetric');
isequal(x1,x2(:,:,1))

ans = logical
   0

norm(x1-x2(:,:,1),'fro')

ans = 0.6526

Now the difference is substantial.

Case 3 - Same as Case 1, zero out the "unneeded" terms.

rng(100)
X1 = fft2(rand(10,12));
x1 = ifft2(X1,'symmetric');  % garden variety 2D case, uses ifftn
X2 = cat(3,X1,X1);
X2(:,8:12,:) = 0;  % zero out the entire right side
X2(7:10,1,:) = 0;  % zero out the top half of the first column
x2 = ifft2(X2,'symmetric');  % ND case, uses ifft twice
isequal(x1,x2(:,:,1))

ans = logical
   0

norm(x1-x2(:,:,1),'fro')

ans = 2.2737

Unsurprsingly, that doesn't work.

I believe that these results follow in that ifft2 can only apply the 'symmetric' flag in one direction when using the two-call sequence to ifft for N-D inputs.

@Matt J

[Matt J]

depending on whether or not the input is a 2D matrix, i.e., ismatrix(X) is true. If true, then ifft2 calls ifftn, otherwise it uses two calls to ifft.

I'd call that a bug. While ifftn(z) is a separable operation (i.e., if can be decomposed into 1D ffts), ifftn(z,'symmetric') is not, and if treated as such will obviously give the wrong result for multi-slice input. Here is a further test that that is what's happening:

n=3;

X = rand(10,12,n) + 1j*rand(10,12,n);

test3D(X)

PercentError = 1.7087e-14

test3D(X,'symmetric')

PercentError = 65.4130

function test3D(X,varargin)
%Tests ifft2 on a 3D input array X

    n=size(X,3);
    
    xcell=cell(1,1,n);
    for i=1:n %loop over 3D slices, apply ifft2 to each slice
        xcell{i}=ifft2(X(:,:,i),varargin{:});
    end
    x_expected=cell2mat(xcell); %expected result
    
    x_observed=ifft2(X,varargin{:});  %observed result
    
    percentError=@(a,b) norm(a(:)-b(:),inf)/norm(b(:),inf)*100;
    
    PercentError=percentError(x_observed,x_expected)

end

[Paul]

Hi Matt,

I suspect that the developers were expecting that the 'symmetric' flag would be used on inputs that are not symmetric only due to accumulation of numerical trash in developing the input to ifft2, and not used on inputs that are fundamentally nonsymmetric as in your example and in my Case 2 and 3 examples.

In case you didn't see already, ifft2 is an ordinary m-file:

which ifft2

/MATLAB/toolbox/matlab/datafun/ifft2.m

If you do file a bug report, would you please post back here with a summary of the response?

[Matt J]

I have filed it, and will post back.

[Matt J]

I suspect that the developers were expecting that the 'symmetric' flag would be used on inputs that are not symmetric only due to accumulation of numerical trash.

Clearly they didn't assume that when ismatrix(X)=true. In that case, even large breaks in conjugate symmetry are not felt in the result.

 X1=rand(5);X1=X1+flipud(fliplr(X1)); 
 X2=X1; X2(1:12)=0; 
 X1=ifftshift(X1), X2=ifftshift(X2)

X1 = 5×5

1.2534 1.1076 1.3020 1.3020 1.1076 0.9483 0.7162 0.4986 0.6450 1.0906 1.3957 0.7819 1.3570 1.2072 1.4175 1.3957 1.4175 1.2072 1.3570 0.7819 0.9483 1.0906 0.6450 0.4986 0.7162

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X2 = 5×5

1.2534 1.1076 1.3020 0 0 0.9483 0.7162 0.4986 0 0 1.3957 0.7819 1.3570 0 0 0 1.4175 1.2072 0 0 0 1.0906 0.6450 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

  ifft2(X1,'sym')==ifft2(X2,'sym')

ans = 5×5 logical array
   1   1   1   1   1
   1   1   1   1   1
   1   1   1   1   1
   1   1   1   1   1
   1   1   1   1   1

[Paul]

I tried my hand at a basic function that computes the the 2D IFFT with ifftn on each slice of an N-D array

rng(100);
x = rand(24,23,32,25);
X = fft2(x);  

x1 = newifft2(X);  
norm(x-x1,'fro')

ans = 1.1307e-13

x2 = newifft2(X,'symmetric');
norm(x-x2,'fro')

ans = 1.1279e-13

t1 = timeit(@() ifft2(X))     % 2x-ifft on N-D array

t1 = 0.0013

t2 = timeit(@() newifft2(X))  % ifftn on each 2D slice

t2 = 0.0194

t2/t1

ans = 15.1649

t1 = timeit(@() ifft2(X,'symmetric'))     % 2x-ifft on N-D array

t1 = 0.0012

t2 = timeit(@() newifft2(X,'symmetric'))  % ifftn on each 2D slice

t2 = 0.0211

t2/t1

ans = 17.4536

function x = newifft2(X,symflag)
if nargin == 1
    symflag = 'nonsymmetric';
end
d = size(X);
indices = arrayfun(@(d) {ones(1,d)},d);
C = mat2cell(X,d(1),d(2),indices{3:end});
x = cellfun(@(C) ifftn(C,symflag),C,'uni',false);
x = cell2mat(x);
end

[Matt J]

Here's another attempt. It seems to improve upon the loop, if nothing else.

clc,rng(100);

X = fft2(rand(41,70,60,10));

x1 = ifft2nLoop(X,'symmetric');
x2=  ifft2n(X,'symmetric');
percentDiff = norm(x1-x2,'fro')/norm(x1,'fro')*100

percentDiff = 3.5485e-14

t0 = timeit(@() ifft2(X,'symmetric') ) %wrong result

t0 = 0.0049

t1 = timeit(@() ifft2nLoop(X,'symmetric'))  % ifftn on each 2D slice

t1 = 0.0940

t2 = timeit(@() ifft2n(X,'symmetric') )  %proposed

t2 = 0.0155

t2/t0

ans = 3.1307

t2/t1

ans = 0.1645

function x = ifft2n(X,symflag)
arguments
   X
  symflag='nonsymmetric';
end

 if strcmp(symflag,'symmetric')

   xsiz=size(X);
   [m,n,p]=size(X);
   mn=m*n; 
 
   I=reshape(1:mn,m,n);

   Q=round(fft2(ifft2(I,'symmetric')));
  
   msk=(Q~=I);

  X=reshape(X,mn,p);
  X=X(Q,:);
  X(msk,:)=conj(X(msk,:));
  X(1,:)=real(X(1,:));

  X=reshape(X,xsiz);
  
 end
 
   x=ifft2( X,symflag);

end



function Y=ifft2nLoop(X,varargin)
%Computes intended ifft2 on a n-D input array X by looping


    [~,~,n]=size(X);
    Y=X;
    for i=1:n %loop over 3D slices, apply ifft2 to each slice
        Y(:,:,i)=ifft2(X(:,:,i),varargin{:});
    end

end

[Paul]

Referring back to this comment:

"Clearly they didn't assume that when ismatrix(X)=true. In that case, even large breaks in conjugate symmetry are not felt in the result."

I think that is basically a bonus. In other words, the developers' (hypothesized) intent is to use 'symmetric' when the input is not symmetric only due to numerical trash. The fact that the ifft2 happens to yield the expected result with symflag 'symmetric' when ismatrix(X) = true even when X is far from symmetric is not in conflict with the (hypothesized) intent.

From ifft2:symflag : "When Y is not exactly conjugate symmetric due to round-off error, ...."

In any case, will be ineresting hear their reponse to your bug report.

[Matt J]

I don't have a response from the yet, probably due to Thanksgiving.

Regardless of developer intent, though, the documentation says we are supposed to get the same result from ifft2(X) as we would from ifft2nLoop(X), when X is n-D. So, at the very least, I think it has to be a documentation bug.