Question using IFFT, try to transform continuous 1 from frequency in time domain

Hey. I have a frequency response of a system. I want to see this frequency response in the time domain, so I want to transform it into this. Here for I want to create a program, where I enter a one-sided frequency spectrum, so the frequency vector (freq) with the associated amplitudes (y_f).

I have now created this program. If I want to transform a pulse (e.g. at 10 kHz, see 1. Example), the time axis and amplitude in the time domain look good.

However, if I try to transform a continuous 1 signal, the result does not look correct anymore. The time axis should be correct, but unfortunately the amplitude is not. As far as I know, the amplitude should be 1 at t = 0, but it has the value 150.000 (which is half of my vector length).

The problem wonders me, since 1. example seems to work.

Does anyone know where my problem is?

Thanks a lot!

Complete program follows:

clear

% signal input for one-sided spectrum

freq = 0:(10/3):500e3;      % preset frequency points

%%%%% 1. Example: dirac impuls at 10 kHz            (expceting cos with T = 100 us)
% for k = 1:length(freq)  % filling y_f with zeros
%     y_f(k) = 0;      
% end
% y_f(3001)=1;   % set sample 3001 (10 kHz) to 1
% figure(1)
% plot(freq,y_f)
%%%%%

%%%%% 2. Example: continuous 1                      (expecting dirac at t = 0)
for k = 1:length(freq)  % filling y_f with zeros
    y_f(k) = 1;      
end
figure(1)
plot(freq,y_f)
%%%%%


%%% IFFT
y_fQ1 = y_f;                 % 1. Quadrant
y_fQ2 = conj(y_f);           % 2. Quadrant Complex conjugate
y_fQ2 = fliplr(y_fQ2(2:end));% Flip array left to right


y_f2 = [y_fQ2,y_fQ1(1)*1 ,y_fQ1(2:end)];   % two-sided frequency 
y_f2 = y_f2 * 0.5 * length(y_f2);          % Adjust amplitude for ifft
y_f2 = ifftshift(y_f2);                    % shift spectrum for ifft    ...i.I.i...   -->   I.i......i.

y_t = ifft(y_f2,'symmetric');

delta_f = freq(end)-freq(end-1);        % Frequency resolution 
L = 2*length(freq)-1;


Fs = freq(end)*2;        % Sampling frequency

delta_t = 1/Fs;          % time resolution

t =  0 : delta_t : (L-1)*delta_t;   %% create time vector

figure(3)
plot(t, y_t)    % plot result

Answer 1

[David Goodmanson]

Hi Jan,

In the secnd example the dirac value should be quite large, In fact it should be 3e5, the number of points, not half of that. In the style that you are doing things, y(-f) is created from y*(+f) and the entire array is [y(-f) y(0) y(+f)]. So, since

cos(f0*t) = (1/2) (exp(if0*t) + exp(-if0*t))

the amplitude at 10kHz should be 1/2, not 1. Then for scaling, the ifft is multiplied by N, the length of the array, not half the length. In the first example you get a time domain cosine of amplitude 1 like before. In the second example you end up with 3e5 (actually 300001) for the dirac peak height.

The changes are made in the code below. Also, the arrays are created without for loops. You can save time and space by doing those in Matlab fashion instead of using for loops.

As to why the peak is 3e5 instead of 1, one can look at energy conservation (Parseval's theorem). For first example, the cosine wave, there are peaks of height (1/2) at 10000 Hz and -10000 Hz, zero everywhere else. The total frequency domain energy is the sum of squares, (1/2)^2 + (1/2)^2 = (1/2). With the scaling you have chosen (ifft x N) the result is a cosine wave in the time domain of amplitude 1. Since the mean of cos^2 is (1/2), this example shows that in general the equivalent energy in the time domain is going to be the mean square value of the time domain signal.

For the second example, all ones in the frequency domain, the energy is much larger than in the first case: N*(1)^2 = N. In the time domain with its single peak of height N, the mean square value is N^2 / N = N. So it works.

clear
% signal input for one-sided spectrum
freq = 0:(10/3):500e3;                 % preset frequency points

% 1. Example: dirac impuls at 10 kHz   (expceting cos with T = 100 us)
y_f = zeros(size(freq));
y_f(3001) = 1/2;                       % set sample 3001 (10 kHz) to 1
figure(1)
plot(freq,y_f)

% 2. Example: continuous 1             (expecting dirac at t = 0)
% y_f = ones(size(freq));
% figure(1)
% plot(freq,y_f)

% IFFT 
y_fQ1 = y_f;                               % 1. Quadrant
y_fQ2 = conj(y_f);                         % 2. Quadrant Complex conjugate
y_fQ2 = fliplr(y_fQ2(2:end));              % Flip array left to right
y_f2 = [y_fQ2,y_fQ1(1)*1 ,y_fQ1(2:end)];   % two-sided frequency 
y_f2 = y_f2*length(y_f2);                  % Adjust amplitude for ifft         CHANGED
y_f2 = ifftshift(y_f2);                    % shift spectrum for ifft ..I.. -_>  I....
y_t = ifft(y_f2,'symmetric');
delta_f = freq(end)-freq(end-1);           % Frequency resolution 
L = 2*length(freq)-1;
Fs = freq(end)*2;                          % Sampling frequency
delta_t = 1/Fs;                            % time resolution
t =  0 : delta_t : (L-1)*delta_t;          % create time vector
figure(3)
plot(t, y_t)                               % plot result

[Jan Lindenau]

Thank you very much for your help. I answer unfortunately so late, because I wanted to insert this knowledge into my program and see if it works. It does! Thanks again for the explanation.