time domain signal reconstruction from frequency domain

Hi everybody and thanks for taking time on this.

I have a time domain signal vectorized at a sampling frequency of 16kHz. I want to convert it back to time domain by using phase and magnitude, after I have made some modification to its spectrum. The code is, e.g.:

[x, fs]=wavread('file'); l_x=length(x);   
M=round((32*1e-3)*fs); N=M/2; % analysis window 
shift=floor((l_x-M)/N);
G=ones(1,length(R));                

for i=1:shift
    right=fft( x_r((i-1)*N+1:M+(i-1)*N) ); 

      phaseR=angle( right );      
      magniR(i,:) = abs(right(1:N)); %take only first M/2 points

      magniRR(i,:)=magniR(i,:).*G;        %apply spectrum modification

frame_r(i,:)= real(ifft( [magniRR(i,:) magniRR(i,N:-1:1)].*exp(1j*phaseR)' ))

end

% % overlap-and-add syntheses Allen & Rabiner's method
indf = N*[ 0:(shift-1) ]; inds = [ 1:M ]';       
indexes = indf(ones(M,1),:) + inds(:,ones(1,shift));
rr = zeros(1, l_x); wsumR = zeros(1, l_x); window=hann(M);

for m = 1:shift, 
    rr(indexes(:,m)) = rr(indexes(:,m))+ frame_r(m,:).*window';  
    wsumR(indexes(:,m)) = wsumR(indexes(:,m)) + window';
end

rr = rr./wsumR;% divide out summed-up analysis windows
plot(rr)
soundsc(rr,fs);

There is something wrong with this.

Answer 1

[Dr. Seis]

Why are you taking just half of the spectrum?

You need to apply the modification to the entire frequency range (i.e., both positive and negative frequencies). The "fft" needs the amplitudes from both sides of the frequency spectrum to correctly construct the signal in the time domain. The fact that the real parts of the amplitudes in the frequency domain are symmetric about 0 frequency and the imaginary parts of the amplitudes in the frequency domain are anti-symmetric about 0 frequency is what "forces" the data to be 100% real in the time domain. If you are missing half the spectrum, which makes for no symmetry, then your time domain product will be complex (i.e., have both real and imaginary parts).

Just so I understand what is going on, I usually construct and modify things in the frequency domain by taking the following into account:

g = rand(1,512); % define a random time series "g"
dt = 0.1; % define a time increment (seconds per sample)
N = length(g);

Nyquist = 1/(2*dt); % Define Nyquist frequency
df = 1 / (N*dt); % Define the frequency increment

G = fftshift(fft(g)); % Convert "g" into frequency domain and shift to frequency range below
f = -Nyquist : df : Nyquist-df; % define frequency range for "G"

You can check and confirm for yourself that the number of amplitudes in "g" are the same number of amplitudes in "G" - this is a property of the Fourier transform! Now you can do things to the amplitudes according to what frequency they represent. For example, you can take the derivative of the time signal by doing this in the frequency domain:

G_new = G.*(1j*2*pi*f);

Then convert back to the time domain by:

g_new = ifft(ifftshift(G_new));

[joeDiHare]

What I do is taking half the spectrum points, make some modification on the magnitude (e.g. apply a masking G) and then copy/paste the modified half magnitude into the negative frequencies as well, while phase is the original one. (it was not shown in the code, i added it now).

Besides, I found that using only 1 window in the OLA works better than using also a second window in that fft calculation.

Do you think these two steps are conceptually wrong?

(thanks very much)

[Dr. Seis]

The imaginary parts of the amplitudes for the positive frequencies are *not* a mirror image of the negative ones - they are flipped (i.e. a positive imaginary amplitude on the positive frequency side will be a negative imaginary amplitude on the negative frequency side). Only the real parts of the amplitudes are a mirror image.

Answer 2

[Wayne King]

It looks like you are implementing the overlap and add method. This method is an efficient way (if it is coded correctly) of implementing the convolution of a very long signal with an FIR filter by breaking the long convolution up into shorter segments, which is what you are calling frames above.

I'm not sure why you are using the polar form of the DFT coefficients and not just using ifft()

I think this is unnecessary:

frame_r(i,:)=abs(ifft(...
[RR(i,:) abs(right(N+1:end))'].*exp(1j*phaseR)' )).*...
cos(angle(ifft([RR(i,:) abs(right(N+1:end))'].*exp(1j*phaseR)' )) );

See:

http://en.wikipedia.org/wiki/Overlap-add_method

For a general introduction to overlap and add.

Wayne

[joeDiHare]

Hi Wayne, thanks for your answer.

You are rigth, polar form is unnecessary and I am now using:

x = real(ifft( magnitude.*exp(1j*phase) );

The reconstruction back from STFT is however still a bit artifacty; it sounds like it is "jigging" around the transictions of the OLA windows.

I will play a bit with different windows, and different overlaps.

Answer 3

[Wayne King]

You're still using the polar form. All you need to do:

x = randn(8,1);
xdft = fft(x);
x1 = ifft(xdft);
% what you're doing
x2 = ifft(abs(xdft).*exp(1j*phase(xdft)));

[joeDiHare]

I know, but as I stated, I want to modify the magnitude spectrum before reconstruct back, which is why I need the polar form. Certainly, previous formula I used is a bit unecessary (polar form with eurler decomposition of Im part).