time domain signal reconstruction from frequency domain

48 次查看(过去 30 天)
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.

采纳的回答

Dr. Seis
Dr. Seis 2011-10-11
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));
  2 个评论
joeDiHare
joeDiHare 2011-10-11
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
Dr. Seis 2011-10-11
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.

请先登录,再进行评论。

更多回答(2 个)

Wayne King
Wayne King 2011-10-10
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:
For a general introduction to overlap and add.
Wayne
  1 个评论
joeDiHare
joeDiHare 2011-10-10
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.

请先登录,再进行评论。


Wayne King
Wayne King 2011-10-10
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)));
  1 个评论
joeDiHare
joeDiHare 2011-10-10
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).

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Filter Analysis 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by