fft(x) different from fft(ifft(fft(x))?

I normalized the spectrum of a time series and windowed the spectrum, then something strange happened. The code below can run in MATLAB or Octave.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
k = n/10;
w = sin(linspace(0,pi/2,k));
y1(1:k) = y1(1:k) .* w;
y1(m:-1:m-k+1) = y1(m:-1:m-k+1) .* w;
y1(n-m+3:n) = y1(m-1:-1:2);
x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');


As shown in subplot 3, the red line is the difference between x1 and x2, meaning that they are exactly the same, but as shown in subplot 4, their spectrum (black and red, respectively) is different. So, It's weird FFT of IFFT of a spectrum is different from the spectrum itself.

• Check your gains in fft() and ifft(). If I remember correctly, running through a fft()/ifft() round trip in MATLAB gives you a net gain of N, the FFT size. – Jason R May 16 '14 at 16:32
• I'm confused. You're windowing the spectrum? But you expect it to be unwindowed when inverse transformed? – endolith May 16 '14 at 16:42
• @endolith: I agree that what he's trying to do is not very clear, but if you trace through the above code, you effectively have at one point y2 = fft(real(ifft(y1));, so he expects y1 == y2, which is the crux of his question from my reading of it. – Jason R May 16 '14 at 16:52
• @JasonR: Why does it need real()? Is the windowing of the spectrum not correctly symmetrical, which produces small imaginary parts, which when removed cause y1 != y2? – endolith May 16 '14 at 16:57
• @endolith, I found the bug. Changing "y1(n-m+3:n) = y1(m-1:-1:2);" to "y1(n-m+3:n) = conj(y1(m-1:-1:2));" fixes it. – Lee May 18 '14 at 8:11

Hilmar was right, its a problem with the FFT indexing. Actually, you took care of the DC which does not appear in the right-side copy of FFT. Like-wise you've to take care of Fs/2 component.

Here is the modified working code.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
k = n/10;
w = sin(linspace(0,pi/2,k));
y1(1:k) = y1(1:k) .* w;
y1(m-k:m-2) = y1(m-k:m-2).*fliplr(w(2:end));
y1(m-1:m+k-2) = y1(m-1:m+k-2).*w;
y1(n-k+2:n) = y1(n-k+2:n).*fliplr(w(2:end));

x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');


The resulting figures are:

• bingo. perfect. – Lee May 18 '14 at 17:39
• Whether or not should I taper the ends of x before fft? – Lee May 19 '14 at 0:34
• Why would you need to do it? How do you think it will benefit you? – learner May 20 '14 at 3:37
• in fact, I didn't do the tapering. but other people told me it's better to taper the two ends and he didn't tell me why. – Lee May 22 '14 at 18:47

The problem appears to be that your windowing function doesn't preserve complex conjugate symmetry. So ifft(y1) has a significant imaginary part. By discarding this through the real() operation you induce a significant error which results in the discrepancy. To verify try

z = ifft(y1); plot(imag(z));

• y r correct. I was wrong. y1 was symmetrical but not conjugate. – Lee May 18 '14 at 17:35