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 = (rand(41,70,60,10));
   X=complex(X,X);

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

percentDiff = 3.6174e-14

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

t0 = 0.0045

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

t1 = 0.0860

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

t2 = 0.0107

t2/t0

ans = 2.4062

t2/t1

ans = 0.1245

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.

[Matt J]

This seems a bit better:

clc,rng(100);

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

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

percentDiff = 7.4344e-15

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

t0 = 0.0050

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

t1 = 0.0967

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

t2 = 0.0107

t2/t0

ans = 2.1254

t2/t1

ans = 0.1107

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

 if strcmp(symflag,'symmetric')

    X(:,1,:)=ifft(X(:,1,:),[],1,'symmetric');
      w=ceil(width(X)/2);
    X(:,2:w,:)=ifft(X(:,2:w,:),[],1);

   if rem(width(X),2)==0
     X(:,w+1,:)=ifft(X(:,w+1,:),[],1,'symmetric');
   end

    x=ifft(X,[],2,'symmetric');

 else
 
    x=ifft2(X);

 end

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]

Took me a while to figure out what's happening in ifft2nLoop.

I was unaware of that usage of size where the last sz output is returned as the product of that dimension and all trailing dimensions.

And that Y(:,:,i) thing .... I didn't realize that linear indexing can be used like that. Is that documented anywhere?

Much better approach than mine using cellfun.

[Matt J]

And that Y(:,:,i) thing .... I didn't realize that linear indexing can be used like that. Is that documented anywhere?

Generally speaking you can mix subscript and linear indexing as long as it's in the form,

Y(subscript_1,subscript_2,...,subscript_n, linearindex)

I haven't been able to find documentation of this, except of course for the familiar case where n=0.

[Matt J]

Tech Support confirms that it is a bug.

[Paul]

We've talked about two (related) issues in this thread.

(a) with 'nonsymmetric' the output of ifft2 for a 2D input is not exactly the same as the corresponding slice of the output when that same 2D input is a slice of an ND input to ifft2

(b) for an ND input, ifft2 doesn't respect 'symmetric' in both directions, which can lead to incorrect results when the input has a large break in symmetry, e.g., due to zero filling.

Which of those (or both) is the confirmed bug?

Did they say that they plan publish a public bug report?

[Matt J]

(a) with 'nonsymmetric' the output of ifft2 for a 2D input is not exactly the same as the corresponding slice of the output when that same 2D input is a slice of an ND input to ifft2

I assume you meant "with 'symmetric' "? What I reported were my findings from this earlier comment, and they agreed that the large percent error of 65%, seen when symflag='symmetric', was due to a bug. The only thing they said about follow-up action was "we expect it to be addressed in a future release".

(b) for an ND input, ifft2 doesn't respect 'symmetric' in both directions, which can lead to incorrect results when the input has a large break in symmetry, e.g., due to zero filling.

Not sure what (b) was, or when we concluded that zero filling produced unexpected discrepancies. When symflag='symmetric', it is to be expected that large breaks in symmetry will produce discrepancies, because the negative frequency axis data is being expressly ignored.

[Paul]

This has been a long thread and I need a reset.

percentError=@(a,b) norm(a(:)-b(:),inf)/norm(b(:),inf)*100;
rng(12345);

Here are the cases I'm tracking:

Case 1: symmetric transform, call with non-symmetric, 2D doesn't match ND slice, but they're close

x = rand(10);
X = fft2(x); % symmetric transform of real signal
x1 = ifft2(X,'nonsymmetric');
x2 = ifft2(cat(3,X,X),'nonsymmetric');
[percentError(x1,x2(:,:,1)),percentError(x,x1),percentError(x,x2(:,:,1))]

ans = 1×3

1.0e-13 * 0.2234 0.3351 0.2513

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

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

Case 2: symmetric transform, call with symmetric, 2D doesn't match ND slice, but they're close

x1 = ifft2(X,'symmetric');
x2 = ifft2(cat(3,X,X),'symmetric');
[percentError(x1,x2(:,:,1)),percentError(x,x1),percentError(x,x2(:,:,1))]

ans = 1×3

1.0e-13 * 0.2234 0.3351 0.2513

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

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

Case 3: nonsymmetric transform, call with nonsymmetric, 2D doesn't match ND slice, but they're close

x = rand(10) + 1j*rand(10);
X = fft2(x); % nonsymmetric transform of complex signal
x1 = ifft2(X,'nonsymmetric');
x2 = ifft2(cat(3,X,X),'nonsymmetric');
[percentError(x1,x2(:,:,1)),percentError(x,x1),percentError(x,x2(:,:,1))]

ans = 1×3

1.0e-13 * 0.3500 0.3797 0.2884

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

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

Case 4: nonsymmetric transform, call with symmetric, 2D doesn't match ND slice, and they're not close

x1 = ifft2(X,'symmetric');
x2 = ifft2(cat(3,X,X),'symmetric');
[percentError(x1,x2(:,:,1)),percentError(x,x1),percentError(x,x2(:,:,1))]

ans = 1×3

63.2666 81.1045 77.5687

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

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

●

Of course, neither result from Case 4 is anywhere close to the original input to fft2 because that input doesn't have a symmetric transform to begin with.

Does the confirmed bug cover all four cases? Or just Case 4? (Or some other subset of Cases?)

[Matt J]

The bug covers only cases where symflag='symmetric', i.e. Cases 2 and 4. The only reason you don't see significant errors in Case 2 is because you've given it symmetric input data, so the input flag symflag='symmetric' doesn't really have to do any work.

[Paul]

Why would Case 2 be considered a bug but not Cases 1 and 3?

Or maybe Cases 1 and 3 could be considered as a distinctly different issue than Case 2; the issue being that ifft2 doesn't return the same output for the same input being 2D or a slice of ND? Are you o.k. with the behavior of Cases 1 and 3?

BTW, I don't really have a dog in this hunt because, IIRC, I've never used fft2/ifft2 in anger, so I'm quite interested in other's opinions of those functions' functionality.

[Matt J]

So, just to be clear, I don't judge whether a bug is in play based solely on input/output. I'm judging it also based on how we know the data is processed in each of these cases. We haven't really been able to identify anything improper about how the processing is done when symflag='nonsymmetric' within ifft2.m. There are discrepancies in the outputs of the Case 1 and Case 3 tests, but those appear to arise from normal floating point precision noise.

Conversely, we have identified impropriety in the way symflag='symmetric' is processed when X is 3D or higher and Case 4 exposes that. I consider the bug to be in play in Case 2 as well, because the data is being processed in the same incorrect wasy as in Case 4. Granted, the test in Case 2 doesn't expose the bug - you are again getting agreement within float precision. But that's only because, with an already-symmetric X as input, it's an uninformative test. It doesn't mean the processing is being done properly.