Converting frequency domain to Time Domain using IFFT

Hello Experts,

This problem confused me for a long time.

I have a function in the frequency domain k(w). I want to discrete it and use ifft to get correspond time series. I already have time_array and dt, that is [0:dt:length(t_array)*dt-dt]. how can i take w's sample rate, length and use ifft(do i need to use ifftshift) to let the result to correspond to time_array?

Thanks a lot and deeply gratitute!!

Answer 1

[Paul]

Hi luo xj,

Would need more information on what you have.

Here is an example showing how to reconstruct a signal. Maybe it can be adapted to your needs.

Define a fucntion and its Fourier transform.

syms t w real
k(t) = exp(-abs(t-1))*rectangularPulse(-3,3,t);
K(w) = simplify(fourier(k(t),t,w))

K(w) = 

Define the parameter for the frequency domain sampling and the time domain reconstruction

N = 100;
n = -N/2:N/2;
dt = 0.1;
tvals = n*dt;
wvals = n/N*2*pi/dt;

Sample the Fourier Transform (analogous to sampling k(w) in the problem statement

Kvals = double(K(wvals));

Reconstruct samples of the underslying time domain signal

kvals = fftshift(ifft(ifftshift(Kvals/dt)));

Compare the reconstructed samples to the original signal

figure
fplot(subplot(211),real(k(t)));
hold on
stem(tvals,real(kvals))
fplot(subplot(212),imag(k(t)));
hold on
stem(tvals,imag(kvals))
xlabel('Time (sec)')

[luo xj]

Omg!!! Can't believe u could answer this question! Thanks a lot!

I just don't know how to figure it out since I'm a newer. I try what u tips above but can't get great result, coulu u show more tips? Here is my code.

The data i've got:

dt = 1e-08;
N = 943;
time_array = 0:dt:N*dt-dt;

Now I got my frequency func:

syms w real
k(w) = w + 1i*0.5*sqrt((-1i*w))

I want to sample on k(w) and use ifft to get array correspond to time_array.

Here is my code:

fs = 1/dt;
df = fs/N;
omega = df.*(0:N-1);
k_omega = omega + 1i*(1/2)*(-1i.*omega).^(1/2);
k_t = ifft(k_omega);

How ever, the result k_t isn't correspond to time_array.

This question confused me for a long time, could u give more tips? Soooo gratitude Paul~

[luo xj]

I use my way on ur code but got incorrect result lolll, could u erase my confuse plz!

syms t w real
k(t) = exp(-abs(t-1))*rectangularPulse(-3,3,t);
K(w) = simplify(fourier(k(t),t,w));
N = 100;
n = -N/2:N/2;
dt = 0.1;
tvals = n*dt;
wvals = n/N*2*pi/dt;
down = wvals.^2+1;
up1 = -exp(-4-3.*wvals*1i);
up2 = exp(6.*wvals*1i)+exp(2)-2*exp(4+2.*wvals*1i)*1i-wvals.*exp(2)*1i;
Kvals = up1.*up2./down;
Kvals = double(Kvals);
%Kvals = double(K(wvals));
kvals = fftshift(ifft(ifftshift(Kvals/dt)));
figure
fplot(subplot(211),real(k(t)));
hold on
stem(tvals,real(kvals))
fplot(subplot(212),imag(k(t)));
hold on
stem(tvals,imag(kvals))
xlabel('Time (sec)')

[Paul]

As to the second comment, up2 was missing a term and one term was missing a factor of 1i.

syms t w real
k(t) = exp(-abs(t-1))*rectangularPulse(-3,3,t);
K(w) = simplify(fourier(k(t),t,w));
N = 100;
n = -N/2:N/2;
dt = 0.1;
tvals = n*dt;
wvals = n/N*2*pi/dt;
K(w)

ans = 

down = wvals.^2+1;
up1 = -exp(-4-3.*wvals*1i);
%up2 = exp(6.*wvals*1i)+exp(2)-2*exp(4+2.*wvals*1i)*1i                          -wvals.*exp(2)*1i;
up2 =  exp(6.*wvals*1i)+exp(2)-2*exp(4+2.*wvals*1i)   +wvals.*exp(6*wvals*1i)*1i-wvals.*exp(2)*1i;

Kvals = up1.*up2./down;
Kvals = double(Kvals);
%Kvals = double(K(wvals));
kvals = fftshift(ifft(ifftshift(Kvals/dt)));
figure
fplot(subplot(211),real(k(t)));
hold on
stem(tvals,real(kvals))
fplot(subplot(212),imag(k(t)));
hold on
stem(tvals,imag(kvals))
xlabel('Time (sec)')

[Paul]

As to the first comment

syms w real
k(w) = w + 1i*0.5*sqrt((-1i*w))

k(w) = 

The code assumed w is in Hz. Is that correct? Usually that symbol means frequency in rad/sec.

Assuming that k(w) is a valid Fourier Transform to begin with, the underlying time domain signal has infinite energy. The inverse transform of the first term on the right hand side is the derivative of the Dirac delta function. AFAICT, the second term does not have a closed form inverse transform

syms t real
K(t) = ifourier(k(w),w,t)  % assuming w in rad/sec

K(t) = 

Because we are dealing with an infinite energy signal, I don't think ifft can be applied. If we insist on using ifft for this problem, we are essentially applying a window to k(w) and losing all of the high frequency content. Frankly, I'm not even sure if a limiting argument, i.e., the more points we use to extend the sampling of k(w) the closer ifft will reconstruct K(t), can be applied in this case. I mean, we can always apply the code, but that doesn't mean it yields the correct answer. With that said ...

dt = 1e-08;
N = 943;
n = -N/2:N/2;
tvals = n*dt;
wvals = n/N*2*pi/dt;  % assuming argument of k(w) is rad/sec
kvals = wvals + 1i*0.5*sqrt((-1i*wvals));
Kvals = fftshift(ifft(ifftshift(kvals/dt)));

figure
plot(tvals,abs(Kvals))
xlabel('Time (sec)')
ylabel('K(t)')

[luo xj]

It was my carelessness that led to the second low-level mistake. I feel very ashamed of it.

As for first comment, w is rad/sec exactly (It's formula of power law). But as you can see, the time_array after you used ifft is [-N/2*dt: N/2*dt], but what I want time_array is [0,N*dt], maybe i can sample 2N points to solve this problem......

Thank u so much to answer question above, you guided me on the right path!

As for your consideration, em..., maybe I can't get ideal result by ifft, but I need to try in this way caused I don't have any other options lol.

[luo xj]

I'm sorry for giving the wrong example, I carelessly check my code. Here is my working code.

It couldn't run because the data like kgrid is generate by toolbox k-waves. U can see that I take many many stupid attempt....

So your great answer helped me a lot！

%sensor_data: num_sensor_points*N,means nums sensor points get pressure
%message in N*dt.
%init_data = sensor_data;
t_array = kgrid.t_array;
N = length(t_array);
dt = kgrid.dt;
%T = kgrid.Nt*dt;
fs = 1/dt;
df = fs/N;

% get k_omega
alpha = 1/2;
gamma = 1/2;
%omega = linspace(0,n*f-f,n);
%omega = 2*pi/(N*dt)*(0:N-1);
n = -N:N;
omega = n*df*pi;%2*pi*
%omega = fftshift(omega);
%omega = (0:N-1)*df;
%w = pi/(2*size(p,1)*dt) * (-size(p,1)/2:size(p,1)/2-1);

k_omega = omega + 1i*alpha*(-1i.*omega).^gamma;
tau = t_array;
%k_omega = fftshift(k_omega);
%k_omega_tau = exp(1i.*kron(k_omega(:),tau(:).'));
k_omega_tau = zeros(2*N+1,N);
for i = 1:2*N+1
    for j = 1:N
        k_omega_tau(i,j) = exp(1i*k_omega(i)*tau(j));
    end
end

% get k_t_tau
k_t_tau_old = zeros(2*N+1,N);
k_t_tau = zeros(N,N);
for i = 1:N
    k_t_tau_old(:,i) = fftshift(ifft(ifftshift(k_omega_tau(:,i)/dt)));
    k_t_tau(:,i) = k_t_tau_old(N+1:2*N,i);
end

%for i = 1:N
%    for j = i+1:N
%        k_t_tau(i,j) = 0;
%    end
%end


%k_t_tau = ifft(fftshift(k_omega_tau,1),n,1);
%k_t_tau = (ifft(k_omega_tau,n,1));
%test = ifft(k_omega_tau(:,100),n);
%k_t_tau = zeros(length_t,length_t);
%for i = 1:length_t
%    k_t_tau(i,:) = ifftshift(ifft(k_omega_tau(i,:)));
%end
k_t_tau = real(k_t_tau);
%k_t_tau_old = k_t_tau;
%for i = 1:N
%k_t_tau(:,i) = k_t_tau_old(:,N-i+1);
%end

% get attenuated data
nums_points = num_sensor_points;

%q_a = cumsum(sensor_data,2)*dt;
q_a = zeros(num_sensor_points,N);
q_a(:,1) = sensor_data(:,1);
for i = 2:N
    q_a(:,i) = q_a(:,i-1)+sensor_data(:,i);
end
q_a = q_a.*dt;

%q_a(:,1) = sensor_data(:,1); 
q = k_t_tau* q_a';
q = q'/sqrt(2*pi)*dt;
%test = k_t_tau;
%q = q_a*B';
p = zeros(num_sensor_points,N);
%p = [zeros(nums_points,1), diff(q,1,2)/dt];
p(:,1) = q(:,1)/dt;
for i = 2:N
    p(:,i) = (q(:,i)-q(:,i-1))/dt;
end
%atten_data = real(p)*1/sqrt(2*pi);
atten_data = real(p)*1;