% Generate input data: square wave built from the first few frequency/phase harmonics
fmax = 10; % maximum frequency, GHz
Nf = 201; % frequency points, from zero to fmax
data = zeros(1,Nf);
data(1) = 1; % square wave DC offset
cycles = 4; % square wave cycles
phase = 2.1; % square wave phase shift, radians
for k = 1:2:15 % sum several harmonics
data(1 + k * cycles) = 1.27 / k * exp(1j * ((k-1)/2 * pi + k * phase));
end
freq = 0:fmax/(Nf-1):fmax;
subplot(3,1,1); plot(freq, abs(data)); xlabel('Input Spectrum [GHz]');
% Calculate ifft of half-spectrum data
h = [conj(data) fliplr(data(2:length(data)-1))];
h(1) = 2 * real(h(1)); % DC tweak
h(length(data)) = 2 * real(h(length(data))); % Nyquist tweak
y = ifft(h);
% Plot results
N = length(h);
freq = 2 * fmax * (-N/2 : N/2-1) / N;
time = (0 : N-1) / 2 / fmax;
subplot(3,1,2); plot(freq, abs(fftshift(h))); xlabel('Expanded Spectrum [GHz]');
subplot(3,1,3); plot(time, N / 2 * y); xlabel('Time-Domain [nanoseconds]');