1
$\begingroup$

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.

figure1

$\endgroup$
10
  • $\begingroup$ 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. $\endgroup$
    – Jason R
    May 16, 2014 at 16:32
  • $\begingroup$ I'm confused. You're windowing the spectrum? But you expect it to be unwindowed when inverse transformed? $\endgroup$
    – endolith
    May 16, 2014 at 16:42
  • 1
    $\begingroup$ @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. $\endgroup$
    – Jason R
    May 16, 2014 at 16:52
  • $\begingroup$ @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? $\endgroup$
    – endolith
    May 16, 2014 at 16:57
  • 1
    $\begingroup$ @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. $\endgroup$
    – wsdzbm
    May 18, 2014 at 8:11

2 Answers 2

1
$\begingroup$

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: enter image description here

$\endgroup$
4
  • $\begingroup$ bingo. perfect. $\endgroup$
    – wsdzbm
    May 18, 2014 at 17:39
  • $\begingroup$ Whether or not should I taper the ends of x before fft? $\endgroup$
    – wsdzbm
    May 19, 2014 at 0:34
  • 2
    $\begingroup$ Why would you need to do it? How do you think it will benefit you? $\endgroup$
    – learner
    May 20, 2014 at 3:37
  • $\begingroup$ 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. $\endgroup$
    – wsdzbm
    May 22, 2014 at 18:47
3
$\begingroup$

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));
$\endgroup$
1
  • $\begingroup$ y r correct. I was wrong. y1 was symmetrical but not conjugate. $\endgroup$
    – wsdzbm
    May 18, 2014 at 17:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.