N = 512; % Number of times, and of frequencies
minf = -2; % Minimum frequency (Nyquist=1)
maxf = 4; % Maximum frequency (Nyquist=1)
df = (maxf-minf)/N; % Frequency interval
window = hanning(N)*ones(1,N); % Apodization

freq = (0:N-1).*df + minf; % Frequency row vector (Nyquist=1)
t = (0:N-1)'; % Time column vector

freq_mat = ones(N,1) * freq;
   % N x N  array of frequencies, varying along the horizontal.
   
t_mat = t * ones(1,N);
   % N x N  array of times, varying along the vertical.
   
signals = exp(i * pi .* freq_mat .* t_mat);
   % N x N  array of signals, each column being a sinusoid
   % at a particular frequency given by freq_mat.

signals_f = fft( window .* signals ); % FFT of each column
signals_power = abs(signals_f).^2 ; % Power spectrum of each column

% Display the power spectra as an image.
hold off
imagesc(freq, (1:N), -signals_power); % minus sign inverts the image.
colormap(gray); % replaces the ugly default color table.
xlabel('signal frequency (Nyquist=1)','fontsize',20);
ylabel('frequency bin in sampled data','fontsize',20);

% Mark out the range of frequencies |f| < Nyquist
hold on
plot(ones(1,100), (1:(N-1)/99:N), 'k.')
plot(-1.*ones(1,100), (1:(N-1)/99:N), 'k.')