MATLAB: Time domain signal reconstruction from frequency domain

reconstruction fft time domain frequency domain conversion reconstruction

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.

Best Answer

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));
Related Question