time domain signal reconstruction from frequency domain
60 次查看(过去 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.
0 个评论
采纳的回答
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 个评论
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
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
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)));
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Filter Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!