Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

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

share|improve this question
    
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 at 16:32
    
I'm confused. You're windowing the spectrum? But you expect it to be unwindowed when inverse transformed? –  endolith May 16 at 16:42
1  
@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 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 at 16:57
1  
@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 at 8:11

2 Answers 2

up vote 1 down vote accepted

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

share|improve this answer
    
bingo. perfect. –  Lee May 18 at 17:39
    
Whether or not should I taper the ends of x before fft? –  Lee May 19 at 0:34
2  
Why would you need to do it? How do you think it will benefit you? –  learner May 20 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 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));
share|improve this answer
    
y r correct. I was wrong. y1 was symmetrical but not conjugate. –  Lee May 18 at 17:35

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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