2

If I start off with a signal which has only real values, performing an fft and ifft returns the exact signal back with no complex entries as expected. But if I pad the fft with zeros to obtain interpolated values in the time domain, the inverse fft always turns out to be a complex double. I have taken care to perform the fftshift() and then pad on both sides, so that the symmetry is not broken. Following is an example code, that shows this behaviour. Am I looking at this the wrong way, or is the computational error after zero padding a bit too much? How do I overcome this?

Code:

%%%%%%%%%

x=linspace(0,2*pi,200);

y = sin(x)+sin(2*x);

Y = fftshift(fft(y));

n=400;

x1 = linspace(0,2*pi,n);

Y1 = zeros(1,n);

Y1((n-200)/2+1:end-(n-200)/2) = Y;

y2 = ifft(fftshift(2*Y1));

plot(x,y);

hold on;

plot(x1,y2(1:end),'x-');

isreal(y)==isreal(y2)

%%
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
Sachu17
  • 61
  • 5
  • Note: `Y2` is not defined, and why are you comparing at the end vs `Y2`, I guess it should be `y2`? – Irreducible Oct 16 '17 at 12:38
  • For comparison there is a FEX submission where you ca see how it is coded: [Frequency-domain-zero-padding](https://de.mathworks.com/matlabcentral/fileexchange/58854-frequency-domain-zero-padding-resampling--interpolation-) – Irreducible Oct 16 '17 at 12:44
  • my bad, it was a typo, it should be y2, but that is not the issue though – Sachu17 Oct 16 '17 at 13:16
  • @Irreducible, I checked the code on that page. It does exact same thing. Only that, at the end before returning, only real part is considered. But, ignoring the imaginary part, just changes the signal. Am I right? – Sachu17 Oct 16 '17 at 13:22

2 Answers2

2

Probably an off by one error, or not ignoring microscopic rounding-error values. For the result of an IFFT to be strictly real (except for tiny numerical noise quantities), the complex input vector has to be exactly conjugate symmetric around element 0 (or around N/2 for even N).

hotpaw2
  • 70,107
  • 14
  • 90
  • 153
0

There are several, small issues with your code:

  1. The number of samples of the original signal, and its FFT, should be odd. To see why, note that the first entry of the FFT is the zero-frequency component, which has no symmetric frequency. The number of remaining frequencies should be even, so that they can be split in two symmetric parts, and the additional zeros can be inserted between the two parts.

  2. The time axes should not be created with linspace. Since linspace fixes both endpoints, you don't get an exact scale factor of 2. To do that you need to create the axes manually with : (colon).

  3. The second fftshift should be ifftshift. This is not important here, because the zero-padded FFT has even size and so fftshift and ifftshift coincide. But as a general rule you should use ifftshift to undo the action of fftshift.

So, the code becomes:

x = (0:200)/201*2*pi;
y = sin(x)+sin(2*x);
Y = fftshift(fft(y));
x1 = (0:401)/402*2*pi;
Y1 = [zeros(1,101) Y zeros(1,100)];
y2 = ifft(ifftshift(2*Y1)); % or fftshift, since length is even
plot(x,y,'o-');
hold on;
plot(x1,y2,'.:');
isreal(y)==isreal(y2)
max(abs(y-y2(1:2:end))) % should be of the order of eps

Now we can check that

  • Both signals are real:

    >> isreal(y)==isreal(y2)
ans =
  logical
   1

    A better test would be to check that the imaginary component of y1 (which might not be exactly zero due to numerical precision issues) is of the order of eps.

  • Every other sample of y1 coincides with the corresponding sample in y, up to numerical errors of the order of eps:

    >> max(abs(y-y2(1:2:end)))
ans =
     6.661338147750939e-16

enter image description here

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • Just a minor correction; the length of the inverse FFT does not have to be odd for conjugate symmetry to result in a real output. The requirement is simply that Y(2:end) == conj(Y(end:-1:2)). This implies that entry n/2 is real. Take, for instance, `A = rand(3,1) + 1j*rand(3,1); Y = [10; A; 4; conj(flip(A))]; ifft(Y)`. The length here is even, but the result is real. – CKT Oct 18 '17 at 14:26
  • @CKT What I meant is that the length has to be odd _if we want to be able to pad with zeros_. If it is even, the component at normalized frequency pi rad/s would have to be "split in two" to accomodate the zeros in between – Luis Mendo Oct 18 '17 at 14:49